Designing of the N-ethyl-4-(pyridin-4-yl)benzamide based potent ROCK1 inhibitors using docking, molecular dynamics, and 3D-QSAR

Molecular docking

Molecular docking experiments of the compounds C03, C08, C11, C13, C19, C22, C25, C29, C34, and C37 have been carried out using AutoDockTools (AutoDock 4.2, Scripps Research, CA, USA) as described by earlier studies (Keretsu, Bhujbal & Cho, 2019; Keretsu, Bhujbal & Cho, 2020). The high-resolution crystal structure of ROCK1 (PDB ID:6E9W) was retrieved from the protein databank and used as a receptor in the docking study. Ligands and water molecules were removed from the original structure. The missing loop (residue S249-G253 and residue L372-T279) has been modeled with Modeller-10.1 (University of San Francisco, California, USA) and validated by Ramachandran plot analysis on the PROCHECK (DOE-MBI services, UCLA) server (Fig. S1). The pyridinylbenzamide-based inhibitor, compound C08 in our dataset, is already available in the co-crystallized form. Therefore, C08 was treated as the most selective compound in our dataset. Compound C34 showed the highest inhibitory activity and was designated as the most active compound. Compounds C03, C11, C13, C19, C22, C25, C29, and C37 have different chemical subgroups in their positions R₁, R₂, R₃ and R₄ and exhibited a diverse range of activities against ROCK1, which were also included in the docking study. C08 was extracted and re-docked into the binding pocket of ROCK1. In brief, the protein was prepared by adding polar hydrogen and Kollman charges. The ligand was prepared by the addition of missing polar and nonpolar hydrogen atoms and assigned Gasteiger charges. To avoid excessive conformational explosions, the number of rotatable bonds was fixed to six. The AutoGrid program was used to define the grid box size of 50  × 50  × 60 points in the X, Y, and Z directions with the grid spacing of 0.375, in such a way that it could conceal the entire ligand, while the grid center was set at X = −65.0, Y = −8.0 and Z = −45, respectively. AutoDock−4.2 was used for the conformational search using the Lamarckian Genetic Algorithm (LGA). The other docking parameters were set to default, and 100 docking runs were executed. The most suitable protein-ligand docking conformation was selected from the lowest binding energy with the lowest RMSD cluster of 100 runs. Polar, nonpolar, and hydrophobic interactions were also considered to evaluate the final protein-ligand binding pose. This docking protocol was repeated for compounds C03, C11, C13, C19, C22, C25, C29, and C37. Finally, the docking analysis was performed using AutoDockTools, PyMOL (Schrodinger, Inc.), and LigPlot+ (EMBL-EBI, Cambridgeshire, UK).

MD simulation

The MD simulations were run in GROMACS 2019.5 (Lindahl et al., 2020) using CHARMM36 all-atom Force Field (2020). The ligand topology and parameter files with partial charges were generated using the CHARMM General Force Field (CGenFF) (Massachusetts, USA). The protein-ligand docked complex was solvated with TIP3P water in a cubic periodic box model. The minimum distance between the protein and the box walls was maintained at 10 Å in the X, Y, and Z directions. The Na⁺ and Cl⁻ counterions were added adequately to neutralize the system by replacing the solvent molecules. The entire system was subjected to energy minimization using the steepest descent integrator. Thereafter, the system was equilibrated for 200 ps in the constant volume (NVT) ensemble to achieve the temperature of 300 K, followed by 400 ps in the constant pressure (NPT) ensemble to achieve 1 bar of pressure. The protein backbone was kept restrained during the NVT and NPT simulations. The particle mesh Ewald (PME) algorithm was used to deal with long-range electrostatic interactions and a cut-off value of 12 Å was used to control the van der Waals interactions, respectively. The timestep was set to 2 fs. Finally, the production run has been continued from the NPT ensemble by unrestraining the protein backbone for 50 ns. Temperature and pressure were controlled by the modified Berendsen thermostat (V-rescale) and the Parrinello-Rahman barostat. During the production run, the frames were collected at every 50 ps interval for further data analysis. RMSD, H-bonds, H-bond distance, and H-bond angle of each protein-ligand complex were calculated using the functions ‘gmx rms’, ‘gmx hbond’, ‘gmx distance’ and ‘gmx angle’ built in Gromacs. The simulation protocol was applied to all protein-ligand complexes.

