Design of new captopril mimics as promising ACE inhibitors: ADME, pharmacophore, molecular docking and dynamics simulation with MM-PBSA and PCA calculations

Abstract

New pyrrolidine derivatives with more than 50% structural similarity with captopril were designed to get new captopril mimics with superior potential to act on both peripheral and central ACE. Further optimization was carried out through pharmacophoric mapping, then pharmacokinetics of these compounds were analyzed, 42 derivatives were selected for further study, as they exhibited potential to pass through BBB. Molecular docking on ACE using captopril and lisinopril as reference drugs was performed, and Compound 28 (2-Pyrrolidin-2-ylidene-N-thiomorpholin-4-ylmethyl-malonamic acid ethyl ester) showed the best docking scores, proving its superiority over captopril and comparability to lisinopril. Further molecular dynamics simulations and energy calculations demonstrated binding stability and close mimicry to both drugs. The results indicate that Compound 28 is a promising candidate for further investigations as a potential drug to act centrally and peripherally. Compound 28 can be synthesized by reacting Cyano-pyrrolidin-2-ylidene-acetic acid ethyl ester through Mannich reaction with thiomorpholine and formaldehyde.

KEYWORDS:

• Captopril mimics
• Design of ACE inhibitors
• Molecular docking
• Virtual screening
• Angiotensin Converting Enzyme (ACE)
• Pyrrolidin derivatives

1. Introduction

The burden of cardiovascular disorders (ischemic heart disease and stroke) on health is evident; WHO (world health organization) has considered these disorders as common causes of death worldwide, with an estimated 17.8 million deaths annually [1]. The mortality rates over the past decade leading to 2017 present an alarming 21% increase. The burden of cardiovascular diseases is not exclusively derived from the fatalities they cause; morbidity and disability associated with cardiovascular diseases are major aspects of the illness's epidemiology [2].

Generally, management and treatment of such cases begin with clinical manifestations. However, the problem in this manifestation is often sudden death [3]. As a result, the focus is shifted towards the underlying risk factors such as hypertension which affects one-third of the adult world population [4]. As in gestational hypertension, several conditions cause secondary hypertension, which is reversible and resolves eventually after the cause passes [5,6]. Most cases suffer from primary or essential hypertension, which is irreversible and multifactorial [7]. The management of this menace includes; (1) non-pharmacological lifestyle intervention such as weight loss, stress reduction, and physical activity [4]. (2) pharmacological drug therapy, including Ca channel blockers, diuretics, and renin-angiotensin system inhibitors (angiotensin-II receptor blockers and angiotensin-I converting enzyme blockers) [8].

Captopril is the first orally active peptide inhibitor of the angiotensin-I converting enzyme (ACE) used to manage hypertension [9]. It is characterized by rapid oral absorption and an exaggerated fall in blood pressure following the initial dose [10]. Additionally, research has shown an intriguing property resulting from the captopril-induced ACE inhibition in the brain [11]. These data indicated that captopril and other ACE inhibitors possess neuron effects, resulting in the retardation of signs of neurodegenerative disorders such as Alzheimer's disease [12,13]. These features made captopril an intriguing prospect for derivatization to attain more potent hypotensive agents with prospective neuroprotective function, significantly benefiting Alzheimer's disease and other disorders. However, one of the limiting factors of its development is the side effects attributed to its mercapto group, which caused allergic reactions, taste disturbance and indirect weight loss, and the extensive metabolism of the inactive disulphide derivative [14–16].

Our design maintained the pyrrole moiety of captopril as the cornerstone. The acidic COOH was substituted with other moieties, like the nitrile group. Substituents were distributed to either the pyrrole nitrogen atom or in the place of the nitrile group adjacent to the pyrrole, as shown in Figure 1. Furthermore, the compounds were designed without the mercapto group to minimize any possible side effects owing to the group. These compounds were then analyzed through various in silico studies.

Figure 1. Designed captopril mimics.

2. Methods

2.1. Design of captopril mimics

The compounds were designed using captopril substructure, and a total of one hundred and two analogues were obtained. These newly designed mimics will be subjected to several CADD techniques, such as pharmacophore mapping, virtual screening, molecular docking and molecular dynamics to get the best hits to be synthesized.

