Discovery of Benzimidazole Derivatives as Novel Aldosterone Synthase Inhibitors: QSAR, Docking Studies, and Molecular Dynamics Simulation

Citation:  Hong-Mei GUO, Na YU, Le FU, Guang-Ping LI, Mao SHU, Zhi-Hua LIN. Discovery of Benzimidazole Derivatives as Novel Aldosterone Synthase Inhibitors: QSAR, Docking Studies, and Molecular Dynamics Simulation[J]. Chinese Journal of Structural Chemistry, 2022, 41(3): 2203193-2203210.

Discovery of Benzimidazole Derivatives as Novel Aldosterone Synthase Inhibitors: QSAR, Docking Studies, and Molecular Dynamics Simulation

English

• Blood pressure is the pressure of blood circulation on the wall of the body's arteries (the main blood vessels of the body). Hypertension is a chronic disease characterized by persistent high systemic arterial blood pressure (systolic blood pressure ≥ 140 mmHg, diastolic blood pressure ≥ 90 mmHg), which can be accompanied by functional or organic damage of heart, brain, kidney and other organs[1]. Clinically, hypertension can be divided into two categories, primary hypertension and secondary hypertension. About 11.3 billion people around the world are affected by the disease[2, 3]. At the same time, it has also become the "first killer" of Chinese people's health. At present, there are 160 million patients with hypertension in China, so it is also known as "the first disease in China"[4]. About half of the patients in our country show intractable hypertension, and their blood pressure is still high even after receiving a combination of three or more antihypertensive drugs[5]. Patients whose blood pressure is not maintained at normal levels have a significantly increased risk of stroke, heart attack and heart failure[6-8].

At present, the treatment of hypertension mainly depends on the RAAS system. Renin-angiotensin-aldosterone system (RAAS) is a pressor regulation system produced by the kidney in the body, which causes vascular smooth muscle contraction and water and sodium retention, resulting in pressor effect. RAAS inhibitors are widely used in clinical practice, but long-term use of pristine and sartan drugs (more than three months) can cause the phenomenon of "aldosterone escape"[9, 10]. Recent studies have found that aldosterone is not only a deterioration factor of hypertension, but also a risk factor of cardiovascular disease. Aldosterone receptor antagonists include Spironolactone and Eplenone (Fig. 1), but the side effects are hyperkalemia and female breast development. Therefore, there is an urgent need to develop antihypertensive methods for aldosterone receptors[3, 11]. Aldosterone synthase (CYP11B2) is located in the adrenal glomerulus and can reduce the level of aldosterone in the blood[5, 12]. Aldosterone synthase (CYP11B2) is a key enzyme in the synthesis of glucocorticoid aldosterone, which is thought to mediate the last three steps of aldosterone synthesis from 11-deoxycorticosterone to aldosterone[13]. The expression of CYP11B2 in adrenal gland is regulated by renin-angiotensin system (RASS) and plays an important role in maintaining electrolyte metabolism in higher organisms. Abnormal overexpression of CYP11B2 can lead to disruption of mineral balance and may lead to high blood pressure. The effective inhibitor of CYP11B2 should effectively block the biosynthesis of aldosterone, thus reducing the level of plasma aldosterone and finally lowering the blood pressure[14]. The effective treatment strategy for a variety of diseases is to inhibit CYP11B2, which can lead to lower levels of acetone in circulating plasma, including hypertension and heart failure[6, 12]. Considering this principle, inhibition of CYP11B2 may be an attractive way to regulate blood pressure and treat aldosterone-related diseases by reducing plasma aldosterone levels[15].

Figure 1

Figure 1.  Structures of spironolactone and eplenone

For this reason, computer-aided drug design methods have been widely and effectively used, including quantitative structure-activity relationship (QSAR) modeling based on assumptions about compound activity and its structure[16-19]. 3D-QSAR methods, such as comparative molecular similarity index analysis (CoMSIA) and comparative molecular field analysis (CoMFA)[19], are of great significance for the analysis of the relationship between structure and activity of compounds, such as the application of CADD technology in research on hypertension-related targets[20, 21]. A 3D-QSAR model was established based on 40 benzimidazole derivatives to conduct in-depth research on the relationship between their structures and activities[22, 23]. Molecular docking research is applicable to different stages of drug discovery, such as predicting the butt structure of the mating-subject complex and sorting the compound molecules according to the binding. The docking protocol helps clarify the most effective binding posture between the match and the subject[24, 25]. The objective of the current study is tantamount to elucidate the mode of interaction of benzimidazole derivatives with aldosterone synthase (CYP11B2). Therefore, the information provided significant guidance about the design of novel CYP11B2 inhibitors.

40 benzimidazole derivatives of CYP11B2 inhibitors come from the reported literature[22, 23], and their experimental activities in micro-molar and inhibitory activities (–logIC50, pIC50) are listed in Table 1. 40 CYP11B2 inhibitors' 3D structure was constructed and minimized utilizing SYBYL-X 2.1.1 software. With Gasteiger-Huckel charge, Tripos force field and Powll energy gradient method, all molecules were optimized. All parameters were default except the maximum optimization limit and convergence criterion[26]. They were setting as 10,000 times and 0.005 kcal/mol, respectively. We used the structure obtained by the above method as a subsequent 3D-QSAR analysis. Thirty compounds were randomly selected as the training set and the remaining compounds as the test set (marked by the symbol "*")[27]. We selected compound No. 22 (pIC50 = 9.699, Fig. 2) as the template, which has the highest activity. After selecting the common Skeleton (Fig. 2), the superimposed structure of all compounds is shown in Fig. 3.

Figure 2

Figure 2.  Structure of compound No. 22

Figure 3

Figure 3.  Alignment of CYP11B2 inhibitors

SYBYL-X 2.1.1 software was applied in the preparation of CoMFA and CoMSIA models. In partial least squares (PLS) analysis, the CoMFA and CoMSIA descriptors are invoked as independent variables, while the pIC50 value is used as a dependent variable for the development of a 3D-QSAR model. The correlation coefficient (Q2) and the best principal component value (N) of cross validation are determined by leave one method (LOO) for cross validation. We performed non-cross-validation using previously acquired N values to estimate the general determination factor (R2). In addition, the estimated standard error (SEE) and Fischer statistical values (F) were determined[28]. Normally, parameter q2 is greater than 0.5 and the model is acceptable[29, 30]. The value of R2 ranged from 0.5 to 1.0, indicating that the model fit well[31]. The accuracy of the prediction ability of the model is verified by the external validation of the test set[30].