MMPBSA binding free energy calculation

The Molecular Mechanics Poison-Boltzmann Surface Area (MMPBSA) is a popular technique to calculate the binding-free energy between the protein and ligand. We have used the g_mmpbsa (Kumari et al., 2014) package to generate the MMPBSA model of the protein-ligand complex, as described in the earlier study (Keretsu, Ghosh & Cho, 2020). After the total MMPBSA binding free energy was calculated, the residue-wise binding energy decomposition calculation was performed.

The binding-free energy from the MMPBSA method can be expressed as – ${\displaystyle \Delta G_{bind}=\Delta G_{complex}-\Delta G_{protein}-\Delta G_{ligand}=\Delta E_{MM}+\Delta G_{sol}-T\Delta S=\Delta E_{vdW}+\Delta E_{ele}+\Delta G_{GB}+\Delta G_{SA}-T\Delta S}$

In which, the free energy of the protein-ligand complex was expressed by the ΔG_complex. The binding energy of the protein and ligand in the solvent was expressed by ΔG_protein and ΔG_ligand. The ΔE_MM represents the interaction energy between the protein and ligand in the gas phase condition, which was obtained by calculating the van der Waals energy (ΔE_vdW) and electrostatic energy (ΔE_ele). The ΔG_sol denote the free energy solvation which was derived by computing the polar solvation ΔG_GB and non-polar solvation ΔG_SA energy. Finally, TΔ S represents the entropy contribution during the ligand binding to its receptor. Due to the extensive computational process and relatively small impact on the final binding energy calculation, the −TΔS term was ignored in this study.

The interaction energy (IE) was calculated by ‘gmx energy’ function in Gromacs. The IE was estimated by combining the short-range Coulombic interaction (Coul-SR) and short-range Lennard Jones (LJ-SR) energy between the protein and ligand.

IE can be expressed as ${\displaystyle IE=Coul-SR+LJ-SR.}$

Dataset building and molecular alignment

In our current study, we have taken 42 novel pyridinyl-benzamide-based compounds which were reported as ROCK1 inhibitors. The inhibitory activity (IC₅₀) values, which were in µM units, were converted to log IC₅₀ (pIC₅₀) values. The converted pIC₅₀ values covered the activity range of 3 log units (pIC₅₀ = 5.02 to 8.52) for the 3D-QSAR study. The average structure of C08 from the last 1 ns MD trajectory, was taken as a representative compound from the dataset. Based on the structure of the reference molecule, the rest of the compounds were sketched, added hydrogen atoms, and calculated the Gasteiger-Hückel partial atomic charges in SYBYL-X 2.1 (Tripos Inc.) to build the dataset. These compounds were minimized using the tripos force field with a convergence of 0.05 kcal/mol by setting the maximum iteration to 2000 runs. The dataset molecules were further divided into the training set of 32 compounds to build the model and the test set of 9 compounds to evaluate the external predictivity of the 3D-QSAR model. Compound C16 has a nonspecific value and was excluded from the dataset during the 3D-QSAR study.

Molecular alignment of the compounds is a vital step in 3D-QSAR, as it was assumed that, sharing the common scaffold, the compounds would bind to the receptor in a similar orientation (Bang & Cho, 2004). We have used the substructure-based alignment, distill-rigid fit, and database alignment functionality available in SYBYL-X 2.1. After that, the aligned compounds were subjected to the 3D-QSAR study.

3D-QSAR study

CoMFA and CoMSIA are two routinely used 3D-QSAR methods to establish the structure–activity relationships of chemical compounds (Gadhe et al., 2011; Pasha et al., 2007). In both methods, the descriptors are calculated in a three-dimensional cubic lattice with 2 Å grid spacing around the molecules. Steric (S) and electrostatic (E) fields were calculated separately for each molecule in CoMFA, by taking SP³ carbon as a probe with a charge value of +1 and an energy tolerance cutoff of 30 kcal/mol. The Lennard-Jones and Coulomb potential functions were used to calculate the steric and electrostatic fields, respectively. In the CoMSIA model, together with the steric and electrostatic fields, three additional descriptors such as hydrophobic (H), hydrogen bond donor (D), and hydrogen bond acceptor (A) have been calculated. A Gaussian-type function was assigned to all grid points for calculating the similarity indices. The rest of the parameters were kept similar to the CoMFA.