2.2. E-pharmacophore generation and screening

The protein containing angiotensin-1 convertingenzyme (ACE) in complex with captopril was obtained from the protein data bank (PDB ID: 1UZF) containing captopril and ACE at a resolution of 2 Å [17]. The protein was preprocessed, solvent molecules beyond 5 Å of heteroatoms were removed, and missing hydrogens and loops were rectified. After minimization using OPLS3e, the structure was imported into “Develop pharmacophore from protein–ligand complex” for the e-pharmacophore generation in Schrodinger Maestro [18], and the pharmacophore was mapped. The validation step was then performed using two data sets; the first set contained eight known ACE inhibitors (captopril, enalapril, lisinopril, fosinopril, ramipril, quinapril, benazepril, and trandolapril). While the second had four hundred structurally similar inactive dummies generated by the DUD-E website (www.dude.docking.org).

2.3. ADMET descriptors

Multiple descriptors of ADMET properties (absorption, distribution, metabolism, excretion and toxicity) of the compounds were calculated using the “Schrodinger Maestro” suite. “QikProp” protocol utilizes the 3D structure of the compounds to predict relevant properties such as log S (aqueous solubility), log Khsa (human serum binding albumin), SASA (total solvent-accessible area), log HERG (predicted IC₅₀ for HERG K⁺ channel blockage), pKa (Supp. Data, Table 1S) and many other pharmacokinetic features. It provides an intrinsic scoring function (#stars) that compares each calculated descriptor with its optimal range for 95% of known drugs [19–22].

2.4. Molecular docking

The file containing the human angiotensin-converting enzyme in complex with lisinopril was downloaded from the protein data bank (PDB ID: 1O86) [23]. The screened compounds were docked via MOE 2019 software [24]. The software's default approach for structural optimization was used to create and optimize both receptors and compounds. After creating a library of chemicals, the software's usual structure optimization methodology was used to generate both receptors and ligands. Using the AMBER12: EHT force field, energy minimization was carried out, and the active site was assigned to the pocket enclosing the co-crystallized ligands. The “Induced-fit protocol” of MOE DOCKTITE Wizard's induced fit protocol, which uses the GBVI/WSA dG scoring function, was used to carry out the docking. Re-docking lisinopril into the binding pocket and calculating the root mean square deviation (RMSD) between poses were both used to accomplish validation [25–29].

2.5. Molecular dynamics simulations and molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) calculations

Molecular dynamics simulations were performed using GROningen MAchine for Chemical Simulations (GROMACS 2022.3) software [30,31]. Enzyme chains containing active sites were isolated. The enzyme and the ligand topologies were generated and joined using the standard published protocol under the Amber 99SB-ILDN force field [32]. Ligands topology files were successfully generated using the ACPYPE-AnteChamber PYthon Parser interface using General Amber Force Field (GAFF) [33]. The typical workflow of GROMACS enzyme-ligand simulations was applied in which systems were solvated using a TIP3P (three-site transferrable intermolecular potential) water model in a dodecahedron box Na⁺/Cl⁻ ions were sufficiently added for neutralization [34]. The energy minimization (1000.0 kJ/mol/nm) was performed using the steepest descent method, followed by a canonical ensemble (NVT at 300 K) and isothermal-isobaric (NPT at 1 bar) to pre-balance the system successively. All systems were simulated with 50 ns MD simulation. After the MD simulation, the GROMACS package built-in tools were deployed to evaluate the trajectories in terms of RMSD, the radius of gyration (Rg), and root mean square fluctuation (RMSF) after residues numbers were corrected [20,35]. While “gmx covar” and “gmx anaeig” functions were used for the implementation of Principle component analysis (PCA). Gibbs free energy landscapes (FEL) were calculated through “g_sham” using the predicted distribution probability of the top two principal components (PCs) [36].

