Identification of PARP12 Inhibitors By Virtual Screening and Molecular Dynamics Simulations

Poly [adenosine diphosphate (ADP)-ribose] polymerases (PARPs) are members of a family of 17 enzymes that performs several fundamental cellular processes. Aberrant activity (mutation) in PARP12 has been linked to various diseases including inflammation, cardiovascular disease, and cancer. Herein, a large library of compounds (ZINC-FDA database) has been screened virtually to identify potential PARP12 inhibitor(s). The best compounds were selected on the basis of binding affinity scores and poses. Molecular dynamics (MD) simulation and binding free energy calculation (MMGBSA) were carried out to delineate the stability and dynamics of the resulting complexes. To this end, root means deviations, relative fluctuation, atomic gyration, compactness, covariance, residue-residue contact map, and free energy landscapes were studied. These studies have revealed that compounds ZINC03830332, ZINC03830554, and ZINC03831186 are promising agents against mutated PARP12.

1 Introduction

Poly [adenosine diphosphate (ADP)-ribose] polymerases (PARPs) are the nuclear enzymes that regulate fundamental cellular processes involving protein degradation and gene expression (Kim et al., 2011; Bai and Virag, 2012; Langelier et al., 2018). PARPs use nicotinamide adenine dinucleotide (NAD+) for the process of post-translation that modify substrate proteins with ADP-ribose, a vital process referred to as ADP ribosylation (Carter-O'Connell and Cohen, 2015; Griffiths et al., 2020). The majority of PARP family members are monoPARPs as they catalyze the transfer of mono-ADP-ribose (MAR) onto their substrates (MARylation) (Cohen, 2020). In contrast, polyPARPs attach polymers of poly-ADP-ribose (PAR) onto their substrates (PARylation) (Carter-O'Connell and Cohen, 2015; Lu et al., 2019). PARP1 is a well-known cancer target involved in DNA damage-induced cellular stress/genetic mutation/cytotoxic chemotherapy (Gozgit et al., 2021). There are four drugs readily available on the market, and other drugs are at their late-stage development (Schiewer et al., 2012; Schiewer and Knudsen, 2014; Green et al., 2015; Rudolph et al., 2021). PARP12 is also a mono-ART of the PARP family, which controls the regulation of cell survival and regrowth (Catara et al., 2017). The function of PARP12 is the mono-ADP-ribosyltransferase that mediates mono-ADP-ribosylation of target proteins, although the molecular mechanism responsible for cell survival is still unclear (Catara et al., 2017; Grimaldi et al., 2020). Recently, PARP12 involvement in intracellular membrane transport at this point deserves much attention and alternative effective strategies (Gozgit et al., 2021). Several studies also suggested that PARP12 is a novel target that reduces breast cancer resistance to genotoxic stress (Ke et al., 2019). Thus, PARP12 is marked as a potential target that must be utilized further to identify and design potential inhibitors/vaccine candidates at the molecular and genetic levels. The authors identify PARP1 potent targets and various computational studies that have been published to utilize the ZINC-FDA library and traditional Chinese medicine (Costantino et al., 2001; Chen et al., 2014; Maksimainen et al., 2021; Sahin and Durdagi, 2021). Similarly, there are other recently published theoretical studies on PARP2 to find novel benzimidazole derivatives (Venugopal and Chakraborty, 2021).

Bioinformatics tools and databases are wonderful approaches to identify the novel drugs and vaccine candidates as an alternative strategy. Taking the lead from here, the current study aimed to utilize the available small-molecule ZINC database, which is a publicly available and accessible database. We used a comprehensive in silico approach involving literature mining, virtual screening based molecular docking, pharmacological analysis, and large-scale molecular dynamics (MD) simulations to identify potential PARP12 inhibitors from a pool of ZINC-FDA. We analyzed 3,100 chemically diverse pharmacologically active compounds (Sterling and Irwin, 2015). Based on the comparative analysis of docking scores, binding affinity, and energies, we selected four compounds, namely, ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189. Consequently, the compounds were analyzed for ADME/T properties and were potential drug-like candidates. Furthermore, we analyzed the conformational stability of the docked complexes using MD simulations with the help of various parameters such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), free energy landscapes (FEL), H-bond monitoring, principal component analysis (PCA), and density distribution map (Ahamad et al., 2021b). Based on MD simulation results, two compounds (ZINC03830332 and ZINC03831186) were identified as PARP12 inhibitors.

2 Materials and Methods

The work was carried out on a Dell workstation, 256 GB RAM, and NVIDIA graphics card GeForce GTX 1080Ti GDDR5 11 GB. The workstation was installed with several bioinformatics software applications such as GROningen MAchine for Chemical Simulations (GROMACS), XMGrace (GROMACS results analysis), and Schrödinger Molecular Mechanics Generalized Born Surface Area (MMGBSA) and 2D plots. The detailed workflow used in this study is represented in Figure 1.

FIGURE 1

FIGURE 1. Optimized protocol employed in this work.

2.1 Protein Preparation and Grid Generation for Docking Studies

