Biomedicines | Free Full-Text | Targeting Allosteric Site of PCSK9 Enzyme for the Identification of Small Molecule Inhibitors: An In Silico Drug Repurposing Study

SBVS is a computational approach utilized to predict optimal ligand–target interactions and the formation of complexes [44]. SBVS enables the ranking of ligands based on their affinity to the target, with the most promising compounds appearing at the top of the list. This strategy involves selecting compounds from a database and categorizing them according to their affinity to the receptor site. In our study, we employed a structure-based in silico screening approach on a library of FDA-approved drugs consisting of 2992 compounds obtained from the Selleckchem chemical database, aiming to identify novel therapeutic candidates with good binding affinity against the PCSK9 enzyme. The screening process involved the docking of all FDA-approved drugs to the PCSK9 enzyme using the Glide extra-precision (XP) docking mode. From the initial screening, 268 compounds were selected based on a docking score better than −7.5 kcal/mol for further docking using Smina docking software (Version 1.1.2). Subsequently, the molecules were ranked based on their average docking scores obtained from both software and the poses of the first 100 compounds were further rescored and evaluated by CNN binding affinity against the PCSK9 protein using the Gnina docking software. (version 1.0) Due to the variation in docking scores generated by Glide, Smina, and Gnina for each ligand, we selected three molecules, namely amikacin, bestatin, and natamycin, based on their superior average docking scores and higher CNN affinities, for subsequent molecular dynamics (MD) simulation and post-MD simulation studies (Figure 2).

Prior to performing virtual screening of the FDA-approved library, the docking methodology was validated by removing the bound ligand from the protein structure (PDB id 6U2P) and redocking it into the active site of the enzyme. The RMSD value was calculated to assess the superimposition of the docked ligand with the co-crystallized ligand structure. The RMSD values were determined to be 0.601 Å and 0.424 Å for the Glide and Smina docking methodologies, respectively. These values below 2.0 Å indicate that both docking software generated correct poses for the co-crystallized ligand, validating the reliability of the docking process [45]. The average docking scores for the bound ligands, amikacin, natamycin, and bestatin, were found to be −8.90, −10.23, −9.98, and −9.78, respectively (Table 2). These compounds were further evaluated using Gnina docking software, which employs convolutional neural networks (CNNs) as a scoring function. The CNN model in Gnina provides predictions for both pose quality (CNNScore) and binding affinity (CNNaffinity). The CNNScore ranks the poses of the ligands, while the CNNaffinity predicts the affinity of each ligand in ‘pK’ units, where a pK value close to 6 indicates 1 μM affinity. The CNNScore of the bound ligands, amikacin, natamycin, and bestatin, were found to be 0.865, 0.58, 0.35, and 0.53, respectively, while their CNNaffinity values were 6.78, 5.45, 6.62, and 5.49, respectively.

Figure 3 provides a comprehensive depiction of the binding poses of all compounds at the allosteric site of the PCSK9 enzyme. It was observed that all compounds interacted with at least one amino acid residue, including Arg357 and Asp360 in the catalytic domain, and Arg458 and Arg476 in the C-terminal domains, as reported in the literature for their interaction with small molecules [35]. All compounds exhibited specific binding at the allosteric site, encompassing both the catalytic and C-terminal domains, except for bestatin, which, due to its smaller structure, was shifted more towards the catalytic domain. In comparison to the crystal structure’s bound ligand, which showed H-bond interactions with amino acid residues Pro331, Arg357, Asp360, Arg458, and Arg476, the amikacin, with its higher number of hydroxyl groups, displayed hydrogen bond interactions not only with these amino acids but also with Glu332, Cys358, Ile474, and Cys477 (Table 3). This increased the stability of amikacin in the allosteric site of the PCSK9 enzyme. Similarly, natamycin exhibited H-bond interactions with all four crucial amino acid residues in the catalytic and C-terminal domains. However, as bestatin was shifted towards the catalytic site, it did not display hydrogen bonding interactions with Arg476, but it formed H-bonds with amino acid residues Glu332, Val333, Arg357, Asp360, Arg458, and Trp461.

A molecular dynamics (MD) simulation study was conducted to investigate the behavior and stability of the PCSK9 enzyme alone and in combination with bound ligand, amikacin, natamycin, and bestatin. The MD simulations were carried out for 100 ns using the Desmond software (Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY, USA, 2020; Maestro–Desmond Interoperability Tools, Schrödinger, New York, NY, USA, 2020). Figure 4, Figure 5, Figure 6 and Figure 7 illustrate the key parameters utilized to assess the molecular stability of the docked compounds during the simulation. Various parameters, such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), and contact mapping of the ligand–protein complexes, were analyzed to confirm the molecular behavior and stability of the compounds. These parameters were monitored over the entire 100 ns simulation duration to gain insights into the conformational changes that occur when the ligands are complexed with the PCSK9 protein.