The MM-PBSA approach calculates the interaction free energy of protein–ligand complexes (potential, polar-solvation and non-polar solvation energies) for molecular dynamics simulation trajectories. The gmx_MMPBSA (1.5.7) tool was employed for the last 2ns of the simulation processes (48–50 ns) with two frame intervals [37]. The trajectories were prepared by removing PBC conditions, and calculations performed using Amberff19SB and GAFF forcefields were used for proteins and ligands under the dielectric model (ipb = 2) [38]. Binding free energies were then calculated using either of the following two equations: $\begin{array}{rl}\mathrm{\Delta }G(\mathrm{binding})& =G(\mathrm{complex})-G(\mathrm{receptor})-G(\mathrm{ligand})\end{array}$ $\begin{array}{rl}\mathrm{\Delta }G(\mathrm{binding})& =\mathrm{\Delta H}-T\mathrm{\Delta S}\end{array}$ (ΔH corresponds to the enthalpy of binding and −TΔS to the conformational entropy after ligand binding)

In which to calculate each term of the above equation is calculated individually through the application of several procedures to get the energy terms for complex, receptor and ligand as described in the literature [39].

2.6. Density functional theory (DFT) analysis

The Spartan ‘14 programme was used to perform quantum chemistry calculations using the DFT method. Spartan ‘14 was used to display all the data files. The density functional theory (DFT) at 6-311G⁺⁺(d,p) basis set/B3LYP approach was utilized to optimize the organic chemical structure of the compounds under investigation, and Chem3D 15.0 software was used to create the original chemical structure [40].

3. Results and discussion

3.1. Design of captopril mimics

These new mimics (1-102) shown in Figure 2 were designed to resemble the captopril 2D structure, keeping the pyrrolidine moiety as the main scaffold. Several computational techniques will explore these new mimics to find the most hopeful candidates.

Figure 2. Designed captopril mimics.

3.2. E-pharmacophore generation and screening

The generated pharmacophore model showed three features in the form of a triangle with features around 3.00–3.48 Å apart from each other. When superposed on captopril, features were distributed as follows, hydrophobic feature in the pyrrolidine residue, negative ionic feature in the carboxylic region, and acceptor feature in the amide carbonyl group. The area under the accumulated curve (AUAC) and receiver operating characteristic (ROC) were calculated for the hypothesis and yielded 0.74 for both. The library was screened, and all designed compounds fitted the designed pharmacophore model; Figure 3 represents the fitted pharmacophoric features for the target compounds and illustrates that the designed mimics share more than 70% of captopril features.

Figure 3. Shared pharmacophoric features of the designed compounds with captopril.

3.3. Drug-likeness study

Lipinski's rule of five measures the “Drug likeness” of compounds through four properties (molecular weight, lipophilicity, hydrogen bond donors, and hydrogen bond acceptors) [41]. All the designed compounds entirely adhered to the rule of five (Table 1S), and no violations were observed. The # stars count provides a comprehensive overview of the inclination of the compounds toward attaining favourable pharmacokinetic and metabolic profiles. The highest portion of our compounds, about 92.16% (ninety-four compounds), achieved high similarity of the calculated properties to the standard ranges of known drugs. At the same time, 5.88% (six compounds) had one property outside the acceptable 95% range of known drugs. Only 1.96% (two compounds) had a higher violation score of 2. Additionally, the “Rule of three”, determining the probability of compounds’ oral bioavailability when complying with certain predetermined factors, was enforced. Only 77.45% (seventy-nine compounds) conformed to the rule of three, while the other 22.55% (twenty-three compounds) had one descriptor, only violating the rule (Figure 4).

Figure 4. Drug-likeness studies for the designed library.

The blood–brain barrier permeability of the library revealed the ability of almost half of the designed compounds to reach CNS 45.10% (forty-six compounds). Most of the compounds showed good metabolic steps (<8), with three metabolic steps being the highest frequent metabolic order. The compounds showed high safety profiles with no predicted blockage of HERG K+, which can cause QT prolongation, ventricular arrhythmias, and sudden death [42]. Moreover, 98% (one hundred compounds) showed no plasma protein binding, adding to the safety profile of the design as well as access to biological targets.

All the above rules guided us to select the most promising drug-like compounds for further analysis (Figure 5). This resulted in forty-two compounds being promoted for molecular modelling study.

Figure 5. The selected most promising drug-like compounds.

3.4. Molecular docking