The human PARP12 with four-point mutations (Val570, Gly573, His577, and Phe607) has been selected in this study. The PARP12 has a 701 amino acid length, and the crystal structure (PDB = 6V3W) is available with the PARP catalytic domain (496–680 amino acids). To identify the potential surface structural pockets, their shape and volume, internal cavities of the protein, and surface areas, 6V3W was given as the input (Gozgit et al., 2021). The active site and the interactive residues were identified using the online tools CASTp and PDBsum and verified using the literature (Laskowski, 2001; Dundas et al., 2006; Gozgit et al., 2021). The ligands (FDA-ZINC compounds) were prepared using AutoDock tools (ADT) by converting the PDB file into a PDBQT format (Huey and Morris, 2008; Ahamad et al., 2019). The crystal structure of PARP12 was prepared by using the ADT protein preparation wizard. Polar hydrogen (H) bonds and missing H-atoms were added while the water molecules and heteroatoms were deleted from the 3D structure. Energy minimization was performed with a default constraint of 0.3 Å root mean square (RMS), and charges were assigned. Finally, the clean structure was saved as a PDBQT file, and docking was performed by AutoDockVina (version 1.1.2) software (Trott and Olson, 2010). The size of the grid box selected was 106 Å × 108 Å × 106 Å, respectively, and generated around the centroid of the compounds and PARP12 complex. Default parameter settings were used during docking. The 3D structure of the docked complexes was analyzed with the help of the PyMOL visualization tool (DeLano, 2002).

2.2 Structure-Based Virtual Screening

The PDBQT formatted ligands (ZINC-FDA compounds) were screened against PARP12. A total of 3,100 molecules were retrieved from the ZINC-FDA database to identify potential PARP12 inhibitors (Sterling and Irwin, 2015). The compounds were ranked on the basis of their binding energy scores and docking interaction poses, and the top four compounds were selected.

2.3 Molecular Dynamics (MD) Simulations

MD simulations were carried out for the best-docked complexes using GROMACS (version 5.18.3) to determine the PARP12 enzyme behavior in the presence of water (Van Der Spoel et al., 2005; Abraham et al., 2015). The topology of PARP12 was created via the GROMOS9643a1 force field (Pronk et al., 2013). The PRODRG server was used to generate molecular topologies and coordinate files (Schuttelkopf and van Aalten, 2004). All the systems were solvated using a simple point-charge model (SPC/E) in a cubic box (Mark and Nilsson, 2001). To neutralize the system, 0.15 M counter ions (Na⁺ and Cl⁻) were added. The energy minimization of all the neutralized systems was performed using the steepest descent and conjugate gradients (50,000 steps for each). The volume (NVT) regulation and pressure (NPT) were run for system equilibration. Steepest descent followed by conjugate gradient algorithms was utilized on enzyme-ligand complexes. The NVT ensemble at a constant temperature of 300 K and the constant pressure of 1 bar was employed. The SHAKE algorithm was used to confine the H-atoms at their equilibrium distances and periodic boundary conditions.

Moreover, the particle mesh Ewald (PME) method was used to define long-range electrostatic forces (Darden et al., 1993; Lee et al., 2016; Wang et al., 2016). The cut-offs for van der Waals and Columbic interactions were 1.0 nm. The LINC algorithm was used to constrain the bonds and angles (Hess et al., 1997). Finally, the production runs were performed for the period of 500 ns The energy, velocity, and trajectory were updated at the time interval of 10 ps The MD simulation analyses were performed, and trajectoris were found by GROMACS utilities and MDTraj-based Python scripts (McGibbon et al., 2015). The Cα-atom deviations of PARP12 and complexes were utilized to calculate using RMSD and RMSF for relative fluctuations of each amino acid. The compactness was measured by Rg, while the SASA was used to know the electrostatic contributions of molecular solvation. PCA is one of the best techniques that help reduce complexity to extract the intensive motions in MD simulation analysis. The matrix was formed for the MD trajectory after excluding rotational and translational movements. So the essential dynamics protocol was implemented to calculate the eigenvectors and eigenvalues and their projections along with the first two principal components (PCs). The diagonalized covariance matrix, eigenvectors, and eigenvalues were identified. The eigenvector of the matrix gives the multidimensional space and the displacement of atoms in the molecule in each direction. This analysis processed the essential subspace built on knowing each atom’s movements, which were plotted by Cartesian trajectory coordinates using GROMACS utilities (Ahamad et al., 2021a).

2.4 MMGBSA Calculation

The MMGBSA (ΔG_bind) was calculated by the Schrödinger Prime module to compute the ligand binding energies. The prime-MMGBSA was used for rescoring the docked poses.

ΔG is calculated by the following equation:

$\mathrm{\Delta }Gbind = E\_complex \left(minimized\right) - E\_ligand\left(minimized\right)- E\_receptor\left(minimized\right),$

where ΔG_bind specifies the binding free energy, E__complex is the input for the energy minimization of the protein-ligand complexes, E__receptor for the free protein, and E__ligand for free ligands. The sum of molecular mechanical minimized energies including van der Waals interaction, internal energies, electrostatic energies, and solvation free energies were also calculated for the docked conformations.

3 Results