The RMSD is a widely used metric to measure the average distance between atoms in different protein structures, typically focusing on the backbone atoms. It provides a measure of the overall displacement of atoms in one frame relative to a reference frame. The RMSD analysis of the apo protein backbone revealed that the PCSK9 protein without any ligand maintained a stable conformation throughout the MD simulation run, with RMSD values ranging between 2.0 and 2.5 Å. Upon the introduction of the bound ligand structures, the RMSD values increased initially, reaching up to 3.0 Å within the first 45 ns of the simulation, and then stabilized within the range of 2.0 to 2.5 Å. The complex backbone RMSD values for the amikacin–enzyme complex showed slight fluctuations up to first 50 ns MD run, then stabilizing between 3.0 and 3.75 Å. In contrast, both the natamycin–protein and bestatin–protein complexes exhibited greater stability, with backbone RMSD values ranging from 2.4 to 3.0 Å. Apart from the amikacin–protein complex, the other two complexes displayed a relatively stable molecular behavior of the protein with slight fluctuations. This indicates that the natamycin–protein and bestatin–protein complexes achieved a high level of equilibrium and dynamic stability during the MD simulation run (Figure 4). Furthermore, the RMSD values of the tested compound–protein complexes closely matched those of the reference bound ligand–protein complex as well as apo protein, suggesting comparable behavioral stability of the complexes.

Additionally, RMSF analysis was performed on the Cα atoms of the amino acid residues in all five systems (Figure 5). RMSF provides valuable insights into local changes along the protein chain by measuring the fluctuation of each residue throughout the entire simulation. Notably, apart from the loop regions, the Cα atoms in the active site of the enzyme exhibited lower atomic fluctuations, suggesting minimal conformational changes. The RMSF analysis also confirmed the stability of each amino acid within the ligand–protein complexes, indicating that all compounds formed stable complexes with the PCSK9 enzyme with minimal structural alterations. Overall, RMSF analysis showed that the compound–protein complexes closely matched those of the reference bound ligand–protein complex as well as apo protein, suggesting comparable behavioral stability of the complexes. This indicates that the natamycin–protein and bestatin–protein complexes achieved a high level of equilibrium and dynamic stability during the MD simulation run. It is concluded that all compounds interacted with both the catalytic and C-terminal domains of the PCSK9 enzyme through common amino acid residues, including Pro331, Glu332, Val333, Arg357, Cys358, Asp360, Arg458, Trp461, Arg476, and Cys477. It can be concluded that all four compounds form stable complexes with the PCSK9 enzyme with minimal structural alterations. The compounds interact with key amino acid residues of the PCSK9 enzyme, including Asp360 and Arg458, which are known to be involved in the binding of PCSK9 to LDL-R.

We conducted further analysis to investigate the formation and stability of hydrogen bonds under dynamic conditions. Figure 6 represents a comprehensive overview of the individual occupancies of H-bond interactions that occurred during MD simulation between ligands and proteins. Consistent with the docking studies, the protein–ligand contacts plot revealed that all compounds interacted with both the catalytic and C-terminal domains of the PCSK9 enzyme through common amino acid residues, including Pro331, Glu332, Val333, Arg357, Cys358, Asp360, Arg458, Trp461, Arg476, and Cys477. When considering direct H-bond interactions without involving water bridges, the bound -ligand structure showed interactions with amino acid residues Pro331, Asp360, and Arg458, with an interaction fraction of over 0.5, indicating that these interactions were present for more than 50% of the MD simulation. Similarly, amikacin exhibited direct H-bond interactions with amino acid residues Glu332, Thr335, Cys358, Asp360, and Cys477, with an interaction fraction of over 0.5, indicating its precise binding between the catalytic and C-terminal domains of the allosteric site of PCSK9 enzyme. In line with the docking study, bestatin also displayed direct H-bond interactions with amino acid residues Val333, Thr335, Arg357, Asp360, Arg412, and Arg458 throughout the MD simulation. However, bestatin also formed a water bridge with Arg476, which remained stable for more than 50% of the MD simulation, suggesting it’s positioning in the center of the allosteric site contrary to the docking result. Conversely, natamycin demonstrated direct H-bond interactions with amino acid residues Thr335, Cys358, and Asp360, and an ionic interaction with Arg458. This indicates a slight shift of natamycin towards the catalytic site, as the direct H-bond interaction with Arg476 observed in the docking study was replaced by a water bridge with an interaction fraction of over 0.5.