The study was performed on a human angiotensin-converting enzyme (ACE) obtained from the protein databank (PDB ID: 1O86). It contained the enzyme complexed with Zinc ion (Zn701) and lisinopril. Lisinopril, a known inhibitor of this metalloenzyme, was used as a reference for docking validation and for comparative analysis of the activity of our compounds in addition to captopril. After the compounds and the receptor were prepared appropriately for docking through MOE. The first step analyzed the validity of our docking protocol by calculating RMSD between the co-crystallized and docked pose of lisinopril, which resulted in 0.73 Å, indicating successful validation (Figure 6). Afterwards, the compounds’ docking commenced, and the results are presented in Table 1. Although some research highlighted the beneficial outcomes of setting restraints during molecular docking [43], our study was performed without such restraints encouraged by other research studies with similar practical aspects [44] and precision of validation. All compounds achieved favourable binding to ACE with negative scores (Kcal/mole). The two drugs used for comparison, namely captopril and lisinopril, scored −6.54 and −12.47 Kcal/mole, respectively. Nearly all compounds (34 compounds, 80%) were able to produce superior scores than captopril, while the remaining eight compounds produced close submarginal binding scores. The top-scoring compounds, 28 and 18, achieved close results to lisinopril, with one even out-performing scoring −12.69 and −11.81 Kcal/mole, respectively.

Figure 6. Validation of docking protocol showing RMSD of 0.73 Å. (Green is the co-crystallized pose, and pink is the re-docked pose).

Table 1. Docking scores of compounds against human ACE in Kcal/mole.

Compound	Score	Compound	Score	Compound	Score	Compound	Score
1	−7.29	12	−6.99	25	−7.68	88	−7.64
2	−6.77	13	−9.06	26	−7.72	90	−6.91
3	−10.33	14	−8.10	27	−7.40	93	−6.66
4	−6.83	17	−7.54	28	−12.69	95	−7.55
5	−7.75	18	−11.81	29	−8.91	96	−7.45
6	−7.34	19	−8.31	30	−8.32	97	−6.92
7	−6.60	20	−7.66	73	−6.80	98	−8.09
8	−7.77	21	−7.50	76	−7.25	100	−7.54
9	−8.85	22	−7.55	77	−6.48	102	−7.46
10	−8.46	23	−8.14	80	−7.03	Captopril	−6.54
11	−8.98	24	−7.13	85	−7.26	Lisinopril	−12.47

Analysis of their binding interactions (Table 2) illustrated the binding of both drugs to Zn701 and the formation of several hydrogen bonds to HIS353 and TYR520. Additionally, the extended nature of lisinopril offered extra binding ionically with GLU162, ASP377, and LYS511. As well as several hydrophobic interactions with HIS383 and VAL518.

Table 2. Binding interactions of lisinopril, captopril and compounds 28 and 18.

	Interactions
Lisinopril
Captopril
28
18

The highest-scoring compound, 28, interacted with Zn701 through the nitrogen atom of its thiomorpholine ring. Furthermore, its double carboxylic groups acted as hydrogen bond acceptors to imidazole rings of HIS353 and HIS513 and the phenolic hydroxyl group of TYR520. On the other hand, Although 18 possessed similar moieties, the placement of the acetamido group didn't allow for the formation of hydrogen bonds with TYR520, thus decreasing its score, as evident in Figure 7.

Figure 7. Superposition of 18 (pink) and 28 (cyan) binding to ACE.

The interesting properties of 28 encouraged further exploration and analysis of its binding to ACE using subsequent molecular dynamic simulation and energy calculations.

3.5. Molecular dynamic simulations and molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) calculations

The computerized approach known as molecular dynamics (MD) is used to investigate the dynamic behaviour of complex systems where atoms and molecules interact over time, thus, generally increasing the effectiveness of drug discovery. In contrast to standard molecular docking techniques, MD simulations consider the targets’ flexibility. When coupled with estimations of binding energies, a more precise prediction of potential inhibitors can be made [45]. ACE (PDB ID: 1O86) in complexes with captopril, lisinopril and 28 were simulated under normal physiological conditions for 50 ns each. Afterwards, stability, dynamic behaviour, and compactness of protein-drug complexes were assessed using structural characteristics such as RMSD, RMSF, and Rg.