Despite the fact that PARP12 is an intriguing target for treating various ailments, no computational studies on PARP12 have been carried out yet. Based on this notion, we screened 3,100 biologically active molecules from the FDA-approved drug database to find out lead(s) against the PARP12 enzyme. Out of these, four compounds (ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189, Supplementary Figure S1) were selected on the basis of their docking score and interaction with the receptor. Moreover, the four best hits were analyzed for Lipinski and ADME/T properties. All-atom MD simulations were further performed on the mutant PARP12 and four-drug complexes. The results of the studies are discussed herein.

3.1 Molecular Docking

Molecular docking is a powerful tool for lead discovery and to underpin molecular interactions. Figure 2A depicts the binding pockets (A* and B* encircled by black dotted lines) of the PARP12 enzyme and binding poses of the docked ligands. Herein, we have shortlisted four compounds that carried the best binding free energies, as well as polar and non-polar interactions. Molecular docking allows prediction of the preferred pose and binding orientations of a ligand inside a receptor binding site corresponding to the best energy score. The top-scoring shortlisted four compounds were docked in the A* binding pocket. However, B* binding pocket interaction compounds have fewer binding energies. Figures 2B,C and Supplementary Figure S1 show the molecular interaction between the four hits and receptor. Table 1 collects the results of the molecular docking study of four shortlisted ZINC-FDA compounds. It is clear from the figures that ligands interacted with the receptor via various polar and non-polar modes (Figures 2D,E and Supplementary Figure S1). It can be seen that ZINC03830332 formed six H-bonds with Phe563, His564, Asn578, Trp581, Tyr596, and Ser604 residues of the receptor with a docking score of −12.40 kcal/mol. On the other hand, the docking score for the PARP12-ZINC03831186 complex was −11.40 kcal/mol which formed two H-bonds with the receptors His564 and Ser595. The two other ligands ZINC03830554 and ZINC03831189 exhibited relatively higher docking scores (−11.80 kcal/mol and −12.11 kcal/mol, respectively) and interacted via seven (Phe563, His564, His577, Asp580, Tyr596, Ser604, and Lys666) and two (Gly594 and Ser604) H-bonds. The 2D plot analyses of PARP12-ligand complexes revealed that the complexes were also stabilized through hydrophobic interactions from nearby residues.

FIGURE 2

FIGURE 2. (A) Close-up view of the binding pockets (A* and B* encircled by black dotted lines) of the PARP12 enzyme. H-bond interactions between (B) ZINC03830332 (white stick), (C) ZINC03831186 (yellow stick), and PARP12 receptor (in green and blue color cartoon models). (D,E) 2D plot of the complexes.

TABLE 1

TABLE 1. Average RMSD, Rg, SASA, and MMGBSA values of the mutant receptor and complexes. Molecular docking results of four shortlisted ZINC-FDA compounds.

3.2 Physicochemical and ADME/T Studies

The shortlisted best four hits, namely, ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 compounds were used to calculate the pharmacological properties such as MW, dipole, SASA, FOSA, FISA, PISA, WPSA, volume, DonorHB, and AccptHB. We found that all the four molecules obeyed the rule of three (Ro3) (Kö;ster et al., 2011). Preferably, an orally active drug/compound should not cross more than one violation, but various exceptions are available (Lipinski, 2004). Two compounds (ZINC03830332 and ZINC03830554) violated one rule increase in MW, which is still acceptable. The results of the RO3 study revealed that the molecular weight of the compound was in the range of 422–693, with 4–10 H-bond acceptor and 0–14 H-bond donating units. The in silico absorption, distribution, metabolism, excretion (ADME), and toxicity (T) are the most widely used techniques in rational drug discovery as it gives a fair idea about the drug candidateship. In this context, we investigated the ADME/T and RO3 properties of the best four hits, and the results are compiled in Supplementary Table S1. The QikProp module output data generated for the compounds indicated that the blood–brain barrier (BBB) was in the range of −0.38 to −4.28 with QPlogPo/w between −0.55 and 5.73, QPlogKhsa between −0.63 and 1.14, SASA between 714.08 and 1,083.70, and log Kp value ranging from −2.41 to −10.90. Overall, the four hits showed fair enough properties and ability needed for a drug to possess (Lee et al., 2003; Lipinski, 2016; Nazarbahjat et al., 2016; Cordeiro and Kachroo, 2020).

3.3 Molecular Dynamics Simulation Analysis

MD simulations have the ability to uncover various dynamic interactions between a ligand and receptor, their interaction mechanism, and stability. MD simulation was carried out at the 0 to 500 nanosecond (ns) time scale for each system with a total of 2,500 ns. The MD simulations of the mutant systems of PARP12 and ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 complexes were utilized to measure the structural changes and parameters such as RMSD, RMSF, Rg, SASA, H-bonds, and Gibbs free energy, i.e., PCA calculation, covariance matrix, and density distribution map. The results of the findings are discussed as follows.

3.3.1 RMS Deviation and Principal Component Analysis

The MD simulation study result dictates the multidimensional data’s relative distance to reduce the dimensional space. Atom-positional RMSD value C-alpha (C-α) atoms were calculated for the mutant system and complexes (Figure 3A). The MD simulation of mutant PARP12 systems and complexes with ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 compounds provided average RMSD values of ∼0.43 nm, ∼0.32 nm, ∼0.34 nm, ∼0.32 nm, and ∼0.37 nm, respectively (Figure 3B). On the other hand, the mutant system of PARP12 attained equilibrium in 75 ns with a relatively higher average RMSD value due to the mutations in the helix and loop region. However, in the case of complexes, almost all the systems attained equilibrium within 30 ns and the steady-state in all the complexes except ZINC03831189, which possess a slightly greater RMSD value (∼ 0.37 nm). Overall, the complexes involving ZINC03830332, ZINC03831186, and ZINC03830554 possessed a stable binding with PARP12.