In order to ensure the robustness of the established 3D-QSAR model, the opportunity correlation is avoided, the robustness of the model is controlled, and the Y randomization test is carried out[32]. In this verification method, the active values of the molecules studied (pIC50) are rearranged randomly throughout the learning set and a new model is derived. The values of Q2 and R2 calculated by the new model are lower than those of the original model, and the operation is repeated several times. If the values of Q2 and R2 of the model are low, the appropriate calibration results are not due to chance, but rather to the robustness of the QSAR model.

After the QSAR model was established, it should be validated internally and externally to ensure that the constructed model is robust, reliable, stable, and predictable, and then used to predict the activity of novel compounds[33]. A strictly reliable 3D-QSAR model requires both internal and external validation. The credible 3D-QSAR model also must satisfy the following conditions: the predicted r2 (r2 pred > 0.5) and external standard deviation error of prediction (SDEPext) in Cruciani-Baroni's method[34, 35].

 $\begin{array}{c} r_{pred }^2 = 1 - \frac{{\sum _{i = 1}^{{n_{EXT}}} {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{{\sum _{i = 1}^{{n_{EXT}}} {{{\left( {{{\rm{y}}_i} - {{\bar y}_{TR}}} \right)}^2}} }}\;\;\;\;{\rm{SDEP}} = \sqrt {\frac{{\sum\nolimits_{i = 1}^{{n_{EXT}}} {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{{{n_{EXT}}}}} \\ n_{EXT}: {\rm{Sample\ size\ of\ test\ set}}; {\rm{y}}_{i}: {\rm{Real\ value}}; \widehat{{\mathrm{y}}_{i}} : {\rm{Calculated\ value}} {\rm{Real\ mean\ value\ of\ the\ training\ set}} \end{array}$

Although the QSAR model works very well, it is not helpful to evaluate the suitability of all compounds for the model[27]. Therefore, the model application domain (AD) must be identified when using the constructed QSAR model to screen compounds, and the availability of compounds shall be considered only within the scope of the QSAR model application domain. The application domain methods for evaluating the model mainly include lever distance, residual standard deviation and similarity distance[36]. Here we used the leverage distance method to evaluate the application domain of the model. For a regression QSAR model, a straightforward evaluation of whether a compound is out of the model application domain is its leverage distance to the model, hi[37]:

 $h_{i }= x_{i}^{T}(X^{T}X) ^{–1}x_{i }(i = 1,…………., n)$

Wherein, xi is the descriptor row vector of the predictive compound, the matrix Xn*k–1 is the descriptor variable matrix of the compound model k of n training sets, and the superscript T represents the transpose of matrix/vector. The upper limit of H* is usually set as ±2 or 3K/N, where N is the number of samples in the training set and K is the parameter of the model. The lever distance is higher than the upper limit N* value indicates that the predicted response is the additional inference result of the model, so the forecast compound is not within the application domain of the QSAR model[37].

Density functional theory (DFT) is part of the powerful tools for dealing with the electronic structure and geometry of interacting multi-electron systems[38, 39]. In order to evaluate the details of the electronic structure of the energy levels of small molecular orbitals, and accurately calculate various electronic properties, we calculated DFT. Gauss View 6.0 program and density functional theory were used to study the relevant properties of small molecule compounds. Because B3LYP/6-31G** can reach a good balance point in terms of time and accuracy, the B3LYP/6-31 G** method and basis set are used in the calculation[40]. We introduced the two compounds with the highest activity and the highest predicted activity into the Gauss View 6.0 platform, and calculated HOMO (Highest Occupied Molecular Orbitals), LUMO (Lowest Unoccupied Molecular Orbital) and MESP (Molecular Electrostatic Potential). The MESP, HOMO and LUMO of the molecular frontier orbitals of all optimized structures have been studied by GaussView 6.0 program[41]. The electrostatic potential V(r) at the r point of the molecular system, the nuclear charge (ZA) at (RA) and the electronic density ρ(r) are used to test the molecules with functional points as follows[42]:

 $V(r) = \sum _{A = 1}^N {\frac{{{Z_A}}}{{\left| {r - {R_A}} \right|}}} - \int {\frac{{\rho \left( {{r^\prime }} \right){d^3}{r^\prime }}}{{\left| {r - {r^\prime }} \right|}}}$

Where N represents the total number of nucleus in a molecule, $\frac{{Z}_{A}}{\left|r-{R}_{A}\right|}$ is the naked nuclear potential contributed by the nucleus, and $\frac{{\rho \left( {{r^\prime }} \right){d^3}{r^\prime }}}{{\left| {r - {r^\prime }} \right|}}$occurs due to a constant electron charge density.

The complex crystal structure of aldosterone synthase and small molecules comes from the PDB database (https://www1.rcsb.org/, PDB ID: 4DVQ, 2.49 Å resolution). The docking of CYP11B2 inhibitor and aldosterone synthase (4DVQ) was simulated by the Surflex-Dock module in SYBYL-X 2.1.1 software. Before docking, the protein needs to be pre-treated. We used SYBYL-X 2.1.1 software to preprocess 4DVQ, including processing the ends of the protein backbone and repairing the missing side chains, as well as adding hydrogen atoms, removing water molecules, and adding charges to 4DVQ[30]. Next, we constructed a pair of interface bag based on 4DVQ's exhibiting small molecule ligands, as shown in Fig. 4. The original molecule was re-docked to binding site, and RMSD between the redocking molecule and the original ligand was analyzed to verify the reliability of the docking method. After docking, the Surflex-Dock module was utilized to score the interactions between the receptor and molecules. Finally, binding posture with the highest inhibitory activity was selected and the corresponding complex was output and analyzed in Discovery Studio 3.0 software and PyMOL software (https://pymol.org/2/).

Figure 4

Figure 4.  Active pockets produced after protein (4DVQ) pretreatment. The yellow stick represents the original ligand (1CA), and magentas stick represents HEM

Molecular dynamics simulation was executed with the AMBER16 package[43, 44]. Receptor loaded the Amber ff14SB force field[45] and the ligand loaded the general Amber force field (GAFF)[46], and the whole complex system was solvated in TIP3PBOX (Buffer ≥ 10.0 Å) and then the system was electrically neutral by adding sodium and chloride ions.

We optimized the system to avoid unreasonable structural areas in the system, and then gradually heat the system from 0 to 300 K for the heating phase. Finally, the temperature was kept at 300 K in the next stage. The entire MD process took 2 fs as the time step. Periodic boundary conditions were applied to hold constant temperature and pressure. The Langevin dynamics method was used to adjust the temperature at a collision frequency of 2 ps-1 and the pressure was set on 1 atm under each anisotropic pressure scale protocol. The Particle Mesh Ewald (PME) method was employed to deal with long-range electrostatics and the range of actual cutoff spatial interaction was less than 1 nm. All covalent bonds involving hydrogen atoms were restricted by the SHAKE method. Next, 100 ns MD simulation was performed, and the trajectory was saved every 2 ps.

The MM/GBSA algorithm was used to process the saved MD simulation trajectories and calculate the binding energy with crystal complexes of different ligands[47, 48]. A total of 100 frames were collected from the last 90 to 100 ns simulation to calculate the average binding energy. Formula is as follows:

 $ΔG_{Bind} = G_{complex} – (G_{protein} + G_{ligand}) = ΔG_{sol} + ΔG_{gas}$
 $ΔG_{sol} = ΔE_{GB} + ΔE_{SURF}; ΔG_{gas} = ΔE_{ele} + ΔE_{vdw}$

Among them, ΔGBind is the binding free energy, Gcomplex, Gprotein, and Gligand are related free energies, and ΔGsol represents the total mechanical energy of molecules in vacuum, which can be further divided into electrostatic contribution (ΔEele) and van der Waals force (ΔEvdw). The term can be calculated using molecular mechanics methods. ΔESOL is the solvation energy, including the polar solvation energy (ΔEGB) calculated by the generalized birth (GB) approximation model and the non-polar part (ΔESURF) obtained by fitting the solvent accessible surface area (SASA) and paired linear combination overlap (LCPO) model. In addition, the energy of each residue is broken down into the main chain and side chain atoms. The energy breakdown can be analyzed to determine the contribution of key residues to the binding.

The three-dimensional QSAR model based on atoms was successfully established by using the partial least-two-multiplication regression analysis, and iso-potential maps were used to analyze the favorable and unfavorable regions, respectively. The visualization of the equipotential maps was helpful for us to analyze which group changes can increase the activity of the compound used to design new compounds. The results of the 3D-QSAR model are presented in Table 2, which has achieved an excellent correlation coefficient (R2 > 0.9) of 0.992 and 0.994 with standard error of estimate (SEE) of 0.07, and a predicted coefficient (Q2 > 0.5) of 0.877 and 0.848, with Fisher values of 223.911 and 343.945. The predicted model summary of statistical data is listed in Table 1. The scatter plot of the XY axis of the correlation between actual pIC50 and predicted pIC50 is represented in Fig. 5(a, b) for the CoMFA and CoMSIA models. The external validation results are shown in Table 2, and we regard the 3D-QSAR model as responsible with good statistical significance.

Figure 5

Figure 5.  pIC50 of predicted versus actual activity of the training and test sets. (A) CoMFA (B) CoMSIA

The COMFA model is exhibited in Figs. 6a and 6b. In the steric field (proportion of 44.6%), the green contours present the contribution to a high steric tendency, while yellow contours represent low steric effect. In Fig. 6a, there is a larger volume of yellow equipotential region at the position of the R3 group, which indicates that the volume of the group decreased here will be beneficial to enhance the activity. In the electrostatic field, the red contour indicates that it is advantageous to increase the negative charge group, while the blue one shows that the increase in the positive charge group is favorable. In Fig. 6b, the red contour around the R1 group exposing a negatively charged substituent at this position might improve the activity and the entire compound is mainly surrounded by red patches. The CoMFA model suggests a little group of the R3 substituent and it should be considered as negatively charged groups. CoMSIA models contour maps of hydrophobic and H-bond acceptor fields are shown in Fig. 6c and 6d. COMSIA contour maps at electrostatic and stereoscopic fields of compound No. 22 showed similar results from COMFA model, and the effect of hydrogen bond to the body field is much smaller than that of hydrophobic field and H-bonded subject field, so only the latter two fields are discussed in this paper. In Fig. 6c, the white figure takes up more than half of the compound, indicating that the compound should have hydrophilic groups to increase its inhibitory activity, especially the R3 group. In Fig. 6d, R3 substituent indicates that the addition of hydrogen bond donor here is not conducive to the improvement of activity. All in all, the introduction of small-volume groups near the R3 substituent, or the introduction to negatively charged groups near R1 replacement bases or the groups capable of accepting hydrogen bonds or the hydrophilic groups near R3, can help to increase the activity of inhibitors.

Figure 6

Figure 6.  Three-dimensional equipotential maps of CoMFA(a, b) and CoMSIA(c, d) with No. 22. (a) stereo, (b) electrostatic, (c) hydrophobic and (d) hydrogen bond acceptor

The robustness of the 3D-QSAR model was verified by the Y-random method, avoiding the chance correlation. In the current case, 60 random tries were executed for the final developed 3D-QSAR models. Y-randomization validation produces models that do not correspond to the initial model. The Q2 and R2 values of the new QSAR model are indeed very low. Therefore, the possibility of random correlations was ruled out. The randomly arranged bioactivities were used for the test. All the results are summarized in the Table 3.

Table 3

The applicability domain (AD) obtained by the leverage distance method is the Williams diagram (Fig. 7), where the standardized residuals (σ ± 2 or 3) and the leverage threshold (h* = 0.146/0.366) are plotted. As Fig. 7b shows, it is obvious that there is one compound outlier of the test set; the chemical is identified as an outlier for the COMSIA model. Other compounds have a standard deviation in the out of ±2 or 3. This mistaken prediction might perhaps be accredited to wrong experimental data or the structures of these compounds.

Figure 7

Figure 7.  Williams plot for CoMFA(a) and COMSIA(b) models

(The leverage threshold h* = 0.146/0.366 and standardized residual σ = ±2 or 3)

The model compound (No. 22) and the new derivative (No. a2) were used for orbital energy level and electron cloud density analysis. The results show these two compounds have unique orbital energy levels and electron cloud properties. In Fig. 8, the MESP diagram shows the MESP curves of the two compounds, in which the negatively charged potential is represented in red and the positively charged region in blue. In drug design, improving the complementarity of electrostatic potential can improve the affinity and selectivity of ligands; and the dominant compounds can be judged by comparing the MESP of different molecules. From these results, obviously the positive and negative charge distributions on the surface of the model compound and the newly designed compound (No. a2) are more uniform and complementary. In addition, it also verifies the credibility of the 3D electrostatic equipotential map constructed by the QSAR model. When the molecules interact in the form of charge transfer complex, the HOMO energy can be used as a measure of the electron-giving ability of the molecule, while the LUMO energy of the electron-receiving ability, that is, the electron transfer from the HOMO of the donor to the LUMO of the acceptor. The results of No. 22 and No. a2, the HOMO's outcome ranges from –5.660 to –6.011 and LUMO –0.802 to –1.084, and frontier orbital energy difference (HOMO-LUMO gap, HLG) of –9.337 and –15.656 are tabularized in Table 4. The energy gap of the compound is the energy difference (HOMO-LUMO) between HOMO and LUMO. It is an important parameter for the excited transition of reactants. The energy levels and electron cloud density distribution of HOMO and LUMO of each molecule were analyzed (Fig. 8). In compounds No. 22 and No. a2, most of the HOMO orbital electrons appear on the common skeleton benzimidazole ring, while the LUMO orbital electrons are concentrated in other parts (except the benzimidazole ring). Therefore, for compounds No.22 and No. a2, the electron transitions are transferred from the benzimidazole ring to other groups (except the imidazole ring), indicating an obvious charge transfer process, which indicates that the benzimidazole ring has strong electrical absorptivity, therefore resulting in an increase in the electron cloud density of the compound as a whole.

Figure 8

Figure 8.  MESP superimposed onto a surface of constant electron density, HOMO and LUMO

Table 4

In the present studies, we employed molecular docking calculations for compound with the highest biological activity (No. 22) and that with the lowest biological activity (No. 04), respectively. Before docking, we calculated the root-mean-square deviation of the redocked conformation of original ligand (1CA, Fig. 9) and the value of RMSD is 0.76Å (< 2.0 Å), which indicated the generated docking scheme is reliable and can be used for subsequent research[49]. The summarized results of docking calculations are presented in Table 5. The highest biological activity (No. 22, pIC50 = 9.699) showed promising docking score (Docking Score = 7.581) surrounded by residues Arg120, Phe130, Ala313, Gly314, Val378, Cys450, Phe487 and HEM. However, compound No. 04 with the lowest bioactivity (pIC50 = 6.217) had only 3.8889 docking scores near Arg110, Phe130, Val378, Gly379, Phe381 and HEM residues. The binding mode and interaction between compound No.22/No.04 and aldosterone synthase are shown in Fig. 10.

Figure 9

Figure 9.  Overlay of the original ligand docked (1CA, shown as cyan sticks) and re-docked (shown as wheat sticks) in the aldosterone synthase active site

Figure 10

Figure 10.  Interaction between protein and molecule. (a) No. 22, (b) No. 04. The small molecule ligand is shown as green sticks and HEM as magentas sticks. π-π interactions are represented by pink dashed lines and hydrogen bond is the yellow dashed lines

As shown in Fig. 10b, compound No. 04 was inserted into the binding cage surrounded by Arg110, Phe130, Val378, Gly379, Phe381 and HEM. The 6-position fluorine atom on the benzimidazole ring of the nucleus forms a hydrogen bond with the guanidine group NH of Arg110, while the nitrogen atom on the R3 substituent makes a H-bond with the amino group of Gly379. It is also seen that the 6-position nitrogen atom on the benzimidazole ring forms a conventional hydrogen bond with Gly379. At the same time, an unfavourable bump is formed between the carbon atom on the benzimida-zole loop and the amino acid residue Phe381. The aromatic ring of compound No. 04 forms a stable π-π stacking interaction of Phe130, and forms π-π shaped interaction and π-alkyl interaction with HEM, including benzimidazole ring and pyridyl. In contrast to compound No. 22 in Fig. 10a, the 6-position fluorine atom makes hydrogen bond with the NH of Arg120. And benzimidazole loop produces amide-π stacking interaction with Ala313 and Gly314, while two methyl groups at pyridyl make π-alkyl and alkyl interactions with Val378 and Phe487, respectively. The cyclopropyl of compound No. 22 forms a hydrophobic interaction with the aromatic ring of Phe130. The difference is perhaps that the R3 substituent forms π-sulfur interaction with Cys450, and only cyclopropyl and HEM produce weaker π-alkyl and alkyl interactions. It is possible that the weaker interaction with HEM and the difference in amino acids caused the two compounds to have such diverse activities. This also explains why hydrophobic interactions can enhance the potency. At the same time, this means that the hydrophobic interaction between the inhibitor and the protein can make it integrate well into the human receptor and form the most durable energy for the drug receptor[50, 51].

To begin our investigation into novel CYP11B2 inhibitors featuring a benzimidazole template, we modified some groups of compound No. 22 with various substituent groups while keeping the benzimidazole ring structure according to the information obtained from molecular docking and QSAR model. Table 6 lists the structure and predicted pIC50 value of the novel CYP11B2 inhibitor compounds. The results show that the five newly compounds designed have better pIC50 values than their model, of which No. a2 is most likely to become the lead compound, and the result of molecular docking further improves the reliability of the obtained compound.

Table 6

ADME/T, including absorption, distribution, metabolism, elimination and toxicity, is a critical parameter that is commonly used in clinical trials and the selection of the development of drugs[52]. We used the online tool SwissADME (http://www.swissadme.ch/index.php) to evaluate the synthe-tic availability of new compounds[53]. For uploaded compounds, the web server can predict the percentage of absorption of the drug in the human intestine. Radar maps of compounds that indicate drug similarity are important for identifying compounds with poor absorption and permeability and for addressing parts of the barrier range in which compounds must demonstrate their drug ability. To accurately examine molecules, we used drug sample indices to compare Lipinski 5 criteria for the selection of oral bio-use drugs (Table 7). Each designed compound complies with Lipinski rules. In addition to the five rules, the MW independent rule accurately predicts the oral bio-utilization of compounds. According to the radar map (Fig. 11), all compounds show good drug similarity.

Figure 11

Figure 11.  Radar maps of compound No. 22 and its derivative

To further verify whether newly designed inhibitors can bind well to aldosterone synthase, we ran 100 ns MD simulation on the template compound No. 22 and the most likely candidate compound molecule No. a2. The root mean square deviation (RMSD) and the root mean square fluctuation (RMSF) of the complex, the ligand and the binding pocket (defined as residues within 5 Å around the ligand) were calculated to research the stability of the structure and residues in the system. Fig. 12 shows RMSD and RMSF after MD simulation. As shown in Fig. 12a, the RMSD of the protein-compound No. 22/No. a2 complexes fluctuated around ~1.75 Å after 40 ns simulation. Although there is a large fluctuation in compound No. a2 (Fig. 12a), the complex of No. a2 and protein also stabilized at ~1.75 Å after 50 ns. The RMSD of two complexes fluctuated within a range of less than 2.75 Å in the late simulation stage. The fluctuating trend of amino acid residues measured by calculating the RMSF parameters is given in Fig. 12b, indicating that the position changes of the residues are relatively stable with the extension of the simulation time. Taking the template compound as a reference, motion intensity of the newly designed molecule tends to be stable and the fluctuation is small. This phenomenon indicates that the binding of the novel molecule to the protein is more stable.

Figure 12

Figure 12.  RMSD and RMSF of complexes docked with compound No. 22/a2 for 100 ns simulations

Then as showed in Figs. 13 and 14, we compared the interaction before and after MD simulation. Compounds No. 22 and No. a2 both form a hydrophobic interaction between similar amino acid residues, which are marked with a dotted line. Figs. 13 and 14 show the protein-ligand interaction diagrams of compound No.22/No.a2 before and after molecular dynamics. It can be seen from the figures that the amino acid residues of the interaction between the small molecule ligand and the protein are roughly similar, which also confirms the result of molecular docking, that is, the interaction between the aldosterone synthase inhibitor and the protein is mainly based on the hydrophobic interaction between the molecules. Compound No. 22 is located in the cleft between protein and HEM (Fig. 13), in a pocket surrounded by residues Arg110, Phe130, Ala313, Phe381, Leu382, Glu383, and Phe487 after 100 ns of molecular dynamics. As shown in Fig. 13b, the difference lies in that compound No. 22 and Arg110 form π-cation interaction instead of a hydrogen bond. The 6-position fluorine atom forms a hydrogen bond with Glu383, while the methyl group on the pyridyl group only generates alkyl interaction with HEM.

Figure 13

Figure 13.  Interaction between 4DVQ and compound No. 22 before (A) and after (B) 100ns MD simulation

Figure 14

Figure 14.  Interaction between 4DVQ and compound No. a2 before (A) and after (B) 100 ns MD simulation

The protein-ligand interaction diagram of the new compound No. a2 designed based on the 3D-QSAR model and the results of molecular docking before and after molecular dynamics are shown in Fig. 14, a diagram of the interaction between compound No. a2 and protein. Fig. 14 depicts that compound No. a2 is surrounded by amino acid residues Trp116, Phe130, Met238, Phe381, Cys450, Phe487 and Ile488, and generates a hydrophobic bond of different lengths with each amino acid residue. The molecular docking effect of the new derivative compound No. a2 is better than that of the template compound No. 22 (Table 5). The protein-ligand interaction diagrams in Figs. 13 and 14 show that the former has more important amino acid residues than the template compound, such as Ile488. Moreover, Ile488 can also form π-π hydrophobic interactions with small molecule ligands. This also explains the importance of the amino acid residues Phe130, Ala313, Phe487 and Ile488.

After MD simulation, the MM/GBSA method was used to decompose the binding free energy to find the key residues in the binding process of 4DVQ and inhibitors. The binding free energy and various energy items in the MD process are shown in Fig. 15 and Table 8. Compared with compound No. 22, compound No. a2 has better binding energy, forms key residues, and has more amino acids, indicating that it may have better activity, which is consistent with its larger molecular docking score.

Figure 15

Figure 15.  Contributions of free energy calculated by MM/GBSA for key residues of 4DVQ (A) No. 22 and (B) No. a2

Table 8

To find novel and highly effective anti-hypertension drugs CYP11B2 inhibitors, based on the CoMFA and CoMSIA methods we established a 3D-QSAR model by using molecular docking. The interaction between compounds and proteins was analyzed, and the activities of 5 novel CYP11B2 inhibitors were designed and predicted. From 3D-QSAR model's contour map analysis, the contribution value of the stereo field of the molecule is larger than that of the electrostatic field, and the effects of hydrophobic and hydrogen bond receptor fields on inhibitor molecular activity cannot be ignored. The introduction to small-volume groups near the R3 substituent or the introduction to negatively charged groups near R1 replacement bases or the groups capable of accepting hydrogen bonds or the hydrophilic groups near R3, can help to increase the activity of CYP11B2 inhibitor compounds. In addition, we also found that the interaction mode of compounds and proteins is mainly hydrogen bonding interaction. Some amino acids played an important role in combining protein and inhibitor, such as Phe130, Ala313 and Phe487. The predicted activity values of the newly designed compounds were higher than the respective templates. In addition, these compounds have shown beneficial results of the synthesis feasibility and ADMET evaluation. In summary, through construction, verification and analysis of the 3D-QSAR model, combined with the analysis of the mechanism of molecular docking, new ideas and directions for the subsequent development of novel anti-hypertension drugs are provided.

1. [1]

Majid, N.; Kemal, P. Automatic classification of hypertension types based on personal features by machine learning algorithms. Math. Probl. Eng. 2020, 28, 1-13.

2. [2]

Mills, K. T.; Bundy, J. D.; Kelly, T. N.; Reed, J. E.; Kearney, P. M.; Reynolds, K.; Chen, J.; He, J. Global disparities of hypertension prevalence and control: a systematic analysis of population-based studies from 90 countries. Circulation 2016, 6, 441-450.

3. [3]

Zhou, B.; Bentham, J.; Di, C. M.; Bixby, H.; Danaei, G.; Cowan, M. J.; Paciorek, C. J.; Singh, G.; Hajifathalian, K.; Bennett, J. E.; Taddei, C.; Bilano, V.; Carrillo-Larco, R. M.; Djalalinia, S.; Khatibzadeh, S.; Lugero, C.; Peykari, N.; Zhang, W. Z.; Lu, Y.; Stevens, G. A.; Riley, L. M.; Bovet, P.; Elliott, P.; Gu, D. F.; Ikeda, N.; Jackson, R. T.; Joffres, M.; Kengne, A. P.; Laatikainen, T.; Lam, T. H.; Laxmaiah, A.; Liu, J.; Miranda, J. J.; Mondo, C. K.; Neuhauser, H. K.; Sundstrom, J.; Smeeth, L.; Soric, M.; Woodward, M.; Ezzati, M. Worldwide trends in blood pressure from 1975 to 2015: a pooled analysis of 1479 population-based measurement studies with 19·1 million participants. Lancet 2017, 10064, 37-55.

4. [4]

Ryo, S.; Wataru, S.; Yuichi, O.; Minami, Y.; Hideki, U.; Yuki, H.; Kanako, S.; Masashi, N.; Yasuhiro, E.; Kei, T.; Hidetoshi, S.; Tomoko, O.; Fumihiko, A. Discovery of novel pyrazole-based selective aldosterone synthase (CYP11B2) inhibitors: a new template to coordinate the heme-iron motif of CYP11B2. J. Med. Chem. 2018, 13, 5594-5608.

5. [5]

Tomaschitz, A.; Pilz, S.; Ritz, E.; Obermayer-Pietsch, B.; Pieber, T. R. Aldosterone and arterial hypertension. Nat. Rev. Endocrinol. 2009, 2, 83-93.

6. [6]

Palo, K.; Barone, N. J. Hypertension and heart failure: prevention, targets, and treatment. Heart. Fail. Clin. 2020, 1, 99-106.

7. [7]

Lehmann, L. H.; Jebessa, Z. H.; Kreusser, M. M.; Horsch, A.; He, T.; Kronlage, M. A proteolytic fragment of histone deacetylase 4 protects the heart from failure by regulating the hexosamine biosynthetic pathway. Nat. Med. 2017, 1, 62-72.

8. [8]

Petrea, R. E.; O'Donnell, A.; Beiser, A. S.; Habes, M.; Romero, J. R. Mid to late life hypertension trends and cerebral small vessel disease in the framingham heart study. Hypertension 2020, 3, 707-714.

9. [9]

Struthers, A. D.; Macdonald, T. M. Review of aldosterone-and angiotensin ii-induced target organ damage and prevention. Cardiovasc. Res. 2004, 4, 663-670.

10. [10]

Yang, T.; Chen, Y. Y.; Liu, J. R.; Zhao, H.; Zhao, Y. Y. Natural products against renin-angiotensin system for antifibrosis therapy. Eur. J. Med. Chem. 2019, 179, 623-633.

11. [11]

Flynn, J. T.; Kaelber, D. C.; Baker-Smith, C. M.; Blowey, D.; Carroll, A. E.; Daniels, S. R.; Ferranti, S. D.; Dionne, J. M.; Falkner, B.; Flinn, Susan K.; Gidding, S. S.; Goodwin, C.; Leu, M. G.; Powers, M. E.; Rea, C.; Samuels, J.; Simasek, M.; Thaker, V. V.; Urbina, E. M. Clinical practice guideline for screening and management of high blood pressure in children and adolescents. Pediatrics 2017, 3, e20171904-e20171978.

12. [12]

Papillon, J.; Lou, C.; Singh, A. K.; Adams, C. M.; Watson, C. Discovery of n-[5-(6-chloro-3-cyano-1-methyl-1h-indol-2-yl)-pyridin-3-ylmethyl]-ethanesulfonamide, a cortisol-sparing CYP11B2 inhibitor that lowers aldosterone in human subjects. J. Med. Chem. 2015, 23, 9382-9394.

13. [13]

Wang, J. G.; Li, Y. Characteristics of hypertension in Chinese and their relevance for the choice of antihypertensive drugs. Diabetes 2013, s2, 67-72.

14. [14]

Ito, R.; Sato, I.; Tsujita, T.; Yokoyama, A.; Sugawara, A. An ubiquitin-proteasome inhibitor bortezomib suppresses the expression of CYP11B2, a key enzyme of aldosterone synthesis. Biochem. Biophys. Res. Commun. 2017, 1, 21-28.

15. [15]

Cui, C. X.; Chen, H.; Li, S. J.; Zhang, T.; Lan, Y. Mechanism of ir-catalyzed hydrogenation: a theoretical view. Coord. Chem. Rev. 2020, 412, 213251-213272.

16. [16]

Ban, F.; Dalal, K.; Li, H.; Leblanc, E.; Rennie, P. S.; Cherkasov, A. Best practices of computer-aided drug discovery: lessons learned from the development of a preclinical candidate for prostate cancer with a new mechanism of action. J. Chem. Inf. Model. 2017, 5, 1018-1028.

17. [17]

Kholod, Y.; Hoag, E.; Muratore, K.; Kosenkov, D. Computer-aided drug discovery: molecular docking of diminazene ligands to DNA minor groove. J. Med. Chem. 2018, 5, 882-887.

18. [18]

Shekhar, C. In silico pharmacology: computer-aided methods could transform drug development. Chem. Biol. 2008, 5, 413-414.

19. [19]

Tantillo, D. J.; Siegel, J. B.; Saunders, C. M.; Palazzo, T. A.; Painter, P. P.; Brien, T. E.; Nuñez, N. N.; Nouri, D. H.; Lodewyk, M. W.; Hudson, B. M.; Hare, S. R.; Davis, R. L. Computer-aided drug design for undergraduates. J. Chem. Educ. 2019, 5, 920-925.

20. [20]

Morales, B. A.; Matute, R. A.; Caballero, J. Understanding the comparative molecular field analysis (CoMFA) in terms of molecular quantum similarity and DFT-based reactivity descriptors. J. Mol. Model. 2015, 6, 156-166.

21. [21]

Zhang, C.; Feng, L. J.; Huang, Y.; Wu, D.; Li, Z.; Zhou, Q.; Wu, Y.; Luo, H. B. Discovery of novel phosphodiesterase-2a inhibitors by structure-based virtual screening, structural optimization, and bioassay. J. Chem. Inf. Model. 2017, 2, 355-364.

22. [22]

Hoyt, S. B.; Park, M. K.; London, C.; Xiong, Y.; Ali, A. Discovery of benzimidazole CYP11B2 inhibitors with in vivo activity in rhesus monkeys. ACS. Med. Chem. Lett. 2015, 5, 573-578.

23. [23]

Sparks, S. M.; Danger, D. P.; Hoekstra, W. J.; Leesnitzer, T.; Becherer, J. D. Development of highly selective pyrimidine-based aldosterone synthase (CYP11B2) inhibitors. ACS. Med. Chem. Lett. 2019, 7, 1056-1060.

24. [24]

Rarey, M.; Kramer, B.; Lengauer, T.; Klebe, G. A fast flexible docking method using an incremental construction algorithm. J. Mol. Biol. 1996, 3, 470-489.

25. [25]

Muthiah, I.; Rajendran, K.; Dhanaraj, P. In silico molecular docking and physicochemical property studies on effective phytochemicals targeting GPR116 for breast cancer treatment. Mol. Cell. Biochem. 2021, 476, 883-896.

26. [26]

Clark, M.; Iii, R.; Opdenbosch, N. V. Validation of the general purpose tripos 5.2 force field. J. Comput. Chem. 2010, 8, 982-1012.

27. [27]

Gagic, Z.; Nikolic, K.; Ivkovic, B.; Filipic, S.; Agbaba, D. QSAR studies and design of new analogs of vitamin E with enhanced antiproliferative activity on MCF-7 breast cancer cells. J. Taiwan Inst. Chem. E 2016, 33-44.

28. [28]

Bai, Y. B.; Gao, Y. Q.; Nie, X. D.; Tuong, T.; Gao, J. M. Antifungal activity of griseofulvin derivatives against phytopathogenic fungi in vitro, in vivo, and 3D-QSAR analysis. J. Agric. Food Chem. 2019, 22, 6125-6132.

29. [29]

Li, Y. F.; Chang, Y. Q.; Deng, J.; Li, W. X.; Jian, J.; Gao, J. S.; Wan, X.; Gao, H.; Kurihara, H.; Sun, P. H.; He, R. R. Prediction and evaluation of the lipase inhibitory activities of tea polyphenols with 3D-QSAR models. Sci. Rep. 2016, 6, 34387-34387.

30. [30]

Tian, C. W.; Li, P. C.; Xin, Y. H.; Lei, Z.; Wan, P. Identification of potential tubulin polymerization inhibitors by 3D-QSAR, molecular docking and molecular dynamics. RSC. Adv. 2017, 7, 38479-38489.

31. [31]

Li, R.; Du, Y.; Gao, Z.; Shen, J. Molecular modeling studies on carbazole carboxamide based BTK inhibitors using docking and structure-based 3D-QSAR. Int. J. Mol. Sci. 2018, 4, 1244-1268.

32. [32]

Khamouli, S.; Belaidi, S.; Ouassaf, M.; Lanez, T.; Chtita, S. Multi-combined 3D-QSAR, docking molecular and ADMET prediction of 5-azaindazole derivatives as LRRK2 tyrosine kinase inhibitors. J. Biomol. Struct. Dyn. 2020, 1-14.

33. [33]

Lei, B.; Li, J.; Lu, J.; Du, J.; Liu, H.; Yao, X. Rational prediction of the herbicidal activities of novel protoporphyrinogen oxidase inhibitors by quantitative structure-activity relationship model based on docking-guided active conformation. J. Agric. Food Chem. 2009, 9593-9598.

34. [34]

Cruciani, G.; Baroni, T.; Clementi, S.; Costantino, F.; Skagerberg, B. Predictive ability of regression models. Part i: standard deviation of prediction errors (SDEP). J. Chemom. 1989, 3, 499-509.

35. [35]

Tosco, P.; Balle, T. A 3D-QSAR-driven approach to binding mode and affinity prediction. J. Chem. Inf. Model. 2012, 2, 302-307.

36. [36]

Tropsha, A. Alexander Golbraikh. Predictive QSAR modeling workflow, model applicability domains, and virtual screening. Curr. Pharm. Des. 2007, 34, 3494-3504.

37. [37]

Tunkel, J.; Mayo, K.; Austin, C.; Hickerson, A.; Howard, P. Practical considerations on the use of predictive models for regulatory purposes. Environ. Sci. Technol. 2005, 7, 2188-2199.

38. [38]

Eroglu, E.; Türkmen, H. A DFT-based quantum theoretic QSAR study of aromatic and heterocyclic sulfonamides as carbonic anhydrase inhibitors against isozyme, CA-II. J. Mol. Graph. Model. 2007, 4, 701-708.

39. [39]

Maciej, S.; Małgorzata, W.; Ryszard, T.; Jakub, G. Application of artificial neural networks and DFT-based parameters for prediction of reaction kinetics of ethylbenzene dehydrogenase. J. Comput. Aided Mol. Des. 2006, 20, 145-157.

40. [40]

Bai, T. W.; Ni, X. F.; Ling, J.; Shen, Z. Q. Mechanism of Janus polymerization: a DFT study. Chin. J. Polym. Sci. 2019, 10, 94-100.

41. [41]

Frisch, G. W.; Trucks, H. B.; Schlegel, G. E.; Scuseria, M. A.; Robb, J. R.; Cheeseman, G.; Scalmani, V.; Barone, G. A.; Petersson, H.; Nakatsuji, X.; Li, M.; Caricato, A.; Marenich, J.; Bloino, B. G.; Janesko, R.; Gomperts, B.; Mennucci, H. P.; Hratchian, J. V.; Ortiz, A. F.; Izmaylov, J. L.; Sonnenberg, D.; Williams-Young, F.; Ding, F.; Lipparini, F.; Egidi, J.; Goings, B.; Peng, A.; Petrone, T.; Henderson, D.; Ranasinghe, V. G.; Zakrzewski, J.; Gao, N.; Rega, G.; Zheng, W.; Liang, M.; Hada, M.; Ehara, K.; Toyota, R.; Fukuda, J.; Hasegawa, M.; Ishida, T.; Nakajima, Y.; Honda, O.; Kitao, H.; Nakai, T.; Vreven, K.; Throssell, J. A.; Montgomery, Jr. J. E.; Peralta, F.; Ogliaro, M.; Bearpark, J. J.; Heyd, E.; Brothers, K. N.; Kudin, V. N.; Staroverov, T.; Keith, R.; Kobayashi, J.; Normand, K.; Raghavachari, A.; Rendell, J. C.; Burant, S. S.; Iyengar, J.; Tomasi, M.; Cossi, J. M.; Millam, M.; Klene, C.; Adamo, R.; Cammi, J. W.; Ochterski, R. L.; Martin, K.; Morokuma, O.; Farkas, J. B.; Foresman, D. J.; Fox, D. J. Gaussian 09, Revision A. 02, Gaussian, Inc.; Wallingford CT 2016.

42. [42]

Suryanarayanan, V.; Singh, S. K. Assessment of dual inhibition property of newly discovered inhibitors against PCAF and GCN5 through in silico screening, molecular dynamics simulation and DFT approach. J. Recept. Signal. Transduct. Res. 2015, 5, 370-380.

43. [43]

Andreas, W. G.; Mark, J. W.; Dong, X.; Duncan, P.; Scott, L. G.; Ross, C. W. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 1. J. Chem. Theory Comput. 2012, 5, 1542-1555.

44. [44]

Salomon, F. R.; Andreas, W. G.; Poole, D.; Grand, S. L.; Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878-3888.

45. [45]

Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015, 8, 3696-3713.

46. [46]

Lindorff, L, K.; Piana, S.; Palmo, K; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins 2010, 8, 1950-1959.

47. [47]

Huang, K.; Luo, S.; Cong, Y.; Zhong, S.; Zhang, J.; Duan, L. An accurate free energy estimator: based on MM/PBSA combined with interaction entropy for protein-ligand binding affinity. Nanoscale 2020, 12, 10737-10750.

48. [48]

Sun, H. Y.; Duan, L. L.; Chen, F.; Liu, H.; Wang, Z.; Pan, P.; Zhu, F.; Zhang, Z. H.; Hou, T. J. Assessing the performance of MM/PBSA and MM/GBSA methods. 7. Entropy effects on the performance of end-point binding free energy calculation approaches. Phys. Chem. Chem. Phys. 2018, 20, 14450-14460.

49. [49]

Fu, L.; Chen, Y.; Guo, H. M.; Xu, L.; Lin, Z. H. A selectivity study of polysubstituted pyridinylimidazoles as dual inhibitors of JNK3 and p38α MAPK based on 3D-QSAR, molecular docking, and molecular dynamics simulation. Struct. Chem. 2021, 32, 819-834.

50. [50]

Lucas, S.; Negri, M.; Heim, R.; Zimmer, C.; Hartmann, R. W. Fine-tuning the selectivity of aldosterone synthase inhibitors: structure-activity and structure-selectivity insights from studies of heteroaryl substituted 1, 2, 5, 6-tetrahydropyrrolo[3, 2, 1-ij]quinolin-4-one derivatives. J. Med. Chem. 2011, 7, 2307-2319.

51. [51]

Yin, L.; Hu, Q.; Hartmann, R. W. 3-pyridyl substituted aliphatic cycles as CYP11B2 inhibitors: aromaticity abolishment of the core significantly increased selectivity over CYP1A2. PLoS One 2012, 11, e48048-e48057.

52. [52]

Lipinski, C. A. Drug-like properties and the causes of poor solubility and poor permeability. J. Pharmacol. Toxicol. Methods 2000, 1, 235-249.

53. [53]

Daina, A.; Michielin, O.; Zoete, V. SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 2017, 7, 42717-42729.

• Figure 1  Structures of spironolactone and eplenone

Figure 2  Structure of compound No. 22

Figure 3  Alignment of CYP11B2 inhibitors

Figure 4  Active pockets produced after protein (4DVQ) pretreatment. The yellow stick represents the original ligand (1CA), and magentas stick represents HEM

Figure 5  pIC50 of predicted versus actual activity of the training and test sets. (A) CoMFA (B) CoMSIA

Figure 6  Three-dimensional equipotential maps of CoMFA(a, b) and CoMSIA(c, d) with No. 22. (a) stereo, (b) electrostatic, (c) hydrophobic and (d) hydrogen bond acceptor

Figure 7  Williams plot for CoMFA(a) and COMSIA(b) models

(The leverage threshold h* = 0.146/0.366 and standardized residual σ = ±2 or 3)

Figure 8  MESP superimposed onto a surface of constant electron density, HOMO and LUMO

Figure 9  Overlay of the original ligand docked (1CA, shown as cyan sticks) and re-docked (shown as wheat sticks) in the aldosterone synthase active site

Figure 10  Interaction between protein and molecule. (a) No. 22, (b) No. 04. The small molecule ligand is shown as green sticks and HEM as magentas sticks. π-π interactions are represented by pink dashed lines and hydrogen bond is the yellow dashed lines

Figure 11  Radar maps of compound No. 22 and its derivative

Figure 12  RMSD and RMSF of complexes docked with compound No. 22/a2 for 100 ns simulations

Figure 13  Interaction between 4DVQ and compound No. 22 before (A) and after (B) 100ns MD simulation

Figure 14  Interaction between 4DVQ and compound No. a2 before (A) and after (B) 100 ns MD simulation

Figure 15  Contributions of free energy calculated by MM/GBSA for key residues of 4DVQ (A) No. 22 and (B) No. a2

计量
• PDF下载量:  0
• 文章访问数:  53
• HTML全文浏览量:  3
文章相关
• 发布日期:  2022-03-01
• 收稿日期:  2021-08-02
• 接受日期:  2021-10-15
通讯作者: 陈斌, bchen63@163.com
• 1.

沈阳化工大学材料科学与工程学院 沈阳 110142

/

• 分享
• 用微信扫码二维码

分享至好友和朋友圈