Discovery of Novel Inhibitors Targeting HDM2: Virtual Screening, Molecular Dynamics Simulation, Binding affinity predication and Cell Assay

doi:10.21203/rs.3.rs-1410429/v1

PDF

Research Article

Discovery of Novel Inhibitors Targeting HDM2: Virtual Screening, Molecular Dynamics Simulation, Binding affinity predication and Cell Assay

https://doi.org/10.21203/rs.3.rs-1410429/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

HDM2(human double-minute 2 protein) is an important drug target for tumor treatment. The main physiological function of HDM2 is to negatively regulate the tumor suppressor protein p53 and inhibit the normal biological function of p53, which leads to the occurrence and development of tumors. Therefore, the discovery of small molecules that can inhibit the activity of HDM2 is a promising cancer therapy. In this study, we reported a virtual screening as an effective strategy to discover potential HDM2 inhibitors from the Specs database. Then 100 ns all-atom molecular dynamics (MD) simulation and binding free energies calculated were carried out on the selected compounds, and cell assay showed that Hit 3 was the most active against SMMC-7721 cell line with an IC₅₀ of 14.14μM. In conclusion, we reported seven prospective candidates with novel skeletons among them Hit 3 is a potential HDM2 inhibitor and provided an insight into the mechanism of interaction of the compounds to HDM2 target.

HDM2 inhibitors

Virtual screening

Molecular dynamics simulation

Binding free energy

Cell assay