3.5.1. Root mean square deviation (RMSD)

To determine whether or not ACE is stable in complex with the chosen compounds, we measured their RMSD of the protein backbone. The average RMSD of captopril, lisinopril and 28 complexes produced close values of 0.14±0.01, 0.16±0.01 and 0.16±0.03 nm, respectively. As shown in Figure 8, the similarity in the behaviour of ACE across the three simulations illustrates the similarity of 28 in terms of effects on the protein to both reference drugs. Furthermore, each compound's RMSD was calculated throughout the simulation in which values of 0.10±0.02, 0.10±0.02 and 0.17±0.03 nm were observed for captopril, lisinopril and 28, respectively (Figure 9). These homogenous data provide evidence of the stability of the three ACE complexes.

Figure 8. RMSD chart of ACE in complexes with captopril, lisinopril and 28 for 50 ns simulation.

Figure 9. RMSD chart of captopril, lisinopril and 28 for 50 ns simulation.

3.5.2. Root mean square fluctuations (RMSF)

The RMSF analysis was performed, and the values were plotted against the residue number to compute the residual mobility and behavioural similarity of 28 to the two reference drugs. Overall, the RMSF shows consistent trends in all three simulations, with residues involved in interactions mentioned before showing low fluctuations, such as HIS513 and TYR520 (Figure 10).

Figure 10. RMSF chart of ACE residues in complex with captopril, lisinopril and 28.

3.5.3. Radius of gyration (Rg) and solvent accessible surface area (SASA)

The radius of gyration (Rg) and solvent-accessible surface area (SASA) are useful tools for gaining insight into the folding and compactness of protein–ligand complexes. The ability of a medication to cause conformational changes in proteins can be evaluated using Rg and SASA. The Rg value indicates how loosely packed the molecules are within a protein, while SASA predicts the extent of conformational changes that occur during binding. When simulating a stable ligand–protein interaction, the Rg and SASA values plateau [36]. Figures 11 and 12 show that the simulated structures of all complexes are relatively conservative and stable at sizes of around 2.41 nm and ∼245 nm² for Rg and SASA, respectively.

Figure 11. The radius of gyration of ACE in complex with captopril, lisinopril and 28.

Figure 12. The solvent-accessible surface area of ACE in complex with captopril, lisinopril and 28.

3.5.4. Principal component analysis (PCA)

Throughout the simulation, PCA was used to evaluate the linked motions. Both eigenvalues and eigenvectors were determined as projections of the proteins’ Cα atoms. Eigenvectors define motion orientation, whereas eigenvalues represent motion magnitude. The projected motions of ACE when complexed with 28 express high similarity to when the metalloenzyme is complexed with either reference drug, as shown by the overlapping in eigenvectors in Figure 13. These overlaps suggest that when attached to 28 and reference inhibitors, the proteins act similarly, creating persistent clusters and occupying the same subspace, indicating the stability of the created complexes. Individual examination of the correlated motion of the first ten eigenvectors (which account for 80% of the variations) reveals that 28 displayed moderate relative motion with residues lower than captopril and slightly greater than lisinopril (Figure 13).

Figure 13. Eigenvector and eigenvalues of ACE in with captopril, lisinopril and 28.

The proteins’ conformational space was analyzed in terms of energy and time using principle component (PC) analysis. PC1 and PC2 were then utilized to determine Gibbs's free energy in each case, and the Free Energy 3D Landscape (FEL) was plotted. These energy graphs depict high and low energies in red and blue colours. Evaluation spatial positioning of the of atoms in the system using Gibbs FEL demonstrates that complexes formed metastable formations that cross multiple energy barriers indicating high stability. Compound 28, captopril and lisinopril all have comparable energy surfaces with canonical-shaped energy minima indicating the presence of stable minima (Figure 14).

Figure 14. Gibbs free energy landscapes for ACE when complexed with captopril, lisinopril and 28.

Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) is one of the most often used methods for calculating binding free energy. The lower the predicted binding free energy of a ligand–protein complex, the more stable that complex is expected to be and the higher the ligand's activity and potency. The energy values of captopril and lisinopril produced consistent values of −5.36±2.75 and −5.04±1.64 kcal/mole, respectively. While the complex with compound 28 resulted in a slightly lower value of −6.39±1.50 Kcal/mole. Consistent with previous results, the values obtained exhibit relatively similar behaviour in both cases (Table 3).