FIGURE 3

FIGURE 3. (A) Mutation induced by PARP12 structures, (B) comparison of the RMSD plot, and (C) dynamic motion of projection in the eigenvector 1 vs. 2 and the plot two generated for the MD complexes showing conformational space of Cα-atoms.

The 2D projection of trajectories on eigenvectors 1 and 2 is the part of essential dynamics (ED) and reflect an overall expansion of the structural dynamics. The large-scale average motion indicated higher atomic fluctuations, and its flexibility measures the atomic mobility of the mutant and complexes. The eigenvector values of the docked complexes ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 were calculated. We found that the essential subspace for the mutant PARP12 of eigenvector 1 in the range from −1.5 to −2.5 nm, vs. eigenvector 2 in the range from −1.5 to 3 nm of a larger cluster for the mutant PARP12 (Figure 3C). Similarly, we extensively investigated the global motion for the complex of the ZINC03830332 compound. We found that the eigenvector 1 in the range from −1.5 nm to −2 nm vs. eigenvector 2 in the range from −1 to 1.5 nm occupied the smallest subspace throughout the MD simulation. However, compounds ZINC03830554 and ZINC03831186 had eigenvector 1 in the range −2 to 2 nm vs. eigenvector 2 in the range −1.5 to 2.5 nm and eigenvector 1 in the range −1 to 3 nm vs. eigenvector 2 in the range −1 to 1 nm, respectively. The compound ZINC03831189 had eigenvector 1 in the range −3 to 1.5 nm vs. eigenvector 2 in the range −1.5 to 2 nm. Out of four complexes, we found a more restricted subspace that leads to well-defined internal motion with ZINC03830332, ZINC03830554, and ZINC03831186 complexes.

The rigid regions and flexibility could be identified by RMSF analysis. The plot’s negative drift indicates the increased movement of the Ca atoms of relative fluctuations. It was found that the compound ZINC03830332 showed stable positive amino acid fluctuations at F563-T566, V570, H577-W581, V583-596, Q613, I660-I663, K666, and Y673. A similar trend was found for the other two compounds with a high degree of flexibility (Figures 4A–D).

FIGURE 4

FIGURE 4. Root mean square fluctuations (RMSFs) of mutant systems of PARP12 and (A) ZINC03830332, (B) ZINC03830554, (C) ZINC03831186, and (D) ZINC03831189.

To gain more insight into the flexibility, the Rg plot, which measures the mass-weighted RMS distance of a cluster of atoms from their center of mass, was determined. The Rg analysis of C-α atoms of ligand–receptor docked complexes average values was calculated between 0 and 500 ns. Figure 5A illustrates the Rg plot, and the results indicated a compact conformation of the complexes, as compared to the native enzyme. For instance, average values ∼1.57 nm, ∼1.59 nm, ∼1.58 nm, ∼1.59 nm, and ∼1.58 nm were found for PARP12 and ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 complexes, respectively (Table 1). Overall, the result suggested mutant expanded conformation while compact conformation for the complexes.

FIGURE 5

FIGURE 5. (A) Rg and (B) SASA fluctuations per residue variation plot analysis of the mutant system and in the presence of compounds with PARP12 at total simulation time.

To understand the behavior of the PARP12 mutant system and all the complexes, we also performed a SASA analysis which revealed the PARP12 surface that can be accessed in the solvent (water). The average SASA value of the mutant PARP12, complexes with ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 were found to be 101.70, 102.17, 104.02, 103.26, and 103.24 nm², respectively, (Figure 5B and Table 1).

3.3.2 Hydrogen Bond Vetting

Hydrogen bonding plays a significant role in knowing the complexes’ stability. To determine the binding affinity of the selected compounds toward the PARP12 receptor, MD trajectories were analyzed, and an H-bond investigation was carried out to determine the total number of bond formations. The ZINC03830332 and ZINC03830554 formed 0–14 H-bonds consistent for both the compounds’ complexes with PARP12 throughout the 500 ns MD simulation. Similarly, the compound ZINC03831186 complex has quite less 0–5 number of H-bonds, whereas ZINC03831189 docked complexes showed H-bonds 0–4 (supplementary Figures S3A–D). The intermolecular H-bond examination revealed that the complex containing ZINC03830332 possessed a minimum of fourteen H-bonds during the MD simulation. The number of H-bond interactions also increased in the case of ZINC03830554 during the simulation. Contrarily, complexes containing ZINC03831186 and ZINC03831189 were unable to increase the number of H-bonds; however, both compounds showed consistency in the H-bonds’ profile.

3.3.3 Free Energy Landscape Analysis