PLS calculation, CoMFA, and CoMSIA model building and validation

The partial least squares (PLS) method was used to derive the statistically significant CoMFA and CoMSIA models, followed by the cross-validation (q²) and non-cross validation (r²) of the training set. In the PLS method, the pIC₅₀ values (dependent variables) of the training set were correlated with their respective CoMFA and CoMSIA descriptor fields (independent variables). Leave-one-out (LOO) was used to obtain the cross-validation coefficient, q², the standard error of prediction (SEP), and the optimal number of components (ONC). Subsequently, a non-cross-validated correlation coefficient, r², was generated for statistical indexing of the predictive power of the CoMFA and CoMSIA models along with the Fisher’s statistics (F-value), standard error of estimation (SEE), and field contribution values of each descriptor. The column filtering value was set to 2.0 kcal/mol for both cross and non-cross validation methods. In CoMSIA, different combinations of descriptors (S, E, H, A, and D) were used to obtain the best statistically significant model (San Juan & Cho, 2007).

Sometimes, higher values of q² and r² may not be sufficient to evaluate the predictive outcomes of CoMFA and CoMSIA. Thus, the external predictive correlation coefficient (r²_pred) was carried out to assess the predictive power of the models. The predictive correlation coefficient is expressed by the following equation: ${\displaystyle r_{pred}^2=\left(\text{SD-PRESS}\right)∕\mathrm{SD}}$ SD stands for the sum of the squared deviation between the pIC₅₀ values of the test set and the mean activity values of the training set. PRESS is the sum of the squared deviations between the observed and calculated activity values of each test set of compounds.

Contour map analysis

The CoMFA and CoMSIA result was represented as a StDev*Coeff contour map to visualize the field effect of the chemical descriptors around the most selective compound C08. The green and yellow contours represented the sterically favorable and unfavorable regions, while the blue and red contours represented the favorable substitution for positively or negatively charged chemical groups, respectively, in CoMFA and CoMSIA (Gadhe, Kothandan & Cho, 2012). Likewise, the favorable and unfavorable positions of the hydrophobic, H-bond acceptor and H-bond donor were represented by yellow-white, magenta-red, and purple-cyan by color scheme, respectively.

New compound design, ADME/Tox, and synthetic accessibility prediction

Based on the CoMFA and CoMSIA contour maps, we have designed forty new compounds and predicted their activity (pIC₅₀) by the CoMFA and CoMSIA models. Compounds with higher pIC₅₀ values were subjected to docking, MD simulation, and MMPBSA analysis.