Cancer has always afflicted humans with increasing incidence and mortality, with an estimated 9.6 million deaths in 2018 and the number expected to rise to more than 14 million annually over the next decade[1, 2]. In 2020, there are 4.569 million new cancer cases and 3.003 million cancer deaths in China for the Cancer Today website (https://gco.iarc.fr/today). Cancer has become the leading cause of death in China and is a major public health problem.

The tumor suppressor protein p53 is known as the “guardian of the genome” and plays a central role in preventing tumor development through biological functions such as cell cycle arrest, DNA repair, cell senescence, and cell apoptosis[3]. In about half of cancer cases, p53 is inactivated by mutations or deletions. In many other cancer cases, it is inactivated by overexpression of HDM2[4].

HDM2 (also known as MDM2) is the main negative regulator for p53. In normal cells, HDM2 and p53 form a negative feedback loop that maintains both p53 and HDM2 at very low levels[5]. In this process, HDM2 is transcriptionally activated by p53, and HDM2 regulates p53 activity through three main mechanisms. (i) HDM2 inhibits p53 transcriptional activity by binding to the transactivation domain of p53 and thereby inhibits p53-mediated transactivation. (ii) HDM2 exports p53 out of the cell nucleus, leaving p53 inaccessible to targeted DNA and reducing its transcriptional ability. (iii) HDM2 functions as an E3 ubiquitin ligase that promotes the degradation of p53 proteasome[3, 6, 7]. Consequently, HDM2 is considered as an effective inhibitor of p53. HDM2 is abnormally up-regulated in human tumors through gene amplification, increased transcription level and enhanced translation[8], and the overexpression of HDM2 can impair the stability and activities of p53 and have carcinogenic effects. An analysis of nearly 4000 samples from tumors or xenografts from 28 tumor types shows that the frequency of HDM2 amplification in these human tumors was 7%[8]. Therefore, the HDM2 as a drug target that the inhibition of HDM2 can restore the activity of p53.

The HDM2-p53 interaction has been confirmed to be the interaction between the first 120 N-terminal residues of HDM2 and 30 N-terminal residues of p53[9]. Unlike many other protein-protein interfaces, the p53-HDM2 interface forms a narrow and deep cleft[10]. In 1996, the first cocrystal structure of HDM2(residues 17 − 125) with p53(residues 15 − 29) was reported (PDB code 1YCR) [4]. The interaction between HDM2 and p53 is mainly that the three hydrophobic residues of p53 (Phe19, Trp23 and Leu26) interact with the surface of HDM2 to form three hydrophobic cavities. The key residues of these three mimic ligands of p53 occupy the hydrophobic cleft of HDM2 in a “thumb − index − middle finger” conformation.[11]. And some of the other ligands were designed for ubiquitin ligase. The binding modes of all ligands are highly conservative and occupying three hydrophobic groups with a strict arrangement.

In recent years, some highly effective and selective HDM2 inhibitors have been successfully discovered, such as RG7112[12], RG7338[13], SAR405838[14], HDM201[15], AMG-232[16] and CGM097 from Novartis [17]. Other compounds such as MI-219[18] also have strong inhibitory effects on HDM2, as well. However, there are still no relevant medicines used in clinic.

Therefore, designing small molecule inhibitors to block HDM2-p53 interaction has become a possible cancer treatment strategy. In order to develop new HDM2 inhibitors, we have adopted the computer-aided drug design method in this work. Compared with the experimental method, it has become a reliable, cost-effective, and time-saving strategy that can be used to discover new lead compounds. In this study, considering the structural diversity and availability of compounds, the specs database (http://www.specs.net) was first optimized according to the Lipinski and Veber rules, then the specs database was used to virtual screening. Moreover, in order to further understand the binding mode of the HDM2 and the potential modulators, molecular dynamics simulations and MM-PBSA was carried out to calculate the binding free energy[19]. Finally, compounds were tested for tumor inhibition by performing in vitro cell assay.

The workflow for this study which included Activity models, Pharmacophore modeling, Molecular docking, Molecular dynamics (MD) simulations, Binding free energies calculation and cell assay is shown in Fig. 1.

2.1 Data Preparation

For the screening database, the specs database(http://www.specs.net) containing 212,736 compounds were first filtered by Lipinski and Veber rules and 75671 compounds obeyed the Lipinski's cut-off values[20]. Subsequently, the remaining molecules were prepared with LigPrep module provided in the Maestro 12.1 (www.schrodinger.com) using the OPLS_2005 force filed and the protonation states were generated with Epik at PH of 7.4. The prepared database was used to generate 3D conformational set for each molecule using CAESAR in Discovery Studio 4.0 (http://www.accelrys.com)[21] and saved them as Database 1.

For the known ligand database, in order to have a reliable and representative set of active ligands with wide structural diversity, we collected known ligands from literature and patent[11, 22–31]. In this study, compounds whose IC₅₀ values (here reported as pIC₅₀) were determined by the same procedure (fluorescence polarization assay) were considered. Moreover, compounds with defined chirality only were selected in order to have a univocal match between structure and activity. The final the known ligand database was constituted by 48 compounds with pIC₅₀ values in the range 4.0 − 10.0 and sufficient structural heterogeneity. The structures of the collected known ligands are shown in Fig. 2.

For the known ligand database, the molecules in the database need to be properly prepared. Therefore, prepared the structures by washing, adding hydrogen atoms, and energy minimizing in MOE (2019.01). In Database Wash, set the Protonation Field to Dominant (Dominant protonation state at specified pH) and set the Coordinates Field to Rebuild 3D. Other settings with default values. After preparation, the molecules were stored in the mol file with properly protonation. The known ligand data set after preparation was named KLD.

2.2 Preparation of the 3D-Structure of HDM2

Induced-fit theory emphasizes that the ligand can induce the conformational change of the receptor during the binding process. Therefore, it is very important to select a good receptor structure. On this basis, HDM2 protein with a PDB ID of 5OC8 (https://www.rcsb.org/structure/5OC8) [15] was selected as the receptor of this study according to the good resolution. The crystal structure resolution is 1.56 Å and binding of protein HDM2 with small molecule ligand NVP-HDM2, which is a second generation p53–HDM2 inhibitor that had entered phase I clinical trial. Then Discovery Studio 4.0 (DS 4.0) was employed to the prepare protein by adding missing residues, hydrogen atoms as well as deleting water molecules, and then the structure was minimized by using SYBYL-X software (version 2.0) with AMBER99SB force field[34].

2.3 Activity Models

Active models such as QSAR and fingerprint were constructed based on known compounds can be used to filter the initial database.

2.3.1 QSAR Modeling

The QSARs are statistical models that attempt to correlate the structure and experimental properties (activity, selectivity, toxicity, etc.) with computed properties such as molecular descriptors and fingerprints. The QuaSAR suite of applications in MOE is used to build numerical models of the data for prediction and interpretation purposes. For instance, if we get a series of ligands with known activity, we can build a numerical models from the training set to predict the activity of molecules not in the training set.

In this study, the QSAR modeling calculations were performed by using QuaSAR Model in MOE. First, QSAR models used three molecular descriptors to calculate the molecular descriptors for 48 compounds with already known activity in KLD. Then, the AutoQSAR function was employed for the model building using the Partial Least Squares (PLS) regression models. PLS is often used in the main regression analyses[35, 36]. A reasonable QSAR model requires R² > 0.8 and XR² > 0.5. If the model was unreasonable, we need to remove outliers from the database, then rebuilt and cross-validated the model. Finally, the optimized QSAR model was used to filter the database.

2.3.2 Fingerprint Models

A vector of Fingerprint keys contains information on the presence or absence of certain properties or molecular features. According to the design philosophy, different fingerprint schemes emphasize different molecular attributes. In MOE, the fundamental idea of the Fingerprint Model is to encapsulate certain properties directly or indirectly in the fingerprint and then use the fingerprint as a surrogate for the chemical structure. Then comparisons between molecules are reduced to comparing sets of features and measuring the degree of overlapping of the sets.

MOE-Fingerprints can be used for: (i) Fingerprint Clustering to a group of molecules based on the similarity of their fingerprint. (ii) Fingerprint Models to study Quantitative Structure-Activity Relationships using fingerprint descriptions of the molecules. (iii) Similarity Searching to find molecules with similar fingerprints. In this study, we used the fingerprints of a set of active compounds to search databases for similar active compounds by building a fingerprint model.

We built a fingerprint model with KLD mentioned above. The code FP: GpiDAPH3 refers to the Pharmacophore Graph Triangle (GpiDAPH3) fingerprints, and these are three-point pharmacophore fingerprints calculated from the 2D-molecular graph. Therefore, we selected the FP:GpiDAPH3 fingerprint model. Set fingerprint score to Maximum (Maximum of similarities scores) and saved the fingerprint model[37].

2.4 Pharmacophore Modeling

MOE's pharmacophore modeling tools determine the pharmacophoric feature characterizes and their spatial arrangement in 3D that are essential to the binding of ligands to their receptors and the drug activity of the ligands. Pharmacophore models can be generated from the structural data of protein-ligand complexes as well as from ligands when no receptor information is available and from the receptor structure when no ligands are available. The generated models can be used to screen virtual compound libraries for potentially active molecules.

In this study, we used Ligand-based and Complex-based Pharmacophore Modeling for virtual screening.

2.4.1 Ligand-Based Pharmacophore

The common feature pharmacophores were performed by using the pharmacophore Elucidation module in MOE. The HDM2 inhibitor dataset KLD contains 48 known active compounds, was divided into high active compounds(pIC₅₀ ≥ 8) and low active compounds(pIC₅₀ < 8). At first, Pharmacophore Elucidation was employed to identify common features necessary for active HDM2 inhibitors. Set conformation to Clustered to reduce the total number of conformations, the conformation method was set to Conformation Import, the Activity Field was set to pIC50, the value was set to 8, and the other parameters of Pharmacophore Elucidation were set to default values. Then, the selected pharmacophore model with the highest alignment score and manually add ligand shape constraints. Finally, the pharmacophore model was saved for subsequent database screening.

2.4.2 Complex-Based

The Protein-Ligand Interaction Fingerprints (PLIF) tool was used for summarizing the interactions between ligands and proteins using a fingerprint scheme. The interaction between the receptor and the ligand can be observed more intuitively, and the pharmacophore obtained in this way is more reliable. First, selected five crystal structures (PDB ID:4HG7, 4QO4, 4WT2, 4ZYF, 5OC8), which the crystal structure resolution was less than 2.0 Å and binding with good HDM2 inhibitor, imported five HDM2 complex structures into MOE, performed QuickPrep, unchecking the Refine checkbox, using default values for the rest of the settings. The QuickPrep application can correct problems in the protein- ligand complex, adding hydrogen and assigning the protonation states. Tethers were then added to the receptor, heavy atoms and distant atoms are fixed.

The Protein Structure Superposition application was used to superpose the set of 3D protein structures to maximize the spatial overlap of a common atom subset, highlighting the regions of conserved structure and modified regions in the protein collection. It can also be used to overlay and study active sites of related proteins. Finally, the PLIFs were generated and analyzed, and a Pharmacophore Query was generated.

2.5 Molecular Docking

Docking is a computer simulation technique used to predict the binding mode of two interacting chemical species. Dock application from MOE can be used to search for favorable binding modes between ligands and macromolecular targets, which is usually a protein. After the docking is completed, for each ligand, a number of poses are generated and scored. The final highest-scoring poses, along with their scores and conformation energies, were recorded to a database where they are ready for further analysis.

The molecular docking in this study was performed using the crystal structure of HDM2 (PDB ID: 5OC8). First, the structure was prepared using the Protonate 3D module in the MOE to add hydrogens to the entire system and distribute ionization states. The structure was then optimized using energy minimization. Then, the active site of the receptor was checked using SiteView in MOE. The docking site was selected as the ligand-binding site in 5OC8, the initial scoring methodology selected London dG, poses were set to 5, the final scoring methodology selected GBVI/WSA dG, poses were set to 1, while initial settings were adopted for the other conditions. The docking simulation output used score S while considering the binding energy.

Some potential inhibitors with diverse skeletons were selected manually after molecular docking, according to the principle of structural diversity.

2.6 Molecular Dynamics Simulations

Molecular dynamics can dynamically analyze microscopic interactions between compounds from the atomic scale. By analyzing the trajectory of the balanced system, a large amount of kinetic data can be analyzed to predict the macroscopic properties of matter.

The potential inhibitors bound to the HDM2 models were investigated by MD simulations for 100 ns with GROMACS 5.2.0[38]. In this study, the size of the octahedral solvent box matches the size of the protein, and the edge distance is 1 nm. The corresponding Cl⁻ are added to the solvation system to neutralize the system. Before the molecular dynamics simulations, the energy of the protein-ligand complex is minimized, so that the potential energy (Epot) of the system is as small as possible.

During the first stage, a simulation was carried out for 200ps at a temperature of 300K in the Nose-Hoover thermostat (NVT) (constant number of particles, volume, and temperature) ensemble with solute heavy atoms restrained by a force constant of 100 kcal· mol^− 1·Å⁻¹. Then the system was kept to 300K during a 200ps MD simulation in the NPT ensemble, with the same force constant applied to the solute atoms. Finally, in the production phase, the relaxed systems were simulated without restrained under NPT ensemble conditions for 100ns. Within this simulation time, the total energy and RMSD of the protein backbone reached a plateau, and the systems were considered equilibrated and suitable for statistical analysis. The Nose-Hoover thermostat (NVT) was utilized to maintain a constant temperature in all simulations, and the Martyna-Tobias-Klein method was used to control the pressure[34].

2.7 Binding Free Energies Calculated By MM/PBSA

In this study, the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) method[39] was used to estimate the binding free energy between protein HDM2 and potential inhibitors and used the g_mmpbsa tool of GROMACS 5.0.2[19]. In MM/PBSA calculation, the solvent is regarded as a uniform continuous medium, and based on the force field and implicit continuous medium model, the dynamic equilibrium trajectory of the protein-ligand complex is sampled, and the free energies of the complex, protein and ligand are calculated respectively.

2.8 Cell assay

MTS is analog of MTT, the full name is 3-(4,5-dimethylthiazol-2-yl)-5(3-carboxymethoxyphenyl)-2-(4-sulfopheny)-2H-tetrazolium, which is a yellow color dye. Succinate dehydrogenase in the mitochondria of living cells can metabolize and reduce MTS to generate soluble formazan compounds. The content of formazan can be measured at 490nm with a microplate reader. Under normal circumstances, the amount of formazan produced is proportional to the number of viable cells, so the number of viable cells can be inferred from the optical density OD value.

First, inoculate cells, prepare a single cell suspension with RMPI1640 medium containing 10% fetal bovine serum, and inoculate 5000 cells per well into a 96-well plate, with a volume of 100 µl per well, and the cells are inoculated and cultured 12–24 hours in advance. Then, the solution of the compound to be tested was added, the compound was dissolved in DMSO, and the compound was re-screened at the concentrations of 40 µM, 8 µM, 1.6 µM, 0.32 µM, 0.064 µM, and the final volume of each well was 200 µl, and each treatment had 3 replicate wells. After culturing at 37℃ for 48 hours, the adherent cells discarded the medium in the wells, and added 20 µl of MTS solution and 100 µl of culture medium to each well; set up 3 blank duplicate wells, and continue to incubate for 2–4 hours, the light absorption value was measured after allowing the reaction to proceed sufficiently. Finally, the wavelength of 492 nm was selected, and the light absorption value of each well was read with multi-function microplate reader (MULTISKAN FC), and the result was recorded to calculate the IC50 value of the compound.

3.1 Activity Models Filtering

3.1.1 QSAR model generation and validation

In this study, three calculating descriptors (TPSA, SMR, SlogP) and PLS were selected to construct a QSAR model. Then self-cross-validation was performed, and the QSAR model was further analyzed (PLS: R² = 0.289, XR² = 0.174).The cross-validation results were not good. Then removed outliers from KLD, rebuilt and cross-validation the model[40], and got a new QSAR model with 3 descriptors (PLS:R² = 0.825, XR² = 0.752). Finally, the model was evaluated using Model Evaluate in MOE on the database 1. The model predictions were put to a new database field PRED_qsar which represented the predicted value of the model. Selected the compound whose PRED_qsar value is greater than 5.5, screened 29700 compounds, and saved them as Database 2.

3.1.2 Fingerprint models filter

We used the fingerprint model to filter the database 2 and Entered $PRED_ GpiDAPH3_MAX as the Fingerprint prediction field. Selected the compound whose $PRED_ GpiDAPH3_MAX value was greater than 0.5, screened 6729 compounds, and saved them as Database 3.

3.2 Pharmacophore-based virtual screening

3.2.1 Ligand-based pharmacophore model generation and validation

The generation of the Ligand-Based Pharmacophore model was based on KLD. Initially, 68 possible pharmacophore features were discovered by using a feature mapping protocol, and three kinds of common features were included. Then, selected the top ten models of the Alignment Score from the generated pharmacophore model for analyses (Table 1). The coverage of active molecules of the ten models were all 11 compounds. By comparing their separation of actives/inactives and the rationality of the molecular arrangement in the model, model 1(Fig. 3) was selected. The optimized pharmacophore contained three Hydrophobic center (H) and one Aromatic center (R) feature. Finally, added Ligand Shape Volumes (no atoms allowed outside the volume). After assessing the pharmacophore query, virtual screening was carried out E against the database 3. 458 compounds were selected in Database 4.

Table 1

The data fields in the Pharmacophore Elucidation output.
mol	label	cover	overlap	accuracy	P	acc1	acc0	query	Qsize	Q:Aro/PiR	Q:Hyd	Q:Acc2
RHHH_4	R H H H	11	6.8904	0.3333	-0.1898	1.0000	0.1351	4 1 2 2 2	4	1	3	0
RHHa_7	R H H a	11	6.8573	0.3125	-0.1368	1.0000	0.1081	4 1 2 2 4	4	1	2	1
RHHH_2	R H H H	11	6.7462	02708	-0.0498	1.0000	0.0541	4 1 2 2 2	4	1	3	0
RHHa_14	R H H a	11	6.7124	0.3958	-0.3823	1.0000	0.2162	4 1 2 2 4	4	1	2	1
RHHa_9	R H H a	11	6.6524	0.2917	-0.0899	1.0000	0.0811	4 1 2 2 4	4	1	2	1
RHHa_12	R H H a	11	6.5604	0.3542	-0.2484	1.0000	0.1622	4 1 2 2 4	4	1	2	1
RHHa_11	R H H a	11	6.5475	0.3750	-0.3125	1.0000	0.1892	4 1 2 2 4	4	1	2	1
RHHH_3	R H H H	11	6.5424	0.3333	-0.1898	1.0000	0.1351	4 1 2 2 2	4	1	3	0
RHHa_2	R H H a	11	6.5021	0.2917	-0.0899	1.0000	0.0811	4 1 2 2 4	4	1	2	1
RHHa_1	R H H a	11	6.5021	0.2917	-0.0498	1.0000	0.0811	4 1 2 2 4	4	1	2	1

3.2.2 Complex-based pharmacophore model generation and validation

The Protein-Ligand Interactions Fingerprints (PLIF) were generated in MOE, and it was a fingerprint that can characterize protein-ligand interactions. The PLIF display modes of ligands were shown in Fig. 4. We were able to visually analyze the interactions and perform statistical analyses to generate pharmacophore. It was observed that the interaction between ligands and proteins mainly includes the H-bond acceptor and Anionic atom. Therefore, the modification of modify the generated pharmacophore was to include two common features (F1: Acc; F2: Ani&Acc) (Fig. 5). Finally, the generated pharmacophore was used to screen database 4, 456 compounds were obtained, which were saved as Database 5. The pharmacophore was obtained in different ways and compounds were screened in the database so that the result was more reliable.

3.3 Docking-based virtual screening

Redocking was performed for the preprocessed receptor using the 5OC8 ligand. The result showed that an almost similar binding style was identified (RMSD = 0.1369 Å). In this procedure, database 5 was docked with the preprocessed 5OC8. Docking simulation was performed using the initial condition in the application and the scoring function was adopted as GB/VI. After 5 trials, compounds with the highest docking scores (-S) were output and 246 compounds were obtained by docking. Finally, visually inspect the top compound molecules in the docking results, the top 16 compounds with the highest docking scores (-S) and good combination model were selected as Database 6.

3.4 Visual screening

After molecular docking, the top 16 molecules were further analyzed to visually inspect the binding mode of ligands and proteins, and pharmacophores. Finally, 7 compounds (Fig. 6) with structural diversity and appropriate binding modes were selected as potential inhibitors for HDM2 inhibitors. The binding mode of ligands and proteins and pharmacophores are shown in Fig. 7 and Fig. 8.

3.5 Molecular Dynamics Simulation Analyses

In this study, we performed 100 ns MD simulation analysis on the seven Hits and NVP-HDM2 complexes.

The backbone RMSD describes the change of the structure compared with the initial structure, and the simulation can be used to view the stability of convergence and composites. The RMSD of the complexes were shown in Fig. 9. Except Hit2 reaches equilibrium at 40ns, other systems reached equilibrium within the simulation time, and the overall fluctuation was small in the simulation process. The RMSD of all the eight systems showed that the binding of the compound to the receptor was stable. The ligands binding led to a small conformational change in the HDM2 structure from its native one.

The Radius of gyration (Rg) is a structural parameter associated with the overall conformation and three-dimensional (3D) structure of a protein, which has been utilized to get insights into its compactness and folding behavior. The Rg of the complexes were shown in Fig. 10 and the average values of Rg for HDM2 in NVP-HDM2 and potential inhibitors complexes were calculated as 1.325nm, 1.325nm, 1.336nm, 1.321nm, 1.329nm, 1.325nm, 1.320nm, and 1.325nm, respectively. No structural swift was observed in HDM2 after compound binding, where Rg attains a stable equilibrium and thus suggested stability of the complexes throughout the trajectory.

With the aim of determining the mobility and flexibility of each residue in the MD simulation, the Root Mean Square Fluctuation (RMSF) values were calculated. The RMSF of the seven potential active compounds and the original ligands were shown in Fig.11. The fluctuations of the HDM2 backbone were compared with each residue after the ligands binding. The fluctuations were found to be minimum at several key residues of THR49, LEU54, LEU57, TYR60, ILE61, GLY83, VAL93, TYR100. The RMSF suggested thatthe fluctuations of these residues were reduced in the region whereligands were binding.

3.6 MM-PBSA Binding Free Energy Calculations

For the seven-candidate compounds and the reference molecule NVP-HDM2, frames were extracted from 40–60 ns trajectory with stable RMSD value to calculate the MM-PBSA. The dielectric constant is set to 2 in the solute and 80 in the solvent. As shown in Table 2, the binding free energies (ΔG_binding) of NVP-HDM2 and the seven potential inhibitors were − 108.231 kJ/mol, -98.286 kJ/mol, -46.231 kJ/mol,

-41.902 kJ/mol, -120.891 kJ/mol, -95.943 kJ/mol, -59.729 kJ/mol, and − 94.296 kJ/mol. The calculated results of the binding free energy showed that Hit 4 was lower than the reference molecule, which means the binding was more stable with HDM2 than NVP-HDM2, and it shows that Hit4 is more likely to be a potential inhibitor. Except the solvation free energy (∆G_polar), the other energies such as van der Waals energy (∆G_vdw), electrostatic interaction energy (∆G_ele), and non-polar solvation free energy (∆Gnonpola) had a significant contribution in the results of binding free energy. The polar interaction (∆G_ele+∆G_polar) of the eight systems were greater than the nonpolar interactions (∆G_vdw +∆G_nonpolar), indicating that the nonpolar interactions were more conducive to the binding of ligands and HDM2.

Table 2

Binding Free Energies (kJ/mol) of each HDM2 inhibitor complex calculated by MM-PBSA method.
Complexes	Energy(kJ/mol)
Complexes	∆G_vdw	∆G_ele	∆G_polar	∆G_nonpolar	ΔG_binding
NVP-HDM201	-154.554	19.180	82.278	-16.024	-108.231
Hit 1	-132.405	-8.179	58.480	-15.277	-98.286
Hit 2	-75.780	-6.268	43.695	-6.985	-46.231
Hit 3	-101.474	-88.152	159.390	-11.007	-41.902
Hit 4	-174.060	-28.957	98.632	-17.081	-120.891
Hit 5	-144.320	-2.903	65.270	-14.715	-95.943
Hit 6	-90.801	18.703	59.518	-9.565	-59.729
Hit 7	-136.103	-4.274	60.805	-14.372	-94.296

3.7 Cell assay results

Seven Hits were used for cell assay on Hepatoma SMMC-7721, and cisplatin (DDP) and taxol were used as the positive control drug. The experimental results are shown in the Table 3. The results showed that Hit 3 has a strong inhibitory activity on Hepatoma SMMC-7721, and the IC₅₀ value is 14.14µM, which is stronger than the positive control drug DDP(Fig. 12). Other Hits have poor inhibitory effect on Hepatoma SMMC-7721.

Table 3

IC₅₀(µM) of compounds against Hepatoma SMMC-7721.
Compounds	IC₅₀ ± SD(µM)
Hit 1	> 40
Hit 2	> 40
Hit 3	14.14 ± 0.80
Hit 4	> 40
Hit 5	> 40
Hit 6	> 40
Hit 7	> 40
DDP	15.16 ± 1.13
Taxol	0.466 ± 0.025

In this study, we first used Lipinski and Veber rules to filter out compounds that do not meet the properties of a drug. Consequently, we used virtual screening strategy as well as visual inspection to target the HDM2 using three different algorithms (Activity Model, Pharmacophore Model, and Molecular docking of MOE). Finally, seven Hits from the Specs database were selected for further evaluation. In addition, to analyze the binding mode between the complexes, we calculated the molecular dynamics simulation and binding free energies for seven hits and reference molecules. The results showed the binding between the seven Hits and HDM2 was strong, and the seven complex systems were stable. Therefore, seven Hits with novel skeletons can serve as potential candidates for inhibiting HDM2. Cell assay showed that Hit 3 was the most active against SMMC-7721 cell line with an IC₅₀ of 14.14µM and Hit 3 is a potential HDM2 inhibitor

In summary, this virtual screening strategy as well as detailed binding mechanisms of complexes reported herein may provide a new direction for the future discovery and development of potent modulators targeting HDM2.

Acknowledgments

The authors thank the Supercomputing Center of Sichuan University of Science & Engineering Provides Computing Resources.

Funding

This work was supported by the Science and Technology Planning Project of Sichuan Province (Grant 2019YFS0095); the State Key Laboratory of Phytochemistry and Plant Resources in West China (Grant P2017-KF07, P2018-KF14).

Conflict of interest

The authors declare no competing interests.

Ethics approval

Not applicable.

Consent to participate

All authors agreed to contribute to the work, critically drafted the manuscript, and revised it.

Consent for publication

The consent for publication was obtained from all participants.

Availability of data and material

It has been reported and cited in the methodology part of this work.

Code availability

Not applicable.

Author contribution

All persons who meet authorship criteria are listed as authors, and all authors (Cao，Chen，Xiong，Li， Liu and Zhang) certify that they have participated sufficiently in the work to take public responsibility for the content, including participation in the concept, design, analysis, writing, or revision of the manuscript.

Wild WE, Stewart CP (2020) BW, World Cancer Report 2020. IARC Press
Office E (2017) World Cancer Day : Fact sheet and promotion, (2017)
Vazquez A, Bond EE, Levine AJ, Bond GL (2008) The genetics of the p53 pathway, apoptosis and cancer therapy. Nat Rev Drug Discovery 7(12):979–987
Kussie, Paul H, Gorina, Svetlana (1996) Structure of the MDM2 oncoprotein bound to the p53 tumor suppressor transactivation domain,Science
Kubbutat MHG, Jones SN, Vousden KH (1997) Regulation of p53 stability by Mdm2. Nature 387(6630):299–303
Chène P (2003) Inhibiting the p53-MDM2 interaction: an important target for cancer therapy. Nat Rev Cancer 3(2):102–109
Freedman DA, Wu L, Levine AJ (1999) Functions of the MDM2 oncoprotein. Cell Mol Life Sci Cmls 55(1):96–107
Jamil M, David J, Sharon W, Joyce N (1998) The MDM2 gene amplification database. Nucleic Acids Res 26(15):3453–3459
Momand J, Wu HH, Dasgupta G (2000) MDM2–master regulator of the p53 tumor suppressor protein. Gene 242(1–2):15–29
Fry DC (2006) Protein–protein interactions as targets for small molecule drug discovery. Biopolymers 84(6):535–552
Tortorella P, Laghezza A, Durante M, Gomez-Monterrey I, Bertamino A, Campiglia P, Loiodice F, Daniele S, Martini C, Agamennone M (2016) An Effective Virtual Screening Protocol To Identify Promising p53–MDM2 Inhibitors,Journal of Chemical Information & Modeling1216
Vu B, Wovkulich P, Pizzolato G, Lovey A, Ding Q, Jiang N, Liu J-J, Zhao C, Glenn K, Wen Y (2013) Discovery of RG7112: A Small-Molecule MDM2 Inhibitor in Clinical Development. ACS Med Chem Lett 4(5):466–469
Discovery (2013) of RG7388, a potent and selective p53-MDM2 inhibitor in clinical development, Journal of Medicinal Chemistry 56(14) 5979–5983
Shaomeng W, Wei S, Yujun, Zhao D, Isabelle McEachern, Meaux (2014) SAR405838: an optimized inhibitor of MDM2-p53 interaction that induces complete and durable tumor regression,Cancer Research
Dose (2018) and schedule determine distinct molecular mechanisms underlying the efficacy of the p53-MDM2 inhibitor HDM201,Cancer Research
Discovery of AMG (2014) 232, a potent, selective, and orally bioavailable MDM2-p53 inhibitor in clinical development. J Med Chem 57(4):1454–1472
Parks DJ, Lafrance LV, Calvo RR, Milkiewicz KL, Marugán JJ, Raboisson P, Schubert C, Koblish HK, Zhao S, Franks CF (2006) Enhanced pharmacokinetic properties of 1,4-benzodiazepine-2,5-dione antagonists of the HDM2-p53 protein-protein interaction through structure-based drug design. Bioorg Med Chem Lett 16(12):3310–3314
Shangary S, Ding K, Qiu S, Nikolovska-Coleska Z, Bauer JA, Liu M, Wang G, Lu Y, McEachern D, Bernard D, Bradford CR, Carey TE, Wang S (2008) Reactivation of p53 by a specific MDM2 antagonist (MI-43) leads to p21-mediated cell cycle arrest and selective cell death in colon cancer. Mol Cancer Ther 7(6):1533–1542
Kumari R, Kumar R, Lynn A (2014) g_mmpbsa–a GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inform Model 54(7):1951–1962
Li Y, Peng J, Li P, Du H, Li Y, Liu X, Zhang L, Wang L-L, Zuo Z (2019) Identification of potential AMPK activator by pharmacophore modeling, molecular docking and QSAR study. Comput Biol Chem 79:165–176
B BLA, Wei ZA, Sheng GA, B ZZA Discovery of novel modulators targeting human TRPC5: Docking-based virtual screening, molecular dynamics simulation and binding affinity predication,Journal of Molecular Graphics & Modelling102
Zhao Y, Aguilar A, Bernard D, Wang S (2015) Small-Molecule Inhibitors of the MDM2–p53 Protein–Protein Interaction (MDM2 Inhibitors) in Clinical Trials for Cancer Treatment. J Med Chem 58(3):1038–1052
Holzer P, Masuya K, Furet P, Kallen J, Valat-Stachyra T, Ferretti S, Berghausen J, Bouisset-Leonard M, Buschmann N, Pissot-Soldermann C (2015) Discovery of a Dihydroisoquinolinone Derivative (NVP-CGM097): A Highly Potent and Selective MDM2 Inhibitor Undergoing Phase 1 Clinical Trials in p53wt Tumors,Journal of Medicinal Chemistry6348–58
Marugan JJ, Al. E (2006) Enantiomerically Pure 1,4-Benzodiazepine‐2,5‐diones as Hdm2 Antagonists,ChemInform37(37)
Hardcastle IR, Liu J, Valeur E, Watson A, Ahmed SU, Blackburn TJ, Bennaceur K, Clegg W, Drummond C, Endicott JA (2011) Isoindolinone inhibitors of the murine double minute 2 (MDM2)-p53 protein-protein interaction: structure-activity studies leading to improved potency. J Med Chem 54(5):1233–1243
Beck HP, Al E (2012) Improvement of the synthesis and pharmacokinetic properties of chromenotriazolopyrimidine MDM2-p53 protein-protein inhibitors,Cheminform43(11) no-no
Yu S, Qin D, Shangary S, Chen J, Wang G, Ding K, Mceachern D, Qiu S, Nikolovska-Coleska Z, Miller R (2009) Potent and orally active small-molecule inhibitors of the MDM2-p53 interaction. J Med Chem 52(24):7970–7973
Miyazaki M, Naito H, Sugimoto Y, Yoshida K, Kawato H, Okayama T, Shimizu H, Miyazaki M, Kitagawa M, Seki T, Fukutake S, Shiose Y, Aonuma M, Soga T (2013) Synthesis and evaluation of novel orally active p53–MDM2 interaction inhibitors. Bioorg Med Chem 21(14):4319–4331
Novel Inhibitors of the MDM (2014) 2-p53 Interaction Featuring Hydrogen Bond Acceptors as Carboxylic Acid Isosteres. J Med Chem 57(7):2963–2988
Zhang Z, Chu XJ, Liu JJ, Ding Q, Zhang J, Bartkovitz D, Jiang N, Karnachi P, So SS, Tovar C (2013) Discovery of Potent and Orally Active p53-MDM2 Inhibitors RO5353 and RO2468 for Potential Clinical Development,ACS Medicinal Chemistry Letters
Andrea V, Guido, Bold, Alain, De P (2014) Joanna, Tetra-substituted imidazoles as a new class of inhibitors of the p53–MDM2 interaction. Bioorg Med Chem Lett 24(9):2110–2114
Li S, Zhang S, Chen D, Jiang X, Liu B, Zhang H, Rachakunta M, Zuo Z (2020) Identification of Novel TRPC5 Inhibitors by Pharmacophore-Based and Structure-Based Approaches. Comput Biol Chem 87:107302
Jiang X, Li S, Zhang H, Wang L-L (2021) Discovery of potentially biased agonists of mu-opioid receptor (MOR) through molecular docking, pharmacophore modeling, and MD simulation. Comput Biol Chem 90:107405
Liu Y, Lei Y, Guo S, Zuo Z (2020) Ensemble-based virtual screening in discovering potent inhibitors targeting Von Hippel-Lindau (VHL) E3 ubiquitin ligase. Life Sci 262(2):118495
Geladi P, Kowalski BR (1986) Partial Least-Squares Regression: A Tutorial. Anal Chim Acta 185(1):1–17
Siavash Riahi MRG, Moghaddam AB, Pourbasheer E, Norouzi P (2009) Development of a New Combined Chemometrics Method, Applied in the Simultaneous Voltammetric Determination of Cinnamic Acid and 3, 4-Dihydroxy Benzoic Acid,Current Analytical Chemistry5(1) -
Askjaer S, Langgrd M (2008) ChemInform Abstract: Combining Pharmacophore Fingerprints and PLS-Discriminant Analysis for Virtual Screening and SAR Elucidation. J Chem Inf Model 48(3):476–488
GROMACS (2005) : Fast, flexible, and free,Journal of Computational Chemistry
Homeyer N, Gohlke H (2012) Free Energy Calculations by the Molecular Mechanics PoissonBoltzmann Surface Area Method. Mol Inf 31(2):114–122
Hawkins DM, Basak SC, Mills D (2003) Assessing model fit by cross-validation. J Chem Inf Comput Sci 43(2):579–586

PDF

Reviews received at journal
15 Mar, 2022
Reviewers invited by journal
10 Mar, 2022
Editor assigned by journal
07 Mar, 2022
First submitted to journal
01 Mar, 2022

You are reading this latest preprint version

Discovery of Novel Inhibitors Targeting HDM2: Virtual Screening, Molecular Dynamics Simulation, Binding affinity predication and Cell Assay

Status:

Version 1

Abstract

Figures

1. Introduction

2. Material And Methods

2.1 Data Preparation

2.2 Preparation of the 3D-Structure of HDM2

2.3 Activity Models

2.3.1 QSAR Modeling

2.3.2 Fingerprint Models

2.4 Pharmacophore Modeling

2.4.1 Ligand-Based Pharmacophore

2.4.2 Complex-Based

2.5 Molecular Docking

2.6 Molecular Dynamics Simulations

2.7 Binding Free Energies Calculated By MM/PBSA

2.8 Cell assay

3. Results And Discussion

3.1 Activity Models Filtering

3.1.1 QSAR model generation and validation

3.1.2 Fingerprint models filter

3.2 Pharmacophore-based virtual screening

3.2.1 Ligand-based pharmacophore model generation and validation

3.2.2 Complex-based pharmacophore model generation and validation

3.3 Docking-based virtual screening

3.4 Visual screening

3.5 Molecular Dynamics Simulation Analyses

3.6 MM-PBSA Binding Free Energy Calculations

3.7 Cell assay results

4. Conclusion

Declarations

References

Status:

Version 1