Insights into the Interaction Mechanisms of Peptide and Non-Peptide Inhibitors with MDM2 Using Gaussian-Accelerated Molecular Dynamics Simulations and Deep Learning

Inhibiting MDM2-p53 interaction is considered an efficient mode of cancer treatment. In our current study, Gaussian-accelerated molecular dynamics (GaMD), deep learning (DL), and binding free energy calculations were combined together to probe the binding mechanism of non-peptide inhibitors K23 and 0Y7 and peptide ones PDI6W and PDI to MDM2. The GaMD trajectory-based DL approach successfully identified significant functional domains, predominantly located at the helixes α2 and α2’, as well as the β-strands and loops between α2 and α2’. The post-processing analysis of the GaMD simulations indicated that inhibitor binding highly influences the structural flexibility and collective motions of MDM2. Calculations of molecular mechanics–generalized Born surface area (MM-GBSA) and solvated interaction energy (SIE) not only suggest that the ranking of the calculated binding free energies is in agreement with that of the experimental results, but also verify that van der Walls interactions are the primary forces responsible for inhibitor–MDM2 binding. Our findings also indicate that peptide inhibitors yield more interaction contacts with MDM2 compared to non-peptide inhibitors. Principal component analysis (PCA) and free energy landscape (FEL) analysis indicated that the piperidinone inhibitor 0Y7 shows the most pronounced impact on the free energy profiles of MDM2, with the piperidinone inhibitor demonstrating higher fluctuation amplitudes along primary eigenvectors. The hot spots of MDM2 revealed by residue-based free energy estimation provide target sites for drug design toward MDM2. This study is expected to provide useful theoretical aid for the development of selective inhibitors of MDM2 family members.


Introduction
The tumor suppressor protein p53 plays a crucial role in regulating the cell cycle, promoting apoptosis, and repairing DNA damage, thereby safeguarding cells against malignant transformation [1,2].The active form of p53 is highly effective in suppressing the development of tumors.Notably, approximately 50% of all human cancers exhibit malfunctions in p53 function due to deletions or mutations in its DNA-binding domain [3].A host of strategies targeting the p53 pathway have emerged due to its inhibitory effects against tumors [4,5].Numerous proteins play critical roles in regulating the function of p53, such as MDM2 [6][7][8][9][10] and MDMX [11][12][13].In fact, inactivated p53 fails to yield appropriate responses to external stimuli, substantially increasing the risk of tumorigenesis.The overexpression of MDM2/MDMX and p53-MDM2 interactions leads to the inactivation of p53, controlled by negative feedback mechanisms or p53-MDM2.Thus, it is of high significance to probe effective approaches to hold back the MDM2-p53 interaction.
Molecules 2024, 29, x FOR PEER REVIEW 2 of 23 interconnected by β-strands and loops that dictate substrate specificity [14,15].The secondary structures of MDM2 bound to the inhibitor and the binding pocket are illustrated in Figure 1A,B, respectively.The hydrophobic sidechains of residues Phe19', Trp23', and Leu26' from p53 play a crucial role in mediating the interaction between p53 and MDM2 [16][17][18].These residues directly disrupt the binding of p53 to MDM2, making them a promising target for potential anticancer therapies.Various drug candidates, including small-molecule inhibitors and peptide inhibitors, have been developed to target the p53-MDM2 interaction.Peptide inhibitors that disrupt the interactions between p53 and its negative regulator MDM2 have the potential to activate p53 [19][20][21][22][23]. Additionally, lots of small-molecule inhibitors have been designed using structure-based approaches [24][25][26][27][28][29].Some of these inhibitors are undergoing clinical evaluation for anticancer treatment [30][31][32][33][34].There is ongoing research focused on developing potent inhibitors that target the MDM2-p53 interaction.Enhancing our understanding of how these inhibitors bind to MDM2 at atomic levels is helpful for the development of potent inhibitors that disrupt the MDM2-p53 interaction, and also provides crucial insight into the structure-affinity relationships of MDM2inhibitor complexes.Various tools are currently available to probe the conformational dynamics of targets, such as conventional molecular dynamics (cMD) [35][36][37][38][39], Gaussian-accelerated molecular dynamics (GaMD) simulations [21,[40][41][42][43], and calculations of binding free energies [44][45][46][47].These methods have been extensively utilized to uncover the molecular mechanisms and free energy bases of target-ligand identification [48][49][50][51][52][53].GaMD simulations, particularly, effectively overcoming energy barriers in protein systems, have shown successes in studying changes in conformational dynamics and the free energy profiles of targets.To extract valuable information from cMD or GaMD trajectories, MD simulations and machine learning (ML) have been used to deepen our understanding of Enhancing our understanding of how these inhibitors bind to MDM2 at atomic levels is helpful for the development of potent inhibitors that disrupt the MDM2-p53 interaction, and also provides crucial insight into the structure-affinity relationships of MDM2-inhibitor complexes.Various tools are currently available to probe the conformational dynamics of targets, such as conventional molecular dynamics (cMD) [35][36][37][38][39], Gaussian-accelerated molecular dynamics (GaMD) simulations [21,[40][41][42][43], and calculations of binding free energies [44][45][46][47].These methods have been extensively utilized to uncover the molecular mechanisms and free energy bases of target-ligand identification [48][49][50][51][52][53].GaMD simulations, particularly, effectively overcoming energy barriers in protein systems, have shown successes in studying changes in conformational dynamics and the free energy profiles of targets.To extract valuable information from cMD or GaMD trajectories, MD simulations and machine learning (ML) have been used to deepen our understanding of the functions of targets and ligand-target binding mechanisms [54][55][56][57][58]. Miao's group developed a GaMD trajectory-based deep learning (DL) approach called the GaMD, DL, and free energy profiling workflow (GLOW) to effectively decode the molecular mechanisms related to GPCR activation and allosteric modulation [59,60].Furthermore, the GaMD trajectory-based DL approach also obtained successes in the ligand-target identification of molecular mech-anisms and the exploration of significant structural domains of targets [61,62].Binding free energy calculations are powerful tools for elucidating the interaction mechanisms between inhibitors and their targets.The MM-GBSA and SIE methods have been identified as effective approaches for calculating the binding free energies of inhibitors to proteins [63][64][65][66].They have been successfully applied to elucidate protein-protein and protein-inhibitor interactions.
In this study, we focus on investigating the binding mechanisms of four inhibitors to MDM2, including two non-peptide inhibitors K23 and 0Y7, and peptide inhibitors PDI6W and PDI.Regarding the two non-peptide inhibitors, K23 and 0Y7, K23 comprises four aromatic groups, efficiently occupying the binding pockets of MDM2 with a median inhibitory concentration (IC50) value of 1.71 µM [67], while 0Y7 is a piperidinone inhibitor that interacts favorably with the N-terminus of human MDM2, with an IC50 value of 50 nM [68].With respect to the two peptide inhibitors, PDI6W and PDI, PDI6W has a residue sequence (LTFEHWWAQLTS) with an IC50 value of 36 nM toward MDM2 [22]; the other peptide inhibitor, PDI, possesses a residue sequence (LTFEHYWAQLTS) showing inhibiting ability on MDM2 with an IC50 of 44 nM [23].The structures of these inhibitors are depicted in Figure 1C-F.Multiple independent Gaussian-accelerated molecular dynamics (MI-GaMD) simulations were conducted to improve conformational sampling, and deep learning (DL) was utilized to identify critical structural domains.Additionally, principal component analysis (PCA) [69][70][71][72] and the construction of free energy landscapes (FELs) were performed to explore conformational spaces and the free energy profile of MDM2.Calculations of binding free energies and energy decomposition analysis were conducted to evaluate the binding difference between the two types of inhibitors of MDM2 and identify hot spots of MDM2 inhibitors.We expect that this work can offer valuable insights into the mechanisms underlying the inhibition of the p53-MDM2 interaction.