The newly designed compound has also proceeded with the ADME/Tox prediction and synthetic accessibility analysis by pkCSM (Pires, Blundell & Ascher, 2015) (http://biosig.unimelb.edu.au/pkcsm) and SwissADMET (Daina, Michielin & Zoete, 2017) (http://www.swissadme.ch) online web server.

Molecular docking analysis

The reliability of the molecular docking was verified by the ligRMSD (Velázquez-Libera et al., 2020) server, as the docking pose was further used in the MD study. The docked structures of ten compounds C03, C08, C11, C13, C19, C22, C25, C29, C34, and C37 were compared by aligning them with the C34 crystal pose (PDB ID: 6E9W). The final docked pose of each compound was selected based on the ECIDALs (Muñoz Gutierrez et al., 2016) norms with an RMSD tolerance of 2 Å. The 1-(4-(pyridin-4-yl)phenyl)ethan-1-one scaffold of the docked compounds and crystal ligand was superimposed at the same position. The RMSD of each compound and their binding energy from the docking study were tabulated in Table-S1. Except for compound C22 (RMSD = 2.89 Å), the other compounds have shown the RMSD value of less than 2 Å from the crystal ligand. The redocking of the C08 showed that it was stabilized by several polar and hydrophobic interactions in the active site. It formed the H-bond interaction with M156, located in the hinge loop, and another H-bond interaction with the catalytic lysine K105. The Gatekeeper residue F368 formed the π- π interaction with the pyrimidine ring. Other residues, such as I82, F87, A103, L107, Y155, formed the π-alkyl and hydrophobic interactions. When docking with the most active compound C34, an increased number of H-bonds was observed. Along with M156 and K105, residues R84 and F87 also participated in the H-bond interaction with C34 through their backbone hydroxyl (-OH) and amine (-NH2) groups. The detailed molecular docking interactions of compounds C08 and C34 have been illustrated in Fig. 2. The 2D molecular docking interactions of the rest of the compounds C03, C11, C13, C19, C22, C25, C29, C37 were shown in Fig. S2. Residues that participated in hydrophobic and H-bond interactions were shown in green and magenta colors, respectively, in Figs. 2A and 2D. The H-bond interactions were shown in red dotted lines in Figs. 2B and 2E. The active site cavity was shown by the hydrophobic surface model based on the Eisenberg scale in Figs. 2C and 2F, where white to red indicated the increase in hydrophobicity. All docked protein-ligand complexes were taken for the MD simulation study.

MD simulation analysis

Since protein-ligand binding is a dynamic process, and a single docking result would have remained inconclusive. Therefore, we performed MD simulations for all ten protein-ligand complexes for 50 ns. All systems were equilibrated within the initial 10 ns of simulation and thereafter maintained a stable plateau until the end of the run. The RMSD plot of each protein-ligand system was shown in Fig. S3. The fluctuation of RMSDs for the protein was observed within the range of 1.0–4.5 Å, whereas the fluctuations of the ligands were found to be in the range of 0–3.0 Å, respectively. To study the in-depth molecular interaction, we compared the binding pose from the MD trajectory (average structure taken from the final 1 ns) and the docking pose of each compound by superimposing them. The morphed complexes were shown in Fig. 3.

The conformations of compounds C03, C11, C19, C34, and C37 (Figs. 3A, 3C, 3E, 3I, 3J) from the MD simulation displayed close similarity with the corresponding binding pose from the docking study. The H-bond interactions with K105 and M156 remained intact until the final step of the trajectory. In contrast, compounds C08, C13, C25, and C29 Figs. 3B, 3D, 3G, 3H) lost the H-bond interaction between the carbonyl group (-C =0) and K105. However, the H-bond interaction with the backbone nitrogen of M156 remained intact. In addition, compound C25 formed the H-bond interaction with the backbone nitrogen atom (N-atom) of R84 by the carbonyl group. For C22, which is the least active compound in the dataset, the two H-bond interactions with K105 and M156 were found to be broken in the MD simulation analysis.

The overlapped structure (Fig. 3F) showed that the distance between the hinge loop and the pyrimidine-2-amine ring of C22 increased and the methoxybenzene moiety was flipped outward from its docked position. However, C22 remained in the active site throughout the production run, mainly by forming polar and hydrophobic interactions with the surrounding residues. Because M156 is the key residue in forming the H-bond interaction, we examined the distance between the N atom of the M156 backbone and the N atom of the pyrimidine ring of the compounds throughout the simulation in Fig. 4. Aside from C22, we found that the H-bond interaction distances from the residue M156 varied between 2–3.3 Å for other compounds. In the case of compound C22, the H-bond interaction distance was increased to ∼4.9 Å. However, we did not see any major changes in the RMSD plot of C22 as observed in the Fig. S3F. One of the key precedents to determine the existence of the H-bond formation is the calculation of the angle between the donor hydrogen and acceptor (D-H-A) atoms, which was expected to be linear or close to 180°. Figure S4 depicted the distribution of the angle formation among the D-H-A atoms for each compound. Furthermore, we have calculated the relative frequency of the total number of H-bond interactions between ROCK1 and selected compounds in Fig. S5. Throughout the production simulation, the frequency of H-bond interactions between ROCK1 and selected compounds was found to be in the range of 0 to 4. C19 and C22 showed a higher frequency of completely broken H-bond interactions (10% and 16%, respectively) compared to other compounds.