To monitor the distinct conformation of binding, FEL analysis was performed. The principal components (PCs) dictate the most dominant internal modes of motion of a corresponding system. In the ligand binding or complex formation in MD, through PCA, we found the most probable and dominant conformation changes of the protein during binding. The first PC1 and second PC2 utilized projected eigenvectors generated based on PCA (Stein et al., 2006). The contour map of the FEL derived conformational changes of the mutant system and complexes of PC1 and PC2 eigenvectors (Figures 6A–E). The Gibbs free energy landscape shows the global energy minima states (Ali et al., 2019; Ahamad et al., 2021b). The energetically favored PARP12 complex conformation is indicated by dark blue spots, while the unfavorable conformations are indicated by yellow spots. The values of FEL ranged from 0 kJ/mol to 19.1 kJ/mol, 18.2 kJ/mol, 18.2 kJ/mol, 18.1 kJ/mol, and 16.7 kJ/mol for the mutant PARP12 and ZINC03830332, ZINC03831186, ZINC03831189, and ZINC03830554 docked complexes, respectively. The FEL results demonstrated several well-defined energy minima, corresponding to the metastable conformational states. From Figure 6, it is clear that the main free energy wells in the global free energy minima area were altered when selected ZINC-FDA compounds were bound with PARP12. The PC analysis of the existing results confirmed that the FEL analysis is complimented to our previous RMSD finding.

FIGURE 6

FIGURE 6. FEL direction of motion and magnitude analysis of the mutant PARP12 (A) and complexes with ZINC03830332 (B), ZINC03830554 (C), ZINC03831186, and (D) ZINC03831189 (E) throughout MD simulations. The color bar denotes the relative free energy value between 0 and 19.1 kcal mol⁻¹.

The g_covar module by PCA monitors the conformational changes of C-α mutant and complexes. In this analysis, every single atom’s collective motion along with directions provides a covariance matrix that uses atomic motion in Cartesian coordinate space. There are two types of motion: 1) correlated and 2) anticorrelated motion. The covariance matrix is constituted by two intense colors; red-colored region, which means atoms are moving together, and the blue-colored region, which means atoms are moving in the opposite direction. The matrix depicts that the reduction of negative correlation is largely exhibited with the ZINC03830332 complex. The covariance matrix results demonstrated that all the complexes were tolerable atomic displacement of ZINC03830332, ZINC03830554, and ZINC03831186 (Supplementary Figure S4). The increment in negative correlation was observed in the case of mutant, while reduction of both negative and positive motion was noted in the case of the ZINC03831189 complex. The overall positional fluctuations on Cα-atoms of all the docked complexes revealed that ZINC03830332, ZINC03830554, and ZINC03831186 were stable with steady atomic displacement and amplitude toward PARP12.

The atomic density distribution was measured using mutant PARP12 and complexes ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189 ranging from 1.35 nm⁻³, 1.89 nm⁻³, 2.3 nm⁻³, 1.88 nm⁻³, and 1.38 nm⁻³, respectively (Supplementary Figure S5). The results showed the stable density area for ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189. The density area of each atom of the complexes ZINC03830332, ZINC03830554, and ZINC03831186 was examined to be in stable distribution; the energy wells displayed three compounds signifying the binding at the same binding site and the stabilized PARP12 enzyme.

4 Discussion

Poly [adenosine diphosphate (ADP)-ribose] polymerases (PARPs) are one of the attractive and promising targets to treat several life-threatening diseases including cancer (Kim et al., 2011; Bai and Virág, 2012; Langelier et al., 2018). The present study was undertaken to identify promising lead(s) against the mutated PARP12 enzyme. We carried out an extensive in silico study, including virtual screening, drug-likeliness, MD simulations, and MMGBSA analyses. To the best of our knowledge, no computational studies on PARP12 have been carried out before. To this end, virtual screening of the ZINC-FDA dataset of 3,100 compounds was performed.

Molecular docking studies revealed the interaction of the compounds within the active site of receptor via multiple amino acids. For instance, other than H-bonds, the compound ZINC03830332 was also stabilized by Gly564, Thr566, Val570, Gly573, Ile574, Cys575, His577, Phe579, Asp580, Trp581, Arg582, Val583, Cys584, Gly585, Val586, Thr589, Gly594, Ser595, Phe597, Ala598, Asp600, Tyr603, His605, Phe607, Ser659, Ile660, and Phe661 residues. Furthermore, the best four hits (ZINC03830332, ZINC03830554, ZINC03831186, and ZINC03831189) were selected based on their docking score and ADME/T profile used for further studies. The results of the drug-likeness study revealed that the selected compounds obeyed Lipinski’s properties very well and have a high propensity to be developed as a drug.

The MD simulation study revealed that mutation in the enzyme altered the enzyme dynamics as revealed by higher fluctuation throughout the MD simulation. However, upon interaction with the ligand, PARP12 gets stabilized as indicated by smooth fluctuations and showed steady-state throughout the MD simulation. The results of the Rg analysis indicated stable folding behavior of PARP12 after binding with the compounds. The results suggested that the Rg values of complexes with four ZINC-FDA compounds stayed strongly bound to the binding site and maintained stable and enhanced compactness of the PARP12 structure better than the mutant system throughout the MD simulation. The H-bonding profile indicated that compounds ZINC03830332 and ZINC03830554 formed stable complexes throughout the MD simulation. The FEL demonstrated several energy minima corresponding to the metastable conformational states. This analysis observed that the main free energy wells in the global free energy minima area were altered by the ZINC-FDA compounds’ complex bound with PARP12.