Differences in the Contacts of Structural Domains Revealed by Deep Learning
The DL procedure involved several steps: (1) the conformations extracted from the GaMD trajectories were converted into images suitable for DL analysis; (2) these images were randomly split into a training set and a validation set to conduct DL; (3) the outcomes of the DL training were visualized, as shown in Figure 2. Figure 2A shows the classification information and Figure 2B depicts the learning curves of the training and validation datasets.

Free Energy Profiles and Structural Dynamics of MDM2
To gain insights into the dynamic behavior of MDM2, PCA was conducted on the GaMD trajectories of inhibitor-bound MDM2, with the Cα atom coordinates saved in the single GaMD trajectory (SGT) integrated by using three independent GaMD trajectories.The function of eigenvalues over eigenvector indexes was estimated (Figure S1).An eigenvalue is commonly employed to characterize the structural fluctuations of proteins along an eigenvector.The first six eigenvalues accounted for 57.27, 67.96, 56.12, and 51.23% of the total movements of the K23-, 0Y7-, PDI6W-, and PDI-bound MDM2 structures, respectively.Notably, the first eigenvalue of 0Y7-MDM2 is much greater than that of the other inhibitor-bound MDM2 structures, suggesting a greater structural fluctuation amplitude of the 0Y7-bound MDM2 structure along the primary eigenvectors relative to the three other complexes (Figure S1).
To analyze the variations in the free energy profiles of MDM2, SGTs were projected onto the first two eigenvectors, and these projections were used as the reaction coordinates (RCs) to construct FELs, which are illustrated in Figure 3.The GaMD simulations recognized two, four, four, and two energy basins (EBs) in the K23-, 0Y7-, PDI6W-, and PDI-bound MDM2 structures (Figure 3A,C,E,G), respectively, with 0Y7 showing the most pronounced impact on the free energy profiles of MDM2.To investigate the structural variances within different EBs, the representative structures falling into the EBs were superimposed together (Figure 3B,D,F,H).The β-strands and loops located between α2 and α1' of the 0Y7-bound MDM2 structure, including β3, L2, β1', and L3, exhibited the most highly obvious deviation among the various energy states; moreover, this deviation induced changes in the binding poses of 0Y7.It is observed that 0Y7 has four different binding poses (Figure 3D).Although the structures of K23-bound MDM2 in EB1 and EB2 are aligned well, the binding poses of K23 yield sliding and rotation, which induces two different binding orientations (Figure 3B).In spite of the four energy basins (EB1-EB4) of PDI6W-bound MDM2, their representative structures hardly produce obvious deviations; moreover, three key residues (Phe19', Trp23', and Leu26') in PDI6W are also aligned well (Figure 3F).In the PDI-MDM2 complexes, loop L1 deviates between the presentative structures EB1 and EB2; correspondingly, two key residues, Trp23' and Leu26' in PDI, generate obvious sliding (Figure 3H).The aforementioned changes certainly impact the binding of the peptide and non-peptide inhibitors to MDM2.
Molecules 2024, 29, x FOR PEER REVIEW 6 of 2 MDM2, but also effectively prevents p53 from going into the binding pocket of MDM2 The changes in the concerted motions of the structural domains and free energy profile due to the presence of the two types of inhibitors can impact inhibitor-MDM2 binding.In order to probe the effects of the two types of inhibitors on the collective motions of MDM2, the first eigenvector was visualized using the VMD program [73] and the results are depicted in Figure 4. Except for the N-and C-terminals of MDM2, the binding of peptide and non-peptide inhibitors has an obvious effect on the conformations of loops L2 and L1 from MDM2.Compared to the two non-peptide inhibitors (Figure 4A,B), the binding of the two peptide inhibitors, PDI6W and PDI, inhibits the collective motions of loop L2, which implies that these two peptide inhibitors possibly yield interactions with L2, β1', and β3 (Figure 4C,D).In comparison with K23 and 0Y7, the binding of PDI6W changes the fluctuation tendency of loop L1 from MDM2 along the first eigenvector, but the binding of PDI hardly affects the concerted motions of this loop (Figure 4D).In addition, the concerted motions of helix α2' are also affected by inhibitor binding (Figure 4).

Structural Property of MDM2
To investigate the structural stability of MDM2, root-mean-square dev (RMSDs) of backbone atoms from MDM2 were computed throughout the GaMD s tions, with reference to the initially optimized structures (Figure 5A,B).The time of the RMSDs indicates that the structures of MDM2 in the four systems are stable ( Based on the above information, the binding of peptide and non-peptide inhibitors exerts different influences on the free energy profiles of MDM2, the binding poses of inhibitors, and the concerted motions of structural domains.The structural information revealed by the PCA and FELs is in basic agreement with the results from the DL approach.Meanwhile, our current findings are also consistent with the previous results revealed by MD simulations [37,38].According to their structural information, peptide and non-peptide inhibitors bind to the hydrophobic cleft of MDM2 and occupy the binding position of p53 in MDM2, which not only leads to conformational changes in the binding cleft of MDM2, but also effectively prevents p53 from going into the binding pocket of MDM2.The changes in the concerted motions of the structural domains and free energy profiles due to the presence of the two types of inhibitors can impact inhibitor-MDM2 binding.

Structural Property of MDM2
To investigate the structural stability of MDM2, root-mean-square deviations (RMSDs) of backbone atoms from MDM2 were computed throughout the GaMD simulations, with reference to the initially optimized structures (Figure 5A,B).The time course of the RMSDs indicates that the structures of MDM2 in the four systems are stable (Figure 5A).The distributions of RMSDs are estimated in Figure 5B.For the two non-peptide inhibitors, the RMSD of K23-MDM2 is distributed at ~2.75 Å, and that of 0Y7-MDM2 at the peaks of 2.75 and 3.43 Å (Figure 5B).As for the two peptide inhibitors, the RMSDs of PDI6W-MDM2 and PDI-MDM2 are populated at the peaks of 1.61 and 2.98 Å, respectively (Figure 5B).Compared to the two non-peptide inhibitors, the binding of peptide inhibitors decreases the RMSDs of MDM2.This result suggests that the binding of peptide inhibitors is more favorable for the structural stability of MDM2 than non-peptide inhibitors, implying that PDI6W and PDI should generate more interaction contacts with MDM2 than K23 and 0Y7.Meanwhile, our results also show that the GaMD trajectories of the four complexes are reliable for the post-processing analyses.is more favorable for the structural stability of MDM2 than non-peptide inhibitors, implying that PDI6W and PDI should generate more interaction contacts with MDM2 than K23 and 0Y7.Meanwhile, our results also show that the GaMD trajectories of the four complexes are reliable for the post-processing analyses.
To understand the inhibitor-mediated effect on structure flexibility, root-meansquare fluctuations (RMSFs) of MDM2 were estimated using the coordinates of Cα atoms (Figure 5C).In Figure 5C, S1 represents loop L1, and S2 corresponds to the β-strands and loops between α2 and α1', while S3 refers to β2' and the loops between α1' and α2'.To understand the inhibitor-mediated effect on structure flexibility, root-mean-square fluctuations (RMSFs) of MDM2 were estimated using the coordinates of Cα atoms (Figure 5C).In Figure 5C, S1 represents loop L1, and S2 corresponds to the β-strands and loops between α2 and α1', while S3 refers to β2' and the loops between α1' and α2'.It was observed that the S1, S2, and S3 regions of MDM2 exhibit high flexibility, particularly in the case of S2.The analysis highlights that the mobility of the S1, S2, and S3 domains is more pronounced in the 0Y7-MDM2 structure compared to the other structures.With reference to K23-MDM2, the binding of 0Y7 largely enhances the structural flexibility of the structural domains S1, S2, and S3, while the presence of the two peptide inhibitors, PDI6W and PDI, obviously weakens the structural flexibility of these two structural domains (Figure 5C,D), implying the different binding abilities of the four inhibitors to MDM2.
To assess alterations in the secondary structures of MDM2 and the peptide inhibitors, a combination of the program CPPTRAJ and DSSP second-structure analysis [74] was used to investigate the changes in the secondary structures in three separate GaMD simulations.The time evolutions of the secondary structures for the K23-, 0Y7-, PDI6W-and PDI-bound MDM2 structures are displayed in Figure 6A-D Based on the above analyses, the binding of peptide and non-peptide inhibitor yields different impacts on conformations of MDM2: (1) the binding of peptide inhibitor leads to more stable structures of MDM2 than non-peptide inhibitors; (2) the presence o peptide inhibitors induces more rigid structures in the S1, S2, and S3 domains.These re sults are in basic agreement with previous studies [35][36][37][38].

Comparative Calculations of Binding Free Energies
The binding free energies of K23, 0Y7, PDI6W, and PDI to MDM2 were assessed Based on the above analyses, the binding of peptide and non-peptide inhibitors yields different impacts on conformations of MDM2: (1) the binding of peptide inhibitors leads to more stable structures of MDM2 than non-peptide inhibitors; (2) the presence of peptide inhibitors induces more rigid structures in the S1, S2, and S3 domains.These results are in basic agreement with previous studies [35][36][37][38].

Comparative Calculations of Binding Free Energies
The binding free energies of K23, 0Y7, PDI6W, and PDI to MDM2 were assessed through MM-GBSA and SIE calculations so as to understand the binding preferences of the two types of inhibitors, and the results are shown in Tables 1 and 2. In our calculations, 400 snapshots were extracted from the equilibrated cMD trajectories.It was found that the ranking of the calculated binding free energies was consistent with that determined by the known experimental data, implying the reliability of our free energy analyses [22,23,35,36,67,68].As illustrated in Table 1, the binding free energies of K23, 0Y7, PDI, and PDI6W to MDM2 were −10.31, −12.61, −21.96, and −24.67 kcal•mol −1 , respectively.Analysis of the energy components revealed that van der Waals energies (∆E vdW ) are the primary favorable contributors to inhibitor binding.Non-polar solvation energies (∆G sur f ) also provide favorable contributions to the binding process.However, the contributions of entropy changes (−T∆S) to free energies greatly impair the binding of the four inhibitors to MDM2.Electrostatic terms (∆E ele ) also favor inhibitor binding, but polar solvation energies (∆G gb ) provide opposite contributions for the binding of inhibitors.On the whole, the sum of ∆E ele and ∆G gb in the current complexes is unfavorable for the binding of inhibitors.Therefore, two favorable forces, ∆E vdW and ∆G sur f , mostly drive the binding of the four inhibitors to MDM2, which agrees well with previous experimental analyses [35][36][37][38].This result suggests that the optimization of van der Waals interactions and non-polar solvation energies may lead to the potent inhibition of MDM2.
According to Table 2, the calculated binding free energies for K23, 0Y7, PDI, and PDI6W to MDM2 using the SIE method are −7.19,−7.25, −10.53, and −10.79 kcal•mol −1 , respectively.More importantly, the ranking of binding free energies predicted by the SIE method also agrees with that of the experimental values, which further verifies the reliability of our current results.The reaction energy (∆G R ) associated with the desolvation of polar groups always impedes inhibitor bindings (Table 2).As seen in Table 2, the unfavorable reaction energy of the polar groups is partially compensated by the favorable intermolecular Coulomb interaction (∆E c ).Additionally, intermolecular VDWIs (∆E vdW ) also provide partial compensation to this unfavorable effect.
Through the above results, the binding free energies predicted by the MM-GBSA and SIE methods show that the four current inhibitors yield tight associations with MDM2, implying that they have strong competitive ability in binding to MDM2 relative to p53.Furthermore, the VDWIs of the inhibitors with MDM2 play key roles in the binding of the two types of inhibitors, which implies that the hydrophobic groups of the peptide and nonpeptide inhibitors, such as alkyls and aromatic rings, are significant molecular structures to be considered in future drug design in relation to p53-MDM2 interactions.Meanwhile, VDWIs should be paid more attention in the development of clinically available inhibitors targeting p53-MDM2 interactions.

Interaction Network of Inhibitors with MDM2
To understand the contributions of individual residues to inhibitor-MDM2 binding, a residue-based free energy decomposition method was utilized to estimate inhibitor-residue interactions and residue-residue interactions (Figures 7, 8, S2 and S3).The probability distributions of the distances associated with crucial inhibitor-residue interactions were analyzed (Figures 9 and S4).Furthermore, the geometric aspects of hydrophobic interactions and hydrogen bonding interactions (HBIs) were visually illustrated for better comprehension (Figures 7C,D, 8A,B and S3A,B).
In the K23-MDM2 complex, K23 establishes interactions of strength exceeding 1 kcal/mol with six residues, specifically Leu54, Leu57, Gly58, Ile61, Val93, and Ile99 (Figures 7A,B and S2A).As illustrated in Figure 7C, the alkyl groups of residues Leu54 and Ile99 are close to the hydrophobic ring R1 of K23, leading to the formation of CH-π interactions between them.The respective interaction energies for these interactions are −3.03 and −1.46 kcal/mol for Leu54 and Ile99 (Figure 7A).The distances for the mass centers of the sidechains of Leu54 and Ile99 away from the mass center of the ring R1 are, respectively, distributed at 4.63 and 5.68 Å (Figure 7E), which verifies the hydrophobic interactions of these two residues with K23.The alkyl groups of Leu57 and the CH group of Gly58 form multiple CH-π contacts with the R2 ring of K23 (Figure 7C), leading to CH-π interaction energies of −1.52 and −1.37 kcal/mol for residues Leu57 and Gly58, respectively (Figure 7A,B).The distances for the carbon atoms of the alkyl groups of Leu57 and Gly58 from the mass center of the ring R2 are distributed at 5.68 and 4.28 Å (Figure 7E), respectively, which verifies the existence of CH-π interactions.The alkyl group of Val93 engages in the CH-π contacts with the hydrophobic ring R3 of K23, while the alkyl moiety of Ile61 forms CH-π contacts with the hydrophobic group R4 of K23 (Figure 7C).As a result, Val93 and Ile61 separately provide the interaction energies of −2.07 and −1.78 kcal/mol for the K23-MDM2 binding (Figure 7A,B).The distance between the center of mass of the alkyl group from Val93 and that of ring R3 is situated at 4.28 Å (Figure 7E), indicating the presence of a CH-π interaction between K23 and Val93.Similarly, the distance between the center of mass of the alkyl group from Ile61 and that of ring R4 is 4.98 Å (Figure 7E), demonstrating the existence of the CH-π interaction between K23 and Ile61.Moreover, K23 establishes a hydrogen bonding interaction with Leu54 with an occupancy exceeding 98.58%, suggesting the stability of this hydrogen bond throughout the GaMD simulations (Table 3).These analyses align with previous statistical examinations and other research findings [35][36][37][38].
a residue-based free energy decomposition method was utilized to estimate inhibitor-re idue interactions and residue-residue interactions (Figures 7, 8, S2 and S3).The probab ity distributions of the distances associated with crucial inhibitor-residue interactio were analyzed (Figures 9 and S4).Furthermore, the geometric aspects of hydrophobic i teractions and hydrogen bonding interactions (HBIs) were visually illustrated for bett comprehension (Figures 7C,D, 8A,B and S3A,B).a The hydrogen bonds are determined by a donor-acceptor atom distance of <3.5 Å and acceptor-H-donor angle of >120 • .b Occupancy is used to evaluate the stability and strength of the hydrogen bond.
In the 0Y7-MDM2 complex, 0Y7 produces interactions stronger than 1.0 kcal/mol with six residues, including Leu54, Leu57, Gly58, Ile61, Val93, and Ile99 (Figures 7A and S2B).The alkyl groups of Leu57, Gly58, ILE61, and ILE99 establish several CH-π contacts with the hydrophobic ring R2 of 0Y7, resulting in interaction energies of −1.12, −1.14, −1.46, and −1.65 kcal/mol for Leu57, Gly58, ILE61, and ILE99, respectively (Figure 7A,B,D).The distances of the mass centers for the alkyl groups of Leu57, Gly58, ILE61, and ILE99 away from the ring R2 are situated at 5.33, 4.63, 4.63, and 4.63 Å, individually, as depicted in Figure 7F.The alkyl group of Leu54 forms CH-π contacts with the ring R1 of 0Y7, and its corresponding interaction energy is −2.04 kcal/mol (Figure 7A,B).The distance between the center of mass of the alkyl group from Leu54 and that of ring R1 is 4.98 Å (Figure 7F).The alkyl of Val93 structurally forms a CH-CH interaction with that of 0Y7, which provides an energy contribution of −1.65 kcal/mol to the 0Y7-MDM2 binding (Figure 7A,D).Moreover, the distance of 3.93 Å between the carbon atom of the alkyl group from Val93 and that of 0Y7 verifies the existence of the CH-CH interaction (Figure 7F).
PDI6W-MDM2 interactions are illustrated with the relative positions between key residues from PDl6W and MDM2 (Figure 8).Seven residues of PDl6W can strongly interact with MDM2.The Thr18' of the PDl6W inhibitor engages in a significant CH-CH interaction of −3.02 kcal/mol with Gln72 of MDM2 (Figure 8A,C), while the distance between the center of mass of the alkyl group from Gln72 and that of Thr18' is 4.34 Å (Figure 9A).The Phe19' of the PDl6W inhibitor exhibits strong interactions with four MDM2 residues, including Ile61, Met62, Tyr67, and Gln72.Specifically, the interaction energy between Phe19' and Tyr67 is −2.40 kcal/mol, aligning with the π-π interaction between the benzene rings of these two residues in the spatial structure (Figure 8A,D).The interaction energies between Phe19' and residues Ile61, Met62, and Gln72 are −1.61,−1.23, and −3.32 kcal/mol, respectively (Figure 8D).These energy contributions primarily arise from the CH-π interactions between the CH groups of residues Ile61, Met62, and Gln72 and the benzene ring of Phe19' (Figure 8A).As illustrated in Figure 9B, the distances between the mass centers of the alkyl groups of Ile61, Met62, and Tyr67 and that of Phe19' are 3.88, 5.88, and 5.13Å, respectively, which demonstrates the existence of the aforementioned CH-π interactions.A hydrogen bond with an occupancy of 76.31% appears between Phe19' and Gln72 (Table 3 and Figure 8B), implying that this hydrogen bond is stable.Structurally, the CH groups of Gln72, Val93, and Lys94 engage in the CHπ interactions with Trp22' of the pDl6W inhibitor (Figure 8A), resulting in interaction energies of −1.24, −1.41, and −2.18 kcal/mol, respectively (Figure 8E).His73 establishes a π-π interaction of −2.47 kcal/mol with Trp22' (Figure 8A,E).The distances for the mass centers of the alkyl groups of Gln72, His73, Val93, and Lys94 from those of Trp22' are separately distributed at 4.88, 5.88, 6.13, and 4.88 Å (Figure 9C).The most robust interaction is observed between Trp23' and Leu54, with an interaction energy of −3.23 kcal/mol (Figure 8F).This energy primarily originates from two sources: (1) the alkyl group of Leu54 forms a CH-π interaction with Trp23' (Figure 8A), and (2) a hydrogen bond with a 97.62% occupancy is formed between Leu54 and Trp23' (Table 3 and Figure 8B).The interaction energies between Trp23' and Leu57, Gly58, Ile61, and Val93 are −1.47,−1.28, −1.14, and −1.61 kcal/mol, respectively (Figure 8F), reflecting the structural agreement of the CH-π interactions between the CH groups of these four residues and the hydrophobic ring of Trp23' (Figure 8A).The distances between the mass centers of the alkyl groups of Leu54, Leu57, Gly58, Ile61, and Val93 and those of Trp23' are populated at 5.13, 5.13, 4.63, 4.38, and 5.13 Å (Figure 9D), which supports the abovementioned CH-π interactions involving Trp23'.The Leu26' residue of the PDl6W inhibitor exhibits significant interactions with five MDM2 residues, including Leu54, Val93, His96, Ile99, and Tyr100.The interaction energies between Leu26' and Leu54, Val93, and Ile99 are −1.33,−1.30, and −1.20 kcal/mol, respectively (Figure 8G), and are predominantly driven by the CH-CH interaction between the alkyl groups of Leu54, Val93, and Ile99 and the alkyl group of Leu26' (Figure 8A).Additionally, the interaction energies between Leu26' and His96 and Tyr100 are −1.94 and −1.89 kcal/mol (Figure 8G), which is consistent with the CH-π interaction of His96 and Tyr100 with the alkyl groups of Leu26' (Figure 8A).In Figure 9E, the distances between the mass centers of the hydrophobic groups of Leu54, Val93, His96, Ile99, and Tyr100 and those of Trp23' are located at 5.13, 5.38, 5.13, 3.88, and 4.63 Å, respectively, which further reveals key residues interacting with Trp23'.The interaction energies of Thr27' in PDI6W with residues Lys51 and Leu54 are −2.21 and −1.28 kcal/mol, coming structurally from the CH-CH interactions (Figure 8A,H).The distances between the mass centers of the alkyl groups of Lys51 and Leu54 from those of Trp26' are 4.88 and 3.88Å, as illustrated in Figure 9F, implying the existence of these key interactions.PDI6W-MDM2 interactions are illustrated with the relative positions between k residues from PDl6W and MDM2 (Figure 8).Seven residues of PDl6W can strongly inte act with MDM2.The Thr18' of the PDl6W inhibitor engages in a significant CH-CH inte action of −3.02 kcal/mol with Gln72 of MDM2 (Figure 8A,C), while the distance betwee the center of mass of the alkyl group from Gln72 and that of Thr18' is 4.34 Å (Figure 9A The Phe19' of the PDl6W inhibitor exhibits strong interactions with four MDM2 residue including Ile61, Met62, Tyr67, and Gln72.Specifically, the interaction energy betwee Phe19' and Tyr67 is −2.40 kcal/mol, aligning with the π-π interaction between the benzen rings of these two residues in the spatial structure (Figure 8A,D).The interaction energi between Phe19' and residues Ile61, Met62, and Gln72 are −1.61,−1.23, and −3.32 kcal/mo respectively (Figure 8D).These energy contributions primarily arise from the CH-π inte actions between the CH groups of residues Ile61, Met62, and Gln72 and the benzene rin of Phe19' (Figure 8A).As illustrated in Figure 9B, the distances between the mass cente of the alkyl groups of Ile61, Met62, and Tyr67 and that of Phe19' are 3.88, 5.88, and 5.13 respectively, which demonstrates the existence of the aforementioned CH-π interaction A hydrogen bond with an occupancy of 76.31% appears between Phe19' and Gln72 (Tab 3 and Figure 8B), implying that this hydrogen bond is stable.Structurally, the CH grou of Gln72, Val93, and Lys94 engage in the CH-π interactions with Trp22' of the pDl6 inhibitor (Figure 8A), resulting in interaction energies of −1.24, −1.41, and −2.18 kcal/mo respectively (Figure 8E).His73 establishes a π-π interaction of −2.47 kcal/mol with Trp2 As seen in Figures S3 and S4, the PDI inhibitor binds to MDM2 in a mode similar to pDl6W.A hydrogen bond with an occupancy of 81.36% appears between Phe19' and Gln72 (Table 3 and Figure S3B).Additionally, Leu54 forms a hydrogen bond with Trp23', and its occupancy is 95.73% (Table 3 and Figure S3B).Compared to the residue-residue interaction in PDI6W, the interaction energy between Tyr22' and Lys94 is lower than 1.0 kcal/mol, which indicates that there is a weak CH-π contact between the alkyl moiety of Lys94 and Tyr22' (Figure S3E).The probability distributions of the distances relating to key inhibitorresidue interactions for PDI-MDM2 are illustrated in Figure S4.The work by Liu et al. showed that D-peptide inhibitors can form favorable interactions with residues Leu54, Leu57, Ile61, Tyr67, Val93, His96, and Tyr100, which agrees well with our current calculated results [75,76].The study of Strizhak et al. revealed that a stapled peptide produces interactions with Leu54, Ile61, Tyr67, Gln72, Val93, His96, and Ile99, supporting our current findings [77].In previous experimental and theoretical works [14], p53 has yielded strong interactions with Leu54, Leu57, Gly58, Ile61, Val93, Ile99, etc., which not only supports our current results well, but also implies that four inhibitors occupy the binding sites of p53 in MDM2, impairing the binding of p53.Moreover, the sidechains of key residues revealed by our calculations play vital roles in the binding of inhibitors to MDM2.More importantly, the CH-π, CH-CH, and π-π interactions between individual residues of MDM2 and the inhibitors drive the binding of K23, 0Y7, PDl6W, and PDI to MDM2, which should be paid special attention in future drug design in relation to p53-MDM2 interactions.Therefore, it is of high significance to rationally optimize the interactions of inhibitors with the sidechains of key residues in MDM2 for the design of efficient inhibitors.

System Preparation
The initial atomic coordinates of the K23-, 0Y7-, PDI6W-, and PDI-MDM2 complexes were taken from the Protein Data Bank (PDB), corresponding to PDB entries 3LBK, 4HBM, 3JZR, and 3G03, respectively.As there are differences in the residue sequences of MDM2, residues of Thr26-Arg105 in MDM2 were utilized to construct our simulation systems.The residues from the two peptide inhibitors, PDI6W and PDI, were canonical amino acids.The protonation states of the MDM2 residues were validated using the H++ 3.0 program [78].Rational protonation states were assigned to each MDM2 residue, and any missing hydrogen atoms in the crystal structures were added using the Leap module in Amber20 [79,80].The structures of two inhibitors, K23 and 0Y7, were optimized at a semi-empirical AM1 level, and subsequently, BCC charges [81] were assigned to each atom of the two inhibitors using the Antechamber module in Amber20 [82].The ff 19SB force field [83] was employed to parameterize MDM2, while the general Amber force field (GAFF2) [84,85] was used to derive force-field parameters for K23 and 0Y7.The systems, comprising K23-MDM2, 0Y7-MDM2, PDI6W-MDM2, and PDI-MDM2, were solvated in an octahedral periodic water box with a 10.0 Å buffer to mimic a solvent environment.The force-field parameters for the water molecules were based on the TIP3P model [86,87].To ensure neutral simulation systems, an appropriate number of sodium ions (Na + ) and chloride ions (Cl − ) at a concentration of 0.15 M NaCl were added to the water box.The parameters for the Na + and Cl − ions were adopted from the work of Joung et al. [88,89].

Multiple Independent Gaussian-Accelerated Molecular Dynamics
To address potential bad contacts between atoms arising from the initialization of the four current MDM2-related systems, each system underwent a two-step minimization process.This included a steepest descent minimization followed by a 10,000cycle conjugate gradient minimization.The optimized systems were gradually heated from 0 to 300 K over 1 ns in the canonical ensemble (NVT), utilizing a weak harmonic restraint of 2 kcal•mol −1 •Å 2 on heavy atoms.Subsequently, the four systems were further equilibrated at 300 K under an isothermal−isobaric ensemble (NPT).A 2 ns NPT simulation was then conducted to maintain the system density at 1.01 g/cm 3 .Finally, three 2.4 ns independent cMD simulations were separately performed on the four systems during the NVT with periodic boundary conditions using the particle mesh Ewald method (PME).In each independent cMD simulation, the initial atomic velocities were randomly assigned with the Maxwell distribution.
The well-equilibrated systems served as the starting points for three independent Gaussian-accelerated molecular dynamics (GaMD) simulations.GaMD simulations employ a harmonic boost potential to reduce free energy barriers in biomolecules and enhance the conformational sampling of systems.In GaMD simulations, if the potential energy V( ⇀ r ) of the system is lower than a threshold energy E, V( ⇀ r ) is revised to V * ( ⇀ r ), according to Equations ( 1) and ( 2) below: ∆V( In the above equations, the parameter k represents the harmonic force constant.The parameters E and k can be adjusted according to the enhanced sampling principles defined in Equations ( 3) and (4), as shown below: If E is designated as the lower bound (E = V max ), then k 0 can be determined using Equation ( 5): On the contrary, if E is set as the upper bound (E = V min + 1 k ), then k 0 is derived from Equation (6): where the three energy parameters V max , V min , and V avg indicate the maximum, minimum, and averaged potential energies of the simulated systems extracted from the previous cMD simulations, respectively.The parameter σ V represents the standard deviation of the system's potential energies, and σ 0 is a user-defined upper limit for accurate reweighting.
In our current study, 1.2 µs GaMD simulations, composed of three independent simulations of 400 ns, were separately performed on four current MDM2-related systems.To facilitate deep learning (DL) and post-processing analyses, three independent GaMD trajectories were combined into a single GaMD trajectory (SGT), and the CPPTRAJ module integrated with Amber was used to extract data insights into the function of MDM2-related systems.A program called PyReweighting, developed by Miao et al. [90], was utilized to accurately reweight and identify the original free energy profiles of our current systems.All cMD and GaMD simulations employed the SHAKE algorithm to constrain the chemical bonds between the hydrogen atoms and heavy atoms [91].The Langevin thermostat, bringing a collision frequency of 2.0 ps −1 , was utilized to tune the temperatures of the four MDM2-related systems [92].Non-bonded interactions were estimated using the particle mesh Ewald (PME) method [93] with a 12 Å cutoff.The simulations were executed using the pmemd.cudaprogram implemented in Amber20 [94,95].

Deep Learning
To investigate the impacts of inhibitor binding on the internal structures of MDM2, DL was utilized to identify differences in residue contacts.The residue contact maps for each snapshot of MDM2 were computed using the Python packages MDTraj and Contact Map Explorer [96].Contact was defined as ≤4.5 Å between any Cα atoms of two proteins.The resulting 80 × 80 residue contacts were converted into grayscale images of 80 × 80 pixels for analysis by a two-dimensional (2D) convolutional neural network (CNN).A total of 160,000 images were generated for each MDM2-related system, with 80% randomly selected for training and the remaining 20% used for validation.The 2D-CNN model was constructed using the PyTorch package, consisting of three convolutional layers with a 1 × 1 kernel size and 16, 32, and 32 filters, followed by three fully connected layers.The first two fully connected layers comprised 512 and 128 filters with a dropout rate of 0.5 each, while the final fully connected layer served as the classification layer for inhibitor-bound MDM2.
Throughout the 2D-CNN architecture, the "ReLu" activation function was employed in all layers, with the "softmax" activation function used at the classification layer.A maximum pooling layer with a 2 × 2 kernel size was added after each convolutional layer.Backpropagation via vanilla gradient-based pixel attribution [97] was utilized to estimate an attention map of the residue contact gradients to aid in discriminating the functional differences in MDM2 induced by inhibitor binding.The residue contact map was represented using the most populated structural cluster of each MDM2-related system.Our DL program was rewritten using PyTorch based on the work of Miao's group [59].

Principal Component Analysis and Dynamic Cross-Correlation Maps
PCA is a crucial technique for deciphering conformational changes in proteins.In this research, PCA was conducted by diagonalizing the covariance matrix C constructed from the Cα atom coordinates of MDM2, as outlined in Equation ( 7): C =< (q i − < q i >) q j − < q j > T > (7) where q i and q j represent the Cartesian coordinates of the ith and jth Cα atoms from MDM2, while <q i > and <q j > denote their average positions across conformational ensembles obtained from MIGaMD simulations.The eigenvector and the eigenvalue generated by the diagonalization of the covariance matrix characterize the concerted movement of the structural domains and the fluctuation amplitude along a given eigenvector, respectively.For this study, PCA was conducted using the CPPTRAJ program [98] in the Amber suite.

Construction of Free Energy Landscapes
To investigate the influences of peptide and non-peptide inhibitors on the free energy profiles of MDM2, projections of GaMD trajectories onto the first two eigenvectors served as RCs for constructing free energy landscapes (FELs).During the reweighting process in GaMD simulations, reweighted free energy F(A) = −k B Tln(ρ A ) is calculated as where F * (A) = −k B Tlnp * (A) the modified free energy obtained from the GaMD simulations, F C denotes a constant, and β = k B T. The probability distribution p * (A) of the selected RCs from the GaMD simulations can be reweighted to match the canonical ensemble distribution ρ A .All free energy reweighting calculations were performed using the PyReweighting program developed by Miao et al.A detailed description for the reweighting is given in the work of Miao et al. [90].

Binding Free Energies
To evaluate the binding of the two types of inhibitors to MDM2, MM-GBSA and SIE methods were adopted to calculate binding free energies.In the MM-GBSA method, enthalpy changes (∆H) and entropy changes (-T∆S) play essential roles in ligand associations.The binding free energies of the non-peptide inhibitors K23 and 0Y7 and the peptide inhibitors PDI6W and PDI to MDM2 were calculated using the MM-GBSA method based on the following Equation ( 9): where ∆H is calculated using Equation ( 10): where electrostatic interactions (EIs, ∆E ele ) and van der Waals interactions (VDWIs, ∆E vdW ) can be estimated using molecular mechanics and the ff 19SB force field.Polar solvation free energies (PSFEs, ∆G gb ) were estimated using the generalized Born (GB) model proposed by Onufriev et al. [99].Non-polar solvation free energies (NPSFEs, ∆G sur f ) were calculated based on the following empirical Equation (11): where the term ∆SASA represents the variation in the solvent-accessible surface area (SASA) mediated by the binding of inhibitors.Entropy changes (-T∆S) were computed using the MMPBSA.pyprogram within the Amber20 software [100].The two parameters γ and β were set as 0.0072 kcal•mol•Å −2 and 0.0 kcal•mol −1 , respectively [101].
In the SIE method, the SIE function [66] to calculate inhibitor-protein binding free energies is expressed as follows in Equation ( 12): where E c and E vdW represent the intermolecular Coulomb and van der Waals interaction energies in the bound state, respectively.∆G R signifies the change in the reaction field energy caused by the binding of an inhibitor, which was determined by solving the Poisson equation using the boundary element method (BRI BEM) [102,103] and a solvent probe with a variable radius of 1.4 Å [104].γ•∆MSA corresponds to the change in the molecular surface area upon binding.The parameters ρ, D in , γ, and C are the Amber van der Waals radii linear scaling coefficient, the solute interior dielectric constant, the molecular surface area coefficient, and a constant, respectively.The parameter α relates to the overall proportionality coefficient associated with the loss of conformational entropy upon binding [105].

Figure 3 .
Figure 3. Free energy profiles and representative structures of inhibitor-bound MDM2: (A,C,E,G correspond to the FELs of K23-MDM2, 0Y7-MDM2, PDI6W-MDM2, and PDI-MDM2, respectively (B,D,F,H) indicate the superimposition of representative structures of MDM2 and inhibitor trapped in different EBs.The free energy is scaled in kcal/mol.

Figure 3 .
Figure 3. Free energy profiles and representative structures of inhibitor-bound MDM2: (A,C,E,G) correspond to the FELs of K23-MDM2, 0Y7-MDM2, PDI6W-MDM2, and PDI-MDM2, respectively; (B,D,F,H) indicate the superimposition of representative structures of MDM2 and inhibitors trapped in different EBs.The free energy is scaled in kcal/mol.
It was observed that the S1, S2, and S3 regions of MDM2 exhibit high flexibility, particularly in the case of S2.The analysis highlights that the mobility of the S1, S2, and S3 domains is more pronounced in the 0Y7-MDM2 structure compared to the other structures.With reference to K23-MDM2, the binding of 0Y7 largely enhances the structural flexibility of the structural domains S1, S2, and S3, while the presence of the two peptide inhibitors, PDI6W and PDI, obviously weakens the structural flexibility of these two structural domains (Figure 5C,D), implying the different binding abilities of the four inhibitors to MDM2.

Figure 5 .
Figure 5.The RMSD and RMSF values of MDM2 in GaMD simulations: (A) the time course of RMSDs for MDM2; (B) the probability distribution of RMSDs for MDM2; (C) the RMSFs of the Cα atoms from MDM2; and (D) the flexibility domains revealed by GaMD simulations.To assess alterations in the secondary structures of MDM2 and the peptide inhibitors, a combination of the program CPPTRAJ and DSSP second-structure analysis[74] was used to investigate the changes in the secondary structures in three separate GaMD simulations.The time evolutions of the secondary structures for the K23-, 0Y7-, PDI6W-and PDI-bound MDM2 structures are displayed in Figure6A-D, individually.The time evo-

Figure 5 .
Figure 5.The RMSD and RMSF values of MDM2 in GaMD simulations: (A) the time course of RMSDs for MDM2; (B) the probability distribution of RMSDs for MDM2; (C) the RMSFs of the Cα atoms from MDM2; and (D) the flexibility domains revealed by GaMD simulations.

2 Figure 6 .
Figure 6.Stability of secondary structures of MDM2 and peptide inhibitors: (A) time evolution o secondary structure for MDM2 in the K23-MDM2 complex; (B) time evolution of secondary struc ture for MDM2 in the 0Y7-MDM2 complex; (C) time evolution of secondary structure for MDM2 i the PDI6W-MDM2 complex; (D) time evolution of secondary structure for MDM2 in the PDI-MDM complex; (E) time evolution of secondary structure for the peptide inhibitor PDI6W; (F) time evolu tion of secondary structure for the peptide inhibitor PDI.

Figure 6 .
Figure 6.Stability of secondary structures of MDM2 and peptide inhibitors: (A) time evolution of secondary structure for MDM2 in the K23-MDM2 complex; (B) time evolution of secondary structure for MDM2 in the 0Y7-MDM2 complex; (C) time evolution of secondary structure for MDM2 in the PDI6W-MDM2 complex; (D) time evolution of secondary structure for MDM2 in the PDI-MDM2 complex; (E) time evolution of secondary structure for the peptide inhibitor PDI6W; (F) time evolution of secondary structure for the peptide inhibitor PDI.

Figure 7 .
Figure 7. Interaction network of inhibitors K23 and 0Y7 with MDM2: (A) the key residues playi important roles in the binding of inhibitors to MDM2; (B) the radar representation of inhibito residue interactions in MDM2; (C) hydrophobic interactions of inhibitors with residues in the K2 MDM2 complex; (D) hydrophobic interactions of inhibitors with residues in the 0Y7-MDM2 com plex; (E) the probability distribution of the distances between K23 and the key residues in MDM and (F) the probability distribution of the distances between 0Y7 and the key residues in MDM2.

Figure 7 .
Figure 7. Interaction network of inhibitors K23 and 0Y7 with MDM2: (A) the key residues playing important roles in the binding of inhibitors to MDM2; (B) the radar representation of inhibitor-residue interactions in MDM2; (C) hydrophobic interactions of inhibitors with residues in the K23-MDM2 complex; (D) hydrophobic interactions of inhibitors with residues in the 0Y7-MDM2 complex; (E) the probability distribution of the distances between K23 and the key residues in MDM2; and (F) the probability distribution of the distances between 0Y7 and the key residues in MDM2.

Molecules 2024 ,Figure 8 .
Figure 8. Interaction network of PDI6W inhibitor with MDM2: (A) relative geometric positions the key residues in the pDl6W-MDM2 complex; (B) hydrogen bonding interactions in the pDl6W MDM2 complex; and (C-H) residue-residue interaction spectrum between the pDl6W inhibitor an MDM2.

Figure 8 .
Figure 8. Interaction network of PDI6W inhibitor with MDM2: (A) relative geometric positions of the key residues in the pDl6W-MDM2 complex; (B) hydrogen bonding interactions in the pDl6W-MDM2 complex; and (C-H) residue-residue interaction spectrum between the pDl6W inhibitor and MDM2.
, 5.13, 3.88, and 4.63 Å, respectively, which further reveals key residues interacting with Trp23'.The interaction energies of Thr27' in PDI6W with residues Lys51 and Leu54 are −2.21 and −1.28 kcal/mol, coming structurally from the CH-CH interactions (Figure8A,H).The distances between the mass centers of the alkyl groups of Lys51 and Leu54 from those of Trp26' are 4.88 and 3.88Å, as illustrated in Figure9F, implying the existence of these key interactions.

Figure 9 .
Figure 9. Probability distributions of the distances relating to key inhibitor-residue interactions for the PDI6W-MDM2 complex: (A) the distance between Thr18' and the key residue in MDM2; (B) the distances between Phe19' and the key residues in MDM2; (C) the distances between Trp22' and the key residues in MDM2; (D) the distances between Trp23' and the key residues in MDM2; (E) the

Figure 9 .
Figure 9. Probability distributions of the distances relating to key inhibitor-residue interactions for the PDI6W-MDM2 complex: (A) the distance between Thr18' and the key residue in MDM2; (B) the distances between Phe19' and the key residues in MDM2; (C) the distances between Trp22' and the key residues in MDM2; (D) the distances between Trp23' and the key residues in MDM2; (E) the distances between Leu26' and the key residues in MDM2; and (F) the distances between Thr27' and the key residues in MDM2.
Note: Standard errors are given in parentheses.a
Note: Standard errors are given in parentheses.a

Table 3 .
The hydrogen bonds formed between key residues and inhibitors of MDM2.