MMPBSA binding energy and IE estimation

We have calculated the MMPBSA binding-free energy and IE from the snapshots extracted from the final 2 ns of the MD trajectory of each system, to estimate the protein-ligand binding affinity. Electrostatic energy, van der Waals energy, and SASA energy terms were the main contributing factors to total BE. Alternatively, the short-range Coulombic (Coul-SR) interaction and the short-range Lennard-Jones (LJ-SR) energy contributed to the total interaction energy. The details of MMPBSA binding free energy and IE were presented in Table 1. The total MMPBSA binding free energies were found to be −107.85 kJ/mol, −109.22 kJ/mol, −114.08 kJ/mol, −102.36 kJ/mol, −76.09 kJ/mol, −95.44 kJ/mol, −123.18 kJ/mol, −108.85 kJ/mol, −130.61 kJ/mol, −106.69 kJ/mol for compounds C03, C08, C11, C13, C19, C22, C25, C29, C34 and C37, respectively. The IE was found to be −253.41 kJ/mol, −247.97 kJ/mol, −279.69 kJ/mol, −194.26 kJ/mol, −215.37 kJ/mol, −210.93 kJ/mol, −284.75 kJ/mol, −226.21 kJ/mol, −284.21 kJ/mol, −233.93 kJ/mol for the compounds C03, C08, C11, C13, C19, C22, C25, C29, C34 and C37, respectively. The most active compound, C34 resulted in the lowest binding-free energy in the MMPBSA model. In addition, with C25, C34 provided the lowest IE values among the other compounds. However, we did not find any correlation between the pIC₅₀ and MMPBSA or IE values. The common residues, which contributed to the noticeable amount of net positive and net negative BE energy were shown in Table S2. The graphical comparison of BE decomposition by these residues was shown in Fig. 5. Residues E89, V90, E124, D126, Y155, D160, D202 contributed the lowest BE, while residues K105 and R125 contributed the highest BE to the ligands. However, we observed a large variation in BE decomposition by residue D216. In the formation of the complex with C03, C11, C25, and C34, it provided net positive BE values in the total binding energy. Conversely, it contributed negative BE in compounds C08, C13, C19, C22, C29, and C37, respectively.

Table 2:

Structure and activity values of N-ethyl-4-(pyridin-4-yl)benzamide based ROCK1 inhibitor.

DOI: 10.7717/peerj.11951/table-2

Molecular alignment and 3D-QSAR studies

The coherence selection of the compounds and the splitting of them into a training set and a test set is an important step toward the development of the 3D-QSAR. The dataset compounds and their corresponding pIC₅₀ values were shown in Table 2. We indexed the compounds into high, medium, and low activity ranges and chose the test set compounds so that they can cover different activity values while maintaining structural variations. All training set compounds were aligned to the representative compound C08 by their common substructure N-isopropylbenzamide. The aligned compounds and their common substructure has been illustrated in Fig. S6.

We have produced the CoMFA and CoMSIA models based on the 32 compounds in the training set. For model validation, we have restricted the statistical parameters of q² >0.6; r² >0.8 and r²_pred >0.5. Initially, we applied different charge fields, namely, Gasteiger, Gasteiger-Hückle, Pullman, MMFF94 charges, to obtain the best possible CoMFA model. In the cross-validated model, we obtained the final q² = 0.774, an ONC value of 6, and a SEP value of 0.425 in the Gasteiger charge field. In the non-cross-validated model, the r² value and the SEE value were found to be 0.965 and 0.167, respectively. We have employed the same charge to build a reliable CoMSIA model with multiple combinations of steric (S), electrostatic (E), hydrophobic (H), and H-bond acceptor (A) and H-bond donor (D) fields.

Only the steric and electrostatic field produced a reasonable CoMSIA model with q² = 0.676 and ONC = 6 from cross-validation and r² = 0.949 from non-cross-validation, respectively. The r²_pred of CoMFA and CoMSIA (SE) was found to be 0.703 and 0.548, which is higher than the set criteria of 0.5. Overall, both CoMFA and CoMSIA exhibited the robustness of the external predictivity. The detailed statistics of the 3D-QSAR statistics were shown in Table 3. The statistics of the different combinations of descriptors in the CoMSIA model have been tabulated in Table S3. The actual pIC₅₀ and predicted pIC₅₀ from CoMFA and CoMSIA with their residuals were shown in the Table S4.