In addition, SASA study results showed that all the compounds achieved stable hydrophobic contacts, establishing the maximum region of PARP12 complexes accessible to the solvent molecules. The docked complexes exhibited higher average SASA profiles and stability than the native ligands.

5 Conclusion

In the present study, we used pharmacoinformatic techniques to identify potential inhibitors of mutant PARP12 receptors. Summarily, we conducted virtual screening, drug-likeliness, MD simulations, and MMGBSA to identify the best lead compounds amongst the ZINC-FDA compound library. Virtual screening–based molecular docking was carried out to select the best compounds capable of binding to the mutant PARP12 binding site, which helped shortlist the four compounds. Furthermore, Lipinski’s rule of three and ADME/T properties checked the best four compounds with various parameters to get probable drug candidates. The RMSD, RMSF, H-bonds, essential dynamics, and FEL wells analysis dictated that out of 3,100 compounds, only three ZINC-FDA compounds were the most stable, which perfectly bind to the binding pocket PARP12. Overall, the results indicated that the compounds ZINC03830332, ZINC03830554, and ZINC03831186 were better than others. These promising candidates could be used in the future and should be further explored clinically.

Data Availability Statement

The publicly available datasets were analyzed in this study. These data can be found at: https://zinc.docking.org/.

Author Contributions

Conceptualization: TA, MS and DY; data curation: TA, IA, and MS; formal analysis: AA, AGA and DY; funding acquisition: TA and MS; investigation: TA, MS and DY; methodology: TA, MYA and MA; project administration: TA, MS, and DY; resources: TA, MS and DY. Analyzed the data: SA; Writing and corrections: JAS.

Funding

This research was funded by the Scientific Research Deanship at King Khalid University and the Ministry of Education in KSA, grant number IFP-KKU-2020/14.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.847499/full#supplementary-material

References

Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 1-2, 19–25. doi:10.1016/j.softx.2015.06.001

CrossRef Full Text | Google Scholar

Ahamad, S., Hema, K., Kumar, V., and Gupta, D. (2021a). The Structural, Functional, and Dynamic Effect of Tau Tubulin Kinase1 upon a Mutation: A Neuro-Degenerative Hotspot. J. Cell Biochem. 122, 1653–1664. doi:10.1002/jcb.30112

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahamad, S., Islam, A., Ahmad, F., Dwivedi, N., and Hassan, M. I. (2019). 2/3D-QSAR, Molecular Docking and MD Simulation Studies of FtsZ Protein Targeting Benzimidazoles Derivatives. Comput. Biol. Chem. 78, 398–413. doi:10.1016/j.compbiolchem.2018.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahamad, S., Kanipakam, H., Kumar, V., and Gupta, D. (2021b). A Molecular Journey to Check the Conformational Dynamics of Tau Tubulin Kinase 2 Mutations Associated with Alzheimer's Disease. RSC Adv. 11, 1320–1331. doi:10.1039/d0ra07659g

PubMed Abstract | CrossRef Full Text | Google Scholar

Ali, S., Khan, F. I., Mohammad, T., Lan, D., Hassan, M. I., and Wang, Y. (2019). Identification and Evaluation of Inhibitors of Lipase from Malassezia Restricta Using Virtual High-Throughput Screening and Molecular Dynamics Studies. Int. J. Mol. Sci. 20, 884. doi:10.3390/ijms20040884

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, P., and Virág, L. (2012). Role of poly(ADP-Ribose) Polymerases in the Regulation of Inflammatory Processes. FEBS Lett. 586, 3771–3777. doi:10.1016/j.febslet.2012.09.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter-O'connell, I., and Cohen, M. S. (2015). Identifying Direct Protein Targets of Poly-ADP-Ribose Polymerases (PARPs) Using Engineered PARP Variants-Orthogonal Nicotinamide Adenine Dinucleotide (NAD+) Analog Pairs. Curr. Protoc. Chem. Biol. 7, 121–139. doi:10.1002/9780470559277.ch140259

PubMed Abstract | CrossRef Full Text | Google Scholar

Catara, G., Grimaldi, G., Schembri, L., Spano, D., Turacchio, G., Lo Monte, M., et al. (2017). PARP1-produced Poly-ADP-Ribose Causes the PARP12 Translocation to Stress Granules and Impairment of Golgi Complex Functions. Sci. Rep. 7, 14035. doi:10.1038/s41598-017-14156-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, K.-C., Sun, M.-F., and Chen, C. Y.-C. (2014). in Silico Investigation of Potential PARP-1 Inhibitors from Traditional Chinese Medicine. Evid. Based Complement. Alternat. Med. 2014:917605. doi:10.1155/2014/917605

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, M. S. (2020). Interplay between Compartmentalized NAD+ Synthesis and Consumption: a Focus on the PARP Family. Genes. Dev. 34, 254–262. doi:10.1101/gad.335109.119

PubMed Abstract | CrossRef Full Text | Google Scholar