Table 3. MMPBSA calculations of ACE in complex with captopril, lisinopril and 28.

	Captopril	Lisinopril	28
ΔEVDWAALS	−21.60±1.93	−38.41±0.55	−34.11±1.30
ΔEEL	−10.27±1.73	−51.16±1.30	−2.83±0.42
ΔEPB	11.74±0.91	58.18±0.78	9.72±0.05
ΔENPOLAR	−17.88±0.18	−35.77±0.23	−23.54±0.34
ΔEDISPER	32.66±0.08	62.12±0.16	44.37±0.52
ΔGGAS	−31.87±2.59	−89.57±1.41	−36.94±1.36
ΔGSOLV	26.52±0.93	84.53±0.83	30.55±0.62
ΔTotal	−5.36±2.75	−5.04±1.64	−6.39±1.50

*E_VDWAALS (Vander Waals), E_EL (Electrostatic energy), E_PB (Polar solvation energy), E_NPOLAR (Non-polar solvation energy), G_SOLV (Total solvation free energy), G_GAS (Total gas phase free energy), and E_DISPER (Energy dispersion).

3.6. DFT analysis

The DFT/B3LYP approach was used in the current work to perform quantum chemical computations to optimize the selected structures. The DFT (B3LYP) method with 6-311G++(d,p) basis set was applied in this test. The optimized structure and its HOMO and LUMO values are represented in Figure 15 and Table 4. The HOMO energy expresses the ability of the compound to give electrons as an electron donor. It is localized mainly on the thiomorpholine moiety of compound 5 and the mercaptomethylpropanoyl residue of captopril.

Figure 15. DFT calculation for compounds 28 and captopril.

Table 4. The molecular properties of compounds 28 and captopril after DFT calculations.

	Energy	Energy solvation	E Humo	E Lumo	Dipole moment	No. of
	(au)	(KJ/Mol)	(ev)	(ev)	(debye)	conformers
28	−1335.77952	−97.36	−5.39	−1.08	6.96	6912
Captopril	−1030.21773	−45.09	−6.18	0.09	4.11	2592

On the other hand, the LUMO energy displayed by a site that has the ability to act as an electron attractive, i.e. electron acceptors due to vacant orbitals and localized at L-proline residues of compound 28 and captopril, respectively. The electrostatic potential maps of the compounds showed areas with electron localization throughout the molecules, with red and blue colours representing electron-rich (negative) and deficient (positive) localization, respectively. Finally, the DFT calculations revealed favourable energetic parameters for the selected compound (5), Table 4.

4. Conclusion

Our journey started with a search for suitable ACE inhibitors that can act peripherally and centrally; a library of one hundred and two captopril analogues was designed. The library was then subjected to several investigations through pharmacophoric screening, drug-likeness, molecular docking, and molecular dynamic simulation. Most of the designed compounds conformed to the generated pharmacophore. The subsequent drug-likeness study resulted in the most promising forty-two compounds afterwards. Molecular docking illustrated the compounds that can occupy the same space as captopril and showed their ability to interact with the metalloenzyme to various extents. The most promising analogue was compound 28, which out-performed captopril and produced comparable results to lisinopril while attaining higher CNS permeability, thus, a higher potential to act centrally. The molecular dynamic simulation was performed for 50 ns to assess the stability of the complexes. All calculations, such as RMSD, RMSF, SASA and Rg, supported the stability of binding and the close mimicry of 28 to both reference drugs. Additionally, MM-PBSA calculations, PCA and FEL also supported our postulations. Finally, DFT expanded on its favourable energetic parameters as well. Thus, the potential of thiomorpholine pyrrolidine derivative (compound 28) to target both central and peripheral ACE encourages further testing and analysis for approving this hopeful captopril mimic as a future antihypertensive drug.

Acknowledgments

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number: IFP22UQU4290565DSR235.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All data are incorporated in the manuscript and supplementary data files.

Additional information

Funding

The Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia has supported this research work through the project number: IFP22UQU4290565DSR235.