CoMFA and CoMSIA contour map and SAR analysis

From the 3D-QSAR study, we have plotted the correlation plots between the actual pIC50 and the predicted pIC₅₀ in Fig. 6. The compounds in the training and test sets have been labeled dark and light green, respectively. The contour maps generated by CoMFA and CoMSIA around the compound C08 are shown in Figs. 6B, 6C, 6E and 6F. From the CoMFA model, we obtained the steric and electrostatic field contributions of 63% and 37%, respectively, and from the CoMSIA model, we found the steric and electrostatic field contribution of 57.2% and 42.8%, respectively. A green contour map indicates a favorable substitution for the bulky steric groups, whereas a yellow contour indicates that the steric substitution would be unfavorable. Similarly, a blue and red contours represent favorable and unfavorable substitutions for positively charged electron-donating groups, respectively, to increase the inhibitory activity of the molecules. The green contour was observed at position R₄ (Figs. 6B and 6E) near F120, which indicated that the steric groups could be favorable at that position. Compounds C25 and C26 with a propylpyrrolidine moiety at the R₄ position exhibited higher activity values than Compounds C01-C04 and C08-C11 with smaller groups present at their R₄ position. A green contour at the R₃ position near residue V90 also indicated a favorable substitution for the steric groups. The yellow contours were observed near the M156 position in R₂ and near the R84 residue in R₁ position in the CoMFA model, which represented the unfavorable position for steric substitution. The contour maps suggested that a steric group would be unfavorable for R₂ substitution. This could be the reason why compounds C34, C36, and C25-27 had higher activities than C20-24. Figures 6C and 6F show a blue contour covering the dioxane ring at the R₁ position, and two red contours near K105 and R84, indicating that the electron-donating and electron-accepting group replacement may increase the effectiveness of the inhibitors.

Designing of new compounds

Based on the CoMFA and CoMSIA contour map analysis, we have acquired an ideal SAR scheme design a more potent ROCK1 inhibitor drawn in Fig. 7, based on compound 08. A steric and an electropositive group would be preferred for the R₂ and R₄ positions. A steric-electropositive group at the R₃ position could increase the activity of existing compounds. As the CoMFA scheme showed higher statistical significance (q² and q²) and predictive power (q²_pred), the designed compounds were evaluated by the CoMFA model. Based on the scheme in Fig. 7, we have designed forty new compounds and predicted their activity values in the CoMFA model, Table S5. Of these, seven designed compounds, D02, D03, D06, D31, D32, D33, D35, have exhibited higher pIC₅₀ values than the most active compound C34, as shown in Table 4. We have utilized the pkCSM and SwissADME server to analyze the ADME/Tox and SA scores. For comparison, we have also taken the compounds C08, C34, and C37 from the dataset. The detailed analysis has been tabulated in Table S6. Except for compound D06, none of the designed compounds showed AMES toxicity. The SA score predicts the synthetic accessibility of small molecules. A SA score of 1 denotes the ease of synthesis and an SA score of 10 expresses the difficulty of synthesizing the chemical compound. The SA score of the designed compound was predicted to be less than 5, which stated that the compounds would have a low to moderate difficulty in the synthesis routes. We have also examined the docking (Table S7) and MD studies on these seven newly designed compounds to estimate their binding interaction. The RMSD graph was shown in Fig. S7. The MMPBSA and the interaction energy were shown in Table 5. The residue-wise binding energy decomposition was shown in Table S9. The catalytic lysine K105 has exhibited a higher positive binding energy decomposition among the selected residues. Finally, we have taken the last 1 ns average structure of the protein and designed ligand complexes to study the key interactions at the residue level. The designed compounds tend to form a higher number of H-bond interactions with the surrounding residues compared to the dataset compounds. The in-depth interactions in the active site pocket were illustrated in Fig. 8.

Table 4:

Newly designed compounds with higher predicted pIC₅₀ values.

DOI: 10.7717/peerj.11951/table-4