Cordeiro, R., and Kachroo, M. (2020). Synthesis and Biological Evaluation of Anti-tubercular Activity of Schiff Bases of 2-Amino Thiazoles. Bioorg Med. Chem. Lett. 30, 127655. doi:10.1016/j.bmcl.2020.127655

PubMed Abstract | CrossRef Full Text | Google Scholar

Costantino, G., Macchiarulo, A., Camaioni, E., and Pellicciari, R. (2001). Modeling of poly(ADP-Ribose)polymerase (PARP) Inhibitors. Docking of Ligands and Quantitative Structure-Activity Relationship Analysis. J. Med. Chem. 44, 3786–3794. doi:10.1021/jm010116l

PubMed Abstract | CrossRef Full Text | Google Scholar

Darden, T., York, D., and Pedersen, L. (1993). Particle Mesh Ewald: AnN⋅Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 98, 10089–10092. doi:10.1063/1.464397

CrossRef Full Text | Google Scholar

Delano, W. L. (2002). Pymol: An open-source molecular graphics tool. CCP4 Newsl. Protein Crystallogr 40 (1), 82-92.

Dundas, J., Ouyang, Z., Tseng, J., Binkowski, A., Turpaz, Y., and Liang, J. (2006). CASTp: Computed Atlas of Surface Topography of Proteins with Structural and Topographical Mapping of Functionally Annotated Residues. Nucleic Acids Res. 34, W116–W118. doi:10.1093/nar/gkl282

PubMed Abstract | CrossRef Full Text | Google Scholar

Gozgit, J. M., Vasbinder, M. M., Abo, R. P., Kunii, K., Kuplast-Barr, K. G., Gui, B., et al. (2021a). PARP7 Negatively Regulates the Type I Interferon Response in Cancer Cells and its Inhibition Triggers Antitumor Immunity. Cancer Cell 39, 1214–e10. e1210. doi:10.1016/j.ccell.2021.06.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, A. R., Caracappa, D., Benhasouna, A. A., Alshareeda, A., Nolan, C. C., Macmillan, R. D., et al. (2015). Biological and Clinical Significance of PARP1 Protein Expression in Breast Cancer. Breast Cancer Res. Treat. 149, 353–362. doi:10.1007/s10549-014-3230-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffiths, H. B. S., Williams, C., King, S. J., and Allison, S. J. (2020). Nicotinamide Adenine Dinucleotide (NAD+): Essential Redox Metabolite, Co-substrate and an Anti-cancer and Anti-ageing Therapeutic Target. Biochem. Soc. Trans. 48, 733–744. doi:10.1042/BST20190033

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimaldi, G., Schembri, L., Monte, M. L., Spano, D., Di Martino, R., Beccari, A. R., et al. (2020). PARP12-catalyzed Mono-ADP-Ribosylation of Golgin-97 Controls the Transport of E-Cadherin. bioRxiv.

Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). LINCS: a Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 18, 1463–1472. doi:10.1002/(sici)1096-987x(199709)18:12<1463::aid-jcc4>3.0.co;2-h

CrossRef Full Text | Google Scholar

Huey, R., and Morris, G. M. (2008). Using AutoDock 4 with AutoDocktools: A Tutorial. USA: The Scripps Research Institute. 8, 54–56.

Google Scholar

Ke, Y., Wang, C., Zhang, J., Zhong, X., Wang, R., Zeng, X., et al. (2019). The Role of PARPs in Inflammation-And Metabolic-Related Diseases: Molecular Mechanisms and beyond. Cells 8. 1047. doi:10.3390/cells8091047

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, M. S., An, C. H., Kim, S. S., Yoo, N. J., and Lee, S. H. (2011). Frameshift Mutations of Poly(adenosine Diphosphate-Ribose) Polymerase Genes in Gastric and Colorectal Cancers with Microsatellite Instability. Hum. Pathol. 42, 1289–1296. doi:10.1016/j.humpath.2010.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Köster, H., Craan, T., Brass, S., Herhaus, C., Zentgraf, M., Neumann, L., et al. (2011). A Small Nonrule of 3 Compatible Fragment Library Provides High Hit Rate of Endothiapepsin Crystal Structures with Various Fragment Chemotypes. J. Med. Chem. 54, 7784–7796. doi:10.1021/jm200642w

PubMed Abstract | CrossRef Full Text | Google Scholar

Langelier, M. F., Eisemann, T., Riccio, A. A., and Pascal, J. M. (2018). PARP Family Enzymes: Regulation and Catalysis of the poly(ADP-Ribose) Posttranslational Modification. Curr. Opin. Struct. Biol. 53, 187–198. doi:10.1016/j.sbi.2018.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Laskowski, R. A. (2001). PDBsum: Summaries and Analyses of PDB Structures. Nucleic Acids Res. 29, 221–222. doi:10.1093/nar/29.1.221

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., et al. (2016). CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field. J. Chem. Theory Comput. 12, 405–413. doi:10.1021/acs.jctc.5b00935

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Lee, I., Kim, H., Chang, G., Chung, J., and No, K. (2003). “The PreADME Approach: Web-Based Program for Rapid Prediction of Physico-Chemical, Drug Absorption and Drug-like Properties. Blackwell Publishing: Boston, MA, USA, 2003; pp. 418–420.

Google Scholar