Figure 7 represents a 2D pose diagram illustrating the ligand–protein contacts after the MD simulation run. The figure provides a clear visualization of the percentage of interactions between specific ligands and various surrounding amino acid residues at the allosteric site of the PCSK9 enzyme during the 100 ns MD simulation. Notably, the amino acid residue Arg458 in the C-terminal domain was found to play a crucial role in the formation of H-bond interactions and/or water bridges, which remained stable for over 70% of the MD simulation with all four compounds. It was also observed that Arg458 consistently acted as a hydrogen bond donor in its interactions with the compounds, predominantly interacting with the oxygen of the carbonyl group. Furthermore, all of the compounds contributed to the stability of the ligand–protein complex by forming H-bond interactions or water bridges with the amino acid residue Asp360, which acted as a hydrogen bond acceptor and primarily interacted with amino groups.

The MM-GBSA analysis revealed a significant correlation between the experimental and predicted binding affinity through Gibbs free energy calculations. To further validate the binding affinity of the ligands against the protein, post-MD simulation MM-GBSA calculations were conducted. The ΔG_bind, which estimates the binding free energy variation, was used as one of the key parameters for evaluating ligand–protein binding. The MM-GBSA ΔG_bind was calculated by comparing the energy difference between the bound and unbound states of the complexes. For the post-dynamic MM-GBSA analysis, 100 frames from the trajectories of each ligand–protein complex were selected at intervals of 10 ns. The binding free energy of these 100 systems in each simulation was computed using Prime software (version 4.0), and the mean values were reported (Figure 4). The calculated average ΔG_bind values ranged from −84.22 to −76.39 kcal/mol, indicating a high binding affinity of the compounds against the allosteric site of the PCSK9 enzyme. It is important to note that the MM-GBSA scoring function is optimized for predicting binding free energies in a congeneric series of molecules and, therefore, the absolute values calculated may not necessarily align with experimental binding affinities. However, the ranking of the ligands based on the calculated ΔG_bind is expected to reasonably correspond to the ranking based on experimental binding affinity, especially for congeneric series. The ΔG_bind values of all three drug compounds were observed near to bound ligand, which further indicates higher binding affinity against the PCSK9 enzyme. As ΔG_bind values represent approximate free energies of binding, a more negative value indicates stronger binding. Interestingly, all four systems exhibited negative ΔG_bind values, indicating strong stability of the systems. Additionally, bestatin demonstrated a better ΔG_bind value compared to the other compounds and the reference compound, suggesting a higher affinity towards the PCSK9 enzyme (Table 4).

The Prime MM-GBSA method utilizes the VSGB 2.0 solvation model, which incorporates the OPLS_2005 force field for both bonded and non-bonded terms, as well as a solvation term and several physics-based correction terms for hydrogen bonding, π-π interactions, self-contact interactions, and hydrophobic interactions. From the Table 5, it is apparent that the non-bonded terms, such as van der Waals (ΔG_vdW, −68.91 to −59.47 kcal/mol), Coulomb (ΔG_Coul, −60.67 to −11.16 kcal/mol), and lipophilic (ΔG_Lipo, −23.85 to −18.93 kcal/mol) energy terms, play a major favorable role in binding ligands to the PCSK9 enzyme. The physics-based energy term, hydrogen bonding (ΔG_Hbond, −4.82 to −2.83 kcal/mol), moderately supported the binding of the compounds. However, it was observed that the electrostatic solvation energy (ΔG_Solv, 16.67 to 59.47 kcal/mol) strongly disfavored the binding. Additionally, the covalent binding (ΔG_Coval, 3.96 to 6.86 kcal/mol) and π-π packing (ΔG_packing, −0.43 to 1.57 kcal/mol) energy terms also moderately disfavored the binding of compounds to the enzyme. The high negative values of ΔG_Coul and ΔG_Lipo indicate that the allosteric site of PCSK9 is lined with both polar and nonpolar residues, and interaction with these residues through hydrogen bonding increases the affinity. Furthermore, the high negative values of ΔG_vdW suggest that the tested compounds were well embedded within the allosteric sites upon binding. Moreover, it is evident that compounds with balanced polar and nonpolar structural properties, along with a defined number of hydrogen bond acceptors, hydrogen bond donors, and hydrophobic groups, such as bestatin, exhibit good binding affinity for the PCSK9 enzyme compared to more polar compounds. Overall, MM-GBSA analysis revealed a significant correlation between the experimental and predicted binding affinity through Gibbs free energy calculations. The calculated average ΔG_bind values ranged from −84.22 to −76.39 kcal/mol, indicating a high binding affinity of the compounds against the allosteric site of the PCSK9 enzyme. Based on the results of the MM-GBSA analysis, it can be concluded that all four compounds have a high binding affinity towards the PCSK9 enzyme. These results also suggest that the compounds have the potential to inhibit the interaction between PCSK9 and LDL-R, leading to a reduction in LDL cholesterol levels.