Lipinski, C. A. (2004). Lead- and Drug-like Compounds: the Rule-Of-Five Revolution. Drug Discov. Today Technol. 1, 337–341. doi:10.1016/j.ddtec.2004.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Lipinski, C. A. (2016). Rule of Five in 2015 and beyond: Target and Ligand Structural Limitations, Ligand Chemistry Structure and Drug Discovery Project Decisions. Adv. Drug Deliv. Rev. 101, 34–41. doi:10.1016/j.addr.2016.04.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, A. Z., Abo, R., Ren, Y., Gui, B., Mo, J. R., Blackwell, D., et al. (2019). Enabling Drug Discovery for the PARP Protein Family through the Detection of Mono-ADP-Ribosylation. Biochem. Pharmacol. 167, 97–106. doi:10.1016/j.bcp.2019.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Maksimainen, M. M., Murthy, S., Sowa, S. T., Galera-Prat, A., Rolina, E., Heiskanen, J. P., et al. (2021). Analogs of TIQ-A as Inhibitors of Human Mono-ADP-Ribosylating PARPs. Bioorg. Med. Chem. 52, 116511. doi:10.1016/j.bmc.2021.116511

PubMed Abstract | CrossRef Full Text | Google Scholar

Mark, P., and Nilsson, L. (2001). Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 105, 9954–9960. doi:10.1021/jp003020w

CrossRef Full Text | Google Scholar

Mcgibbon, R. T., Beauchamp, K. A., Harrigan, M. P., Klein, C., Swails, J. M., Hernández, C. X., et al. (2015). MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophys. J. 109, 1528–1532. doi:10.1016/j.bpj.2015.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Nazarbahjat, N., Ariffin, A., Abdullah, Z., Abdulla, M. A., Shia, J. K. S., and Leong, K. H. (2016). Synthesis, Characterization, Drug-Likeness Properties and Determination of the In Vitro Antioxidant and Cytotoxic Activities of New 1,3,4-oxadiazole Derivatives. Med. Chem. Res. 25, 2015–2029. doi:10.1007/s00044-016-1660-5

CrossRef Full Text | Google Scholar

Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., et al. (2013). GROMACS 4.5: a High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 29, 845–854. doi:10.1093/bioinformatics/btt055

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudolph, J., Roberts, G., and Luger, K. (2021). Histone Parylation Factor 1 Contributes to the Inhibition of PARP1 by Cancer Drugs. Nat. Commun. 12, 736. doi:10.1038/s41467-021-20998-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahin, K., and Durdagi, S. (2021). Identifying New Piperazine-Based PARP1 Inhibitors Using Text Mining and Integrated Molecular Modeling Approaches. J. Biomol. Struct. Dyn. 39, 681–690. doi:10.1080/07391102.2020.1715262

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiewer, M. J., Goodwin, J. F., Han, S., Brenner, J. C., Augello, M. A., Dean, J. L., et al. (2012). Dual Roles of PARP-1 Promote Cancer Growth and Progression. Cancer Discov. 2, 1134–1149. doi:10.1158/2159-8290.CD-12-0120

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiewer, M. J., and Knudsen, K. E. (2014). Transcriptional Roles of PARP1 in Cancer. Mol. Cancer Res. 12, 1069–1080. doi:10.1158/1541-7786.MCR-13-0672

PubMed Abstract | CrossRef Full Text | Google Scholar

Schüttelkopf, A. W., and Van Aalten, D. M. (2004). PRODRG: a Tool for High-Throughput Crystallography of Protein-Ligand Complexes. Acta Crystallogr. D. Biol. Crystallogr. 60, 1355–1363. doi:10.1107/S0907444904011679

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, S. A. M., Loccisano, A. E., Firestine, S. M., and Evanseck, J. D. (2006). Chapter 13 Principal Components Analysis: A Review of its Application on Molecular Dynamics Data. Annu. Rep. Comput. Chem. 2, 233–261. doi:10.1016/s1574-1400(06)02013-5

CrossRef Full Text | Google Scholar

Sterling, T., and Irwin, J. J. (2015). ZINC 15--Ligand Discovery for Everyone. J. Chem. Inf. Model. 55, 2324–2337. doi:10.1021/acs.jcim.5b00559

PubMed Abstract | CrossRef Full Text | Google Scholar

Trott, O., and Olson, A. J. (2010). AutoDock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J. Comput. Chem. 31, 455–461. doi:10.1002/jcc.21334

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A. E., and Berendsen, H. J. (2005). GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 26, 1701–1718. doi:10.1002/jcc.20291

PubMed Abstract | CrossRef Full Text | Google Scholar

Venugopal, P. P., M, S., and Chakraborty, D. (2021). Theoretical Insights into Molecular Mechanism and Energy Criteria of PARP‐2 Enzyme Inhibition by Benzimidazole Analogues. Proteins 89, 988–1004. doi:10.1002/prot.26077

CrossRef Full Text | Google Scholar

Wang, H., Nakamura, H., and Fukuda, I. (2016). A Critical Appraisal of the Zero-Multipole Method: Structural, Thermodynamic, Dielectric, and Dynamical Properties of a Water System. J. Chem. Phys. 144, 114503. doi:10.1063/1.4943956

PubMed Abstract | CrossRef Full Text | Google Scholar