Novel PARP-1 Inhibitor Scaffolds Disclosed by a Dynamic Structure-Based Pharmacophore Approach.

PARP-1 inhibition has been studied over the last decades for the treatment of various diseases. Despite the fact that several molecules act as PARP-1 inhibitors, a reduced number of compounds are used in clinical practice. To identify new compounds with a discriminatory PARP-1 inhibitory function, explicit-solvent molecular dynamics simulations using different inhibitors bound to the PARP-1 catalytic domain were performed. The representative structures obtained were used to generate structure-based pharmacophores, taking into account the dynamic features of receptor-ligand interactions. Thereafter, a virtual screening of compound databases using the pharmacophore models obtained was performed and the hits retrieved were subjected to molecular docking-based scoring. The drug-like molecules featuring the best ranking were evaluated for their PARP-1 inhibitory activity and IC50 values were calculated for the top scoring docked compounds. Altogether, three new PARP-1 inhibitor chemotypes were identified.

Introduction

Poly(ADP-ribose) polymerases (PARPs) comprise a group of enzymes that share the ability to catalyze the attachment of ADP-ribose moieties to specific acceptor proteins and transcription factors, using nicotine adenine dinucleotide (NAD+) as a substrate [1].

PARP-1 is the best characterized isoform among the PARP family members and is responsible for 85%-90% of poly(ADP-ribosylation) activity [2]. It plays an active role in several biological processes, including inflammation, hypoxic response, transcriptional regulation, maintenance of chromosome stability, DNA repair, and cell death [2-6]. The participation of PARP-1 in DNA repair granted it the designation of guardian angel of DNA [7]. This nuclear enzyme recognizes and binds to DNA strand-breaks via an N-terminal region, which promotes a conformational change in the C-terminal catalytic domain. As a result, this domain becomes activated, exposing the activation site to NAD+ and leading to the poly(ADP-ribosylation) of many targets, including histones and PARP-1 itself [3, 8].

The development of PARP-1 inhibitors as a therapy for several pathologies has been pursued, with special relevance in cancer and ischemic diseases [1]. The by-product of NAD+ cleavage, nicotinamide, has been used as the structural basis for the discovery of PARP-1 inhibitors. A large number of nicotinamide/benzamide derivatives have been studied, and some compounds have entered clinical trials as chemopotentiators in combination with anticancer drugs, as well as stand-alone agents in tumors with BRCA 1/2 mutations, taking advantage of synthetic lethality [8-11]. The drug candidate olaparib (LynparzaTM) was recently approved as the first PARP1/2 inhibitor to treat advanced ovarian cancer in women with defects in the BRCA1/2 genes, who were previously treated with three or more chemotherapeutic lines [12]. Nevertheless, a polypharmacological profile has been assigned to PARP-1 drug candidates. The inhibition of other PARP isoforms, or even the interaction with other inter-family targets, was noted for several inhibitors in clinical trials [1, 13]. Moreover, olaparib was reported to act as a substrate of the p-glycoprotein efflux pump, one of the mechanisms that are associated with resistance to PARP inhibitors [8, 14]. Clearly, more in-depth studies of the determinants of the PARP-1 recognition features are needed to develop novel and more selective PARP-1 inhibitors.

Computational methods have emerged as an important tool in drug discovery, as they disclose key features in the ligand-receptor binding interactions and allow the screening of large compound libraries, thus saving time and resources [15]. Moreover, molecular dynamics (MD) simulations have become an important method to solve one of the biggest challenges in drug discovery, i.e., the use of a single crystal structure of a protein to predict the putative ligand-binding site, not considering the target plasticity that is involved in ligand binding [16]. Different studies have combined MD with pharmacophore modelling, taking advantage of receptor flexibility to build structured-based pharmacophore models. In general, a wide array of drug discovery examples based on this approach have shown that they provide a better prediction of truly active compounds compared with inactive ones and are able to find potential leads for different targets under investigation [17-22].

In this work, a dynamic structure-based pharmacophore methodology was pursued to identify new scaffolds with PARP-1 inhibitory activity. A virtual screening of the available compounds databases was performed using the pharmacophore models generated, and the top scoring compounds identified by molecular docking studies were validated through an in vitro PARP-1 inhibition assay.

Materials and Methods

MD simulations

Four inhibitors that bound to the PARP-1 catalytic domain were retrieved from the Protein Data Bank (PDB codes: 2RCW, 3GN7, 3GJW, 3L3L). Crystal structures were processed using the Protein Preparation Wizard tool in Maestro Suite (Release 2013-1-9.4, Schrödinger, LLC, New York, NY, 2013). Water molecules were removed and atom types were assigned.

For each ligand-bound system, MD simulations in explicit water were performed using the Amber package, v12. Amber FF99SB [23] and Generalized Amber Force Field (GAFF) [24] were assigned to the protein and ligands, respectively. Systems were solvated with TIP3P water molecules [25] in a truncated octahedral box, counter ions were added to neutralize the system net charge, and the periodic boundary conditions were applied. The final systems were composed of ~ 33400 atoms.

After minimizations, systems were submitted to an equilibration phase for 1 ns in NVT conditions, in which protein and ligand atoms were position restrained with a constant force of 10 kcal/mol, to allow relaxation of the solvent molecules. A final production phase of 20 ns was performed and trajectory snapshots were saved at every 10 ps, for each system. The Langevin temperature equilibration scheme was used to keep the temperature constant (300 K), and a constant pressure periodic boundary was applied (1 atm). Electrostatic and Lennard-Jones forces were assessed using the Particle Mesh Ewald summation method [26] and a cut-off of 10 Å, respectively. The SHAKE algorithm was applied to constrain bonds that involved hydrogen atoms.

GROMACS [27] was used to perform the trajectory analysis. For each system, a conformational cluster analysis was carried out using a cut-off of 0.06–0.07 nm RMSD (root mean square deviation) between the backbone superposition of different structures. All snapshots saved from each MD trajectory were extracted and used to perform cluster analysis. To characterize the dynamics features of active site-ligand interactions in the PARP-1 catalytic domain, only the residues that were set to 5 Å around the inhibitor were taken into account in the cluster analysis.

Structure-based pharmacophore modelling and validation

Four different pharmacophore models were built based on the protein-ligand interactions observed after MD simulations. For each ligand-bound system, the clusters that represented more than 80% of the protein structural variability for each simulated system were selected to generate structure-based pharmacophores, using The Receptor-Ligand Pharmacophore Generation protocol of Accelrys Discovery Studio v3.5 (DS), Accelrys, San Diego, USA. This protocol uses receptor-ligand interactions to create selective pharmacophore models. Hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic (HY), negative ionizable (NI), positive ionizable (PI), and ring aromatic (RA) features, as well as the excluded volume spheres set to 5 Å around the inhibitor, were considered in the generation of the pharmacophore models. The hypotheses created were validated by a set of known PARP-1 ligands and decoys obtained from Database Useful Decoys Enhanced (DUD_E)—http://dude.docking.org/, using the validation option incorporated in the protocol. For each cluster, the hypotheses were ranked based on specificity and sensitivity, and the one that presented the best accuracy was chosen.

The best hypotheses that were retained for each cluster of a specified complex were superimposed, and the average coordinate point for each feature, including the excluded volumes spheres, was determined.

Four final pharmacophore models were created, one for each complex. As a final validation, the pharmacophore models were screened against the PARP-1 actives and decoys, to evaluate how well they discriminate active molecules from inactive ones. Moreover, the presence of chemical features that were essential for the interaction with key residues in the PARP-1 catalytic domain was taken into account in the validation of the pharmacophore models.

Database preparation and pharmacophore-based virtual screening

The National Cancer Institute (NCI)– https://cactus.nci.nih.gov/download/nci/ and DrugBank (http://www.drugbank.ca/) databases were downloaded. PARP-1 ligands and decoys were downloaded from the DUD_E database. Seven hundred and forty-two actives (affinity ≤1μM) and 30403 decoys (affinity ≥30 μM) were divided and converted into two databases, DUD_PARP1_ligands and DUD_PARP1_decoys, respectively. The “FAST” conformational analysis model of the catDB program was used to build the four databases, and a maximum of 255 conformations were generated for each molecule.

The four pharmacophore models obtained were used to screen the NCI and DrugBank databases using the “fast flexible database search” settings of Catalyst, to search for novel structural scaffolds with an ability to inhibit PARP-1.

The retrieved hits were subjected to different drug-like filters. Lipinski’s rule of five [28] and the modified Veber rule [29] (not more than 7 routable bonds) were applied. A maximum polar surface area was set to 140.

Docking Studies

The docking studies were performed using Glide (version 5.8). Standard precision (SP) and extra precision (XP) modes were applied, using the OPLS-AA force field [30].

The protein retrieved from the crystal structure of A620223 binding to PARP-1 (PDB code: 2RCW) was used to define the binding site. The Preparation Wizard tool was applied and all water molecules were removed from the crystal. A 15×15×15 Å receptor grid centered on the co-crystalized ligand was generated.

The final selected hits, as well as a set of know PARP-1 inhibitors (downloaded from BindingDB database (http://www.bindingdb.org)), were prepared using the LigPrep module (Schrödinger, LLC, New York, NY, 2013). The pH was set to 7.4 and a maximum of 5 stereoisomers per ligand were generated. The lowest energy ring conformation was kept for each stereoisomer.

An initial docking was performed using the SP-mode and 25 poses were kept for each molecule. A cut-off based on the docking score of reference PARP-1 inhibitors was used, and ligands with the highest score were subjected to XP docking.

PARP-1 enzyme assay

PARP-1 inhibition was evaluated using the HT Universal Colorimetric PARP Assay kit (Catalog #4677-096-K; Trevigen, Gaithersburg, MD, USA), in line with the instructions provided by the manufacturer. The assay evaluates the incorporation of biotinylated poly(ADP-ribose) onto histones proteins in a 96-well plate. Briefly, 10 μL of the test compounds were mixed with 15 μL of PARP-1 enzyme (0.5 U) into rehydrated histone-coated wells for 10 min at room temperature. Subsequently, 25 μL of PARP cocktail containing biotinylated NAD, activated DNA, and PARP buffer were added, and the solutions were incubated again for 60 min. After washing the wells, the detection reaction was performed according to the manufacturer’s protocol and absorbance was recorded at 450nm in a synergy HT plate reader. Stock solutions of the test compounds were prepared in dimethyl sulfoxide (DMSO) and serially diluted to the required concentrations with 1× PARP buffer. To assess the effect of the vehicle on enzyme activity, parallel experiments were performed by substituting the test compound with an equivalent volume of DMSO. IC50 values for the most promising hits were determined by plotting the inhibition data of each compound at different concentrations against the log of the concentration of the inhibitor, using the GraphPad Prism software, version 5. At least six different concentrations of the test compounds were used. A minimum of three independent assays were performed for each sample, and the results are displayed as mean ± standard error of the mean (SEM).

Nuclear Magnetic Resonance (NMR) studies

1D and 2D NMR structure elucidation of the NSC86342, NSC121848, and NSC131753 compounds was obtained using a Brucker Digital NMR-Avance 400 spectrometer, with CD3OD as the internal standard.

NSC131753 MD simulations

MD simulations were performed using (R)-NSC131753 and (S)-NSC131753 complexed with the PARP-1 catalytic domain, using top XP Glide poses as input structures. The MD simulations were performed using the protocol described above for the four complexes taken from PDB, with equilibration and production phases of 50 ps and 100 ns, respectively. Three replicas (100 ns) were run for each system with different initial velocities, to increase sampling.

MD trajectory analysis was performed using the GROMACS package.

Results and Discussion

Structural and dynamic characterization of different complexes with known inhibitors

In this work, MD simulations with different known small-molecule inhibitors were carried out to characterize the dynamic features of active site-ligand interactions in the PARP-1 catalytic domain. In this context, the aim of MD simulations was not the full sampling of the events underlying complex formation or the exploration of ligand induced conformational changes, which can be considered as being absent, given the high global similarity of the starting crystal structures, with a maximum RMSD (as calculated on protein backbone atoms) of 0.62 Å (). Rather, a comparative analysis of the trajectories from the different complexes was used to identify the salient features of the dynamic adaptation of PARP-1 to diverse active site inhibitors. Our general goal was to characterize the cross-talk between the ligands and the protein and highlight the binding interactions that were consistently preserved in multiple configurations, in addition to the ones that were immediately evident from crystal structures. Those conserved binding interactions were then used to develop dynamic pharmacophore models aimed at expanding the chemical diversity space of PARP-1 inhibitors.

The analysis of the main clusters revealed that the key interactions that were present in the co-crystal structures of PARP-1 with four different inhibitors were conserved. Such interactions consisted of three stable hydrogen bonds: two between the amide backbone of Gly202 and the amide moiety of the inhibitors and one between the OH group of Ser243 and the carbonyl group of inhibitors, as presented in , as well as a π-π stacking interaction involving Tyr246 and the aromatic core of the ligands. Furthermore, MD simulations showed that the tyrosine residues present in the binding site were involved in different π-interactions. To define the importance of the tyrosine residues, the contacts between Tyr228, Tyr235, and Tyr246 and the ligands were monitored during 20 ns MD run (). Depending on the bound inhibitors, different tyrosine residues were engaged. Tyr228 appeared to be essential for the π-cation interaction with the protonated amine moiety of the 2RCW and 3L3L ligands, with occupancy of 99% and 100%, respectively, during MD trajectories. For the 3GN7 inhibitor, both Tyr228 and Tyr246 were implicated in this type of interactions, with occupancy of 65.45% and 47.5%, respectively. Moreover, Tyr235 and Phe236 were involved in amide-π interactions with the phenyl rings of the 2RCW, 3GN7, and 3GJW ligands. In addition, Ala237 and the alkyl side chain of Lys242 participated in hydrophobic interactions with the phenyl ring of each inhibitor along the MD trajectory for all ligand-bound systems analyzed (). The protonated amine group of each ligand also appeared to be important for the establishment of charge-charge interactions with some charged residues present in the binding site, such as Glu102, Asp105, and Asp109 (). This type of interaction was especially relevant for 3GJW. The 3GJW ligand was involved in charge-charge interactions with Asp105 (:OD1) for 84.65% of the MD run time. The main interactions for each complex along MD trajectories (the first cluster) are shown in

The plot distances involving the tyrosine residues and ligands along MD trajectories.

A) Distance between Tyr246 centroid π ring and the aromatic core of ligands. B) Distance between Tyr228 centroid π ring and the protonate amine moiety of each ligand. C) Distance between the centroid of amide group (formed by the CO of Tyr235 and N of Phe236) and phenyl ring of each ligand.

2D-Ligand interaction diagrams for each ligand complexed with PARP-1 catalytic domain along MD run.

Dashed lines represent interactions between binding site residues and bounded ligands. Green color pointed to hydrogen bond interactions; Orange indicates π-cation interactions; pink denotes π-π stacked; yellow pointed to amide-π interactions; light pink denotes hydrophobic interactions (alkyl and π-alkyl); brown indicates charge-charge interactions, and turquoise residues indicate van der waals interactions. A) 2RCW. B) 3L3L. C) 3GN7). D) 3GJW.

By exploring the diversity and the motion of the ligands, as well as the flexibility of the binding site residues, four structure/dynamics-based pharmacophores were generated based on the ligand-protein interactions that were monitored during the MD trajectories. The conserved hydrogen bonds, as well as the π-π stacking, the π-cation and charge-charge interactions observed to a higher extent during MD, were considered to generate structure-based pharmacophores.

For each complex, seven representative structures (matching more than 80% of the structural variability) were taken into account to build the pharmacophore models.

For all pharmacophore hypotheses, the Receptor-Ligand Pharmacophore Generation protocol in DS pointed key interactions between the PARP-1 catalytic domain and the ligand, and generated excluded volume spheres that were correlated with steric regions in the binding site that may not be engaged by the ligand substituent groups. In this context, the characterization of the mechanisms of the formation/disappearance of pockets around the ligands due to the immediate conformational response of the protein to known inhibitors can aptly indicate the positions at which the addition/modification of specific substituent groups may allow optimal extensions of binding interactions into previously uncharacterized regions.

The comparison of the excluded volumes spheres obtained based on the crystal and the representative structure after MD simulations for each ligand-bound system () revealed that the excluded volume spheres were generally pointed to the same residues, especially in the nicotinamide binding pocket, which comprises residues such as Gly202, Ser243 and Tyr246 (). Analyses of root mean square fluctuation (RMSF) for all four complexes () revealed that the regions that contained nicotinamide binding residues were quite stable. The highest fluctuations were observed in loop regions (60–67; 78–94; and 118–128).

Representation of the PARP-1 binding site with excluded volume spheres.

The excluded volume spheres, set to 5 Å around each inhibitor complexed with the binding site, were obtained from the crystal and the representative structure after MD simulations for 2RCW, 3L3L, 3GN7 and 3GJW complexes.

Pharmacophore model building

SB_Pharma1, which was based on the 2RCW complex, displayed five functional features, including an HBD and an HBA pointed to Gly202, an HY pointed towards Ala237 and Lys242 (alkyl side chain), an RA directed to Tyr246, and a PI pointed to Tyr228. Nineteen excluded volume spheres were identified, which represented an additional two spheres compared with the crystal structure. One of them was directed to Gln98 (side chain) and the other to Thr226 (backbone). SB_Pharma2, which was the pharmacophore model obtained from 3L3L, exhibited four features and 14 excluded volume spheres, pointed to the same residues of the crystal structure. One HBA and one HBD directed to Gly202, an HY pointed towards Ala237 and Lys242 (alkyl side chain), and a PI pointed to Tyr228 were observed. SB_Pharma3 and SB_Pharma4 were generated from 3GN7 and 3GJW, respectively. Both pharmacophore models exhibited three similar features: an HBA and an HBD pointed to Gly202, an HY center directed towards Ala237 and Lys242 (alkyl side chain) and an RA also directed to Tyr246. A PI center was pointed to Tyr228 or even Tyr 246 in SB_Pharma3, and to Asp105 in SB_Pharma4. Excluded volume spheres (15 and 16, respectively) were also identified. SB_Pharma3 displayed an extra excluded volume (compared with the 3GN7 crystal structure) pointing towards the Asp105 side chain that changed side chain orientation during the MD simulation (). Moreover, SB_Pharma4 also showed an additional excluded volume sphere, directed to Tyr228, which side chain exhibited considerable flexibility along the MD trajectory (Figs ). The excluded volume directed to Arg204, which was observed in the 3GJW crystal structure, was not set in the final pharmacophore model obtained after MD. As illustrated in , pointing to this active site residue did not appear to be essential for productive interaction. The final structure-based pharmacophore models (SB_Pharma1, SB_Pharma2, SB_Pharma3, and SB_Pharma4) obtained from the superposition of dominant conformations for each ligand-bound system, as well as those obtained from PDB crystal structures, are elucidated in .

Representation of structure-based pharmacophore models obtained from crystal structures and the dominant conformations after MD simulations.

Green color indicates hydrogen bond acceptor (HBA); magenta denotes to hydrogen bond donor (HBD); cyan shows hydrophobic center (HY); yellow indicates ring aromatic (RA); and red denotes to positive ionizable center (PI).

The analysis of four structure-based pharmacophore models supported the essential role of Gly202 as an HBD and HBA, as well as the presence of important hydrophobic residues, namely Ala237 and Lys242 (alkyl side chain). Moreover, Tyr228, Tyr235, and Tyr246 were shown to be important for the establishment of different types of π-interactions. The aromatic feature directed to Tyr246, for instance, was necessary to mimic the relevant role of stacking interactions in driving effective binding to the PARP-1 catalytic domain. Finally, the presence of charge-charge interactions mediated by charged residues, as exemplified by Asp105, may be important for the identification of additional interactions that increase the binding affinity between the ligand and the protein.

Pharmacophore-based virtual screening and validation

The four pharmacophore models were validated against DUD_PARP1_ligands and DUD_PARP1_decoys, which were generated by Catalyst, and both sensitivity and specificity were calculated (). Sensitivity was related to the fraction of PARP-1 binders that correctly fit the pharmacophore models. Specificity was related to the fraction of molecules that did not fully fit the pharmacophore hypotheses and were identified as decoys. The comparison of the values obtained for the four structure-based pharmacophore models showed that SB_Pharma1 and SB_Pharma4 displayed a better accuracy compared with the already good one characterizing all pharmacophore models, in general. On such bases, all four hypotheses generated were used to screen the NCI and DrugBank databases. However, to increase the ability to distinguish between active and inactive molecules, only hits with fit values above 2.0 and those that were retrieved by more than one pharmacophore model (in which at least one of them displayed the best accuracy (SB_Pharma1 or SB_Pharma4)), were retained for further docking studies. Overall, 915 and 175 hits were obtained from the NCI and DrugBank databases, respectively.

Importantly, inspection of the retrieved hits identified 3 known PARP-1 inhibitors among the 175 molecules that were obtained from the screening of the DrugBank database: DB0372 (FR257517) [31], DB07787 (FR255595) [32], and D08348 (PJ34) [1]. It is worth noting here that these ligands were not part of the initial training set of ligands that was used to start MD simulations and pharmacophore design. The presence of these inhibitors, which have a higher potency regarding the inhibition of PARP-1 activity, constituted a first important validation of the capacity of our pharmacophore models to recapitulate the chemical and stereoelectronic determinants that underlie the activity of drug molecules.

Docking studies

The overall 1090 retrieved hits (from the NCI and DrugBank databases) were docked at the PARP-1 binding site using Glide SP-mode. To validate and optimize the docking parameters, A620223 co-crystalized with the PARP-1 catalytic domain (PDB code: 2RCW), as well as 14 reference PARP-1 inhibitors, were re-docked. The SP docking results showed that the binding pose of A620223 in the crystal could be optimally reproduced, with an RMSD of 0.64 Å (). Furthermore, the top poses of hits retrieved from the NCI and DrugBank databases were inspected, and a docking score cut-off of -7 was applied, based on the docking score range of the PARP-1 inhibitors that were docked ([-7.5;-11.6]). The remaining compounds were subjected to a second docking run, using Glide XP-mode. To select promising hits, a visual inspection of the compounds was performed. The interaction with key residues, such as Gly202 and Tyr246, as well as the structural diversity between the molecules, was taken into account when choosing potential hits. A total of 60 compounds were chosen for further evaluation.

PARP-1 inhibition and structure-activity relationship

The HT Universal Colorimetric PARP Assay Kit was used to screen and to determine the IC50 values of the promising hits obtained. Only 39 compounds among the 60 chosen above were effectively tested, because of commercial availability or solubility problems. After an initial screening at a concentration of 100 μM, seven compounds displayed a PARP-1 inhibition activity >90%. A new screening at 10 μM was performed. The IC50 was determined for the most promising hits (and ).

PARP-1 inhibitory activity and docking score data for the most promising hits.

) PARP-1 inhibition was determined using HT Universal Colorimetric PARP Assay Kit (Cat #4677-096-k). *Docking score values of both possible NSC86342 cis diastereomers.

Dose-response curves of the three most promising hits.

Each data point represents the mean ± SEM of at the least three independent experiments.

Among the promising molecules, three of them exhibited one or more chiral centers. To determine which isomer was acquired from NCI, the NSC86342, NSC131753, and NSC121848 compounds were characterized by NMR. The 1H NMR and 13C NMR spectra showed that only an isomeric form was present for each sample (). Moreover, the NOESY spectrum allowed the characterization of the enantiomeric form obtained for NSC121848. The H9 proton observed in the NOESY spectrum was correlated with both methyl groups at N1, which indicates that the (R)-enantiomer was present (). Similarly, in the NOESY spectrum of NSC86342, a correlation between H10 and the methyl group at N1was observed. This demonstrates that these two groups have the same orientation, which reveals that NSC86342 is a cis diastereomer (). As can be observed, there was a huge structural variability between the most promising compounds, with NSC131753 showing the highest PARP-1 inhibitory activity (IC50 = 0.24 μM). Moreover, PARP-1 inhibition was well correlated with the XP docking scores of three among the top four most promising hits (). Despite the fact that NSC86342 showed the lowest docking score, it displayed π-cation and π-π interactions with key tyrosine residues (Tyr235 and Tyr246), which have been described as being essential for the binding of PARP-1 inhibitors to the catalytic domain. Moreover, some poses revealed π-charge interactions with Glu102, Glu327, and Tyr228. Taken together, these findings may explain the stability of this compound at the binding site and its high PARP-1 inhibitory activity (and .

The binding mode of the most promising PARP-1 inhibitors at the PARP-1 catalytic domain.

The molecular interactions of the top scored poses were displayed. A) NSC86342. B) NSC121848. C) (R)-NSC131753. D) (S)-NSC131753.

As expected, all promising hits were involved in interactions with conserved binding residues, such as Gly202, Tyr246, or even Tyr235. The presence of donor-acceptor aromatic systems appeared to be essential for PARP-1 inhibition, which is in line with the hydrophilic environment that surrounded the NAD+ binding pocket, with a remarkable presence of aromatic residues (). Consequently, it is easy to understand why the best PARP-1 inhibition activities were displayed by molecules with an aromatic polycyclic skeleton with several HBAs or HBDs, such as NSC131753 and NSC121848 (. The latter established hydrogen bonds with a vast number of residues, such as Phe236, Gly227, and Glu327, in addition to the π-cation and π-π interactions with the key Tyr235 and Tyr246 residues. Although the NMR analysis did not determine which NSC131753 enantiomer was evaluated, both (R) and (S) enantiomeric forms may be involved in different types of interactions with the catalytic domain, in spite of the differences in docking scores. Hydrogen bond interactions involving Gly233 and Phe236 were stablished with the (R)-enantiomer (), while Gly227 and Met229 were implicated in this type of interactions with the (S)-enantiomer (). Moreover, Glu327 played an important role in the interaction profile of both enantiomeric forms, by establishing π-cation interactions with (R)-NSC131753 and H-bond with (S)-NSC131753. Further insights into the binding mode of the two enantiomeric forms will be discussed in the description of the NSC131753 MD simulations analysis.

It is worth mentioning that anthraquinone derivatives, of which NSC102534 is an example, have been recently reported as being PARP-1 inhibitors [33]. The polycyclic aromatic core of these compounds was crucial for the interaction with the binding site ().

In addition, it is important to note that the hits with the most promising PARP-1 inhibitory activity, NSC131753, NSC86342, and NSC121848 consist, to the best of our knowledge, in new PARP-1 inhibitor skeletons. These compounds shared the ability to interact not only with conserved nicotinamide-binding pocket residues, such as Gly202, Tyr235, and Tyr246, but also with some residues located on a donor-site loop (Gly215-Gly233), such as Gly227 (NSC121848), Tyr228 (NSC86342), and Met229 ((S)-NSC131753); this, could explain the binding mode stability and the relevant PARP-1 inhibition values observed.

The most promising hit, NSC131753, contains a chiral center. Despite the performance of 1D and 2D NMR studies, it was not possible to identify the enantiomeric form evaluated against the PARP Assay kit. To determine which enantiomer is more stable at the binding site and to attest the interaction profile obtained from the docking studies, 100 ns long MD simulations were carried out for (R)- and (S)-NSC131753 complexed with the catalytic domain. The RMSD was lower for (S)-NSC131753 (around 0.05 nm compared with 0.09 nm for (R)-NSC131753), even though the two enantiomers revealed being quite stable (). However, the average RMSD calculated on the Cα atoms was lower and more stable for the (R)-enantiomer complex, along the three 100 ns MD replicas ()). Conversely, a similar RMSF distribution () was observed for both enantiomeric forms, with the highest fluctuation observed in two loop regions of the catalytic domain (78–94; 118–128). Nevertheless, a highlighted mobility was observed from residues 317 to 322 in the RMSF plot of (S)-NSC131753, which was not observed for (R)-NSC131753. This may indicate a conformational change induced by the (S)-enantiomer. It is also worth noting that the D-loop residues (215–233) presented higher flexibility in the case of (R)-NSC131753 (at least 0.02 nm) compared with (S)-NSC131753 (around 0.015 nm). This difference may be due to the interaction of (S)-NSC131753 with Tyr228 (through π-π and π-cation interactions) and Met229 (hydrogen bond). In fact, the analysis of the interaction profile during MD showed that both enantiomers established an important number of interactions. Moreover, the main interactions proposed by the docking studies were maintained for both enantiomers, with high occupancy. Tyr235 and Tyr246 were involved in π-π interactions with both enantiomers, with occupancies above 65%. Glu327 was essential for the establishment of hydrogen bonding with both enantiomeric forms, and for charge-charge interactions with (R)-NSC131753, as demonstrated by the docking studies described above. MD trajectories analysis also revealed that Glu102 established a hydrogen bond interaction with the (R)-enantiomer for 41.5% of the MD run time, and with the (S)-enantiomer for 35.01%. A similar type of interaction was identified between ASP105 and (S)-NSC131753, with an occupancy of 38.55% along 100 ns MD simulations.

Conformational statistics obtained for (R)- and (S)-NSC131753, along 100 ns MD run.

(A, D) RMSDs of ligands. (B, E) RMSDs of complexes. (C, F) RMSFs of Cα atoms of catalytic domain.

In summary, the data showed that both enantiomers were able to interact with the catalytic domain with relative stability, via different types of interactions with binding site residues, some of which were revealed only during MD simulations. An example of this is the interaction between (S)-NSC131753 and Tyr228, which was important to stabilize the D-loop and may explain the differences in docking scores observed between (R)- and (S)-NSC131753.

Conclusions

A dynamic structure-based pharmacophore strategy was used to identify novel PARP-1 inhibitors. The pharmacophore models based on the interactions between the PARP-1 catalytic domain and four different inhibitors during MD simulations provided new insights in the ligand binding mode, taking into account the flexibility of both the enzyme and the ligand. Subsequently, the validated pharmacophore models were screened against two virtual compound libraries, to retrieve hits with novel chemical scaffolds. After molecular docking studies using Glide, the top scored drug-like molecules were tested against the PARP kit assay to determine PARP-1 inhibitory activity. Structurally diverse hits with important PARP-1 inhibitory activity were found. Moreover, the dynamic structure-based pharmacophore approach applied here led to the identification of three new PARP-1 inhibitor candidates with skeletons that had not been reported previously: NSC86342, NSC131753, and NSC121848. These candidates will be useful for guiding the further development of novel and more potent PARP-1 inhibitors.

Binding pocket residues of 2RCW.

Nicotinamide (A) and adenine-ribose (B) binding residues are displayed in red and blue, respectively. Violet was used to show D-loop residues (C).

(PDF)

Click here for additional data file.

Asp105 side chain orientation in the crystal structure and after MD simulations (representative structure).

(PDF)

Click here for additional data file.

1H-NMR and 13C-NMR data for NSC86342, NSC121848 and NSC131753.

(PDF)

Click here for additional data file.

NOESY spectrum of NSC121848.

(PDF)

Click here for additional data file.

NOESY spectrum of NSC86342.

(PDF)

Click here for additional data file.

The binding mode of NSC153161 (A), NSC102534 (B), NSC65378 (C) and NSC11907 (D) at the PARP-1 catalytic domain. The molecular interactions of the top scored poses were displayed.

(PDF)

Click here for additional data file.

Table 1

Hydrogen bonds (A), hydrophobic interactions, namely alkyl and π-alkyl interactions (B), and charge-charge interactions (C) with greater occupancy during MD trajectories for 2RCW, 3L3L, 3GN7 and 3GJW complexes.

A
Hydrogen bonds
PDB ID	Donnor	Acceptor	%Occupancy
2RCW	Gly202:N	AAI:O1	40.05
	Ser243:OG	AAI:O1	64.80
	AAI:N3	Gly202:O	73.55
3L3L	Gly202:N	L3L:O2	54.75
	Ser243:OG	L3L:N3	38.85
	L3L:N3	Gly202:O	64.35
	L3L:N1	Gly227:O	23.80
3GN7	Gly202:N	3GN:O12	38.80
	Ser243:OG	3GN:O12	44.50
	3GN:N1	Gly202:O	61.70
3GJW	Gly202:N	GJW:O1	27.15
	Ser243:OG	GJW:O1	62.65
	GJW:N3	Gly202:O	74.80
B
Alkyl and π-alkyl interactions
PDB ID	Residue	%Occupancy
2RCW	Ala237	68.00
	Lys242	97.70
3L3L	Ala237	48.00
	Lys242	84.70
3GN7	Ala237	52.20
	Lys242	75.70
3GJW	Ala237	66.65
	Lys242	64.55
C
Charge-charge Interactions
PDB ID	Residue		%Occupancy
2RCW	Glu102	OE1	29.40
		OE2	32.50
	Asp105	OD1	30.30
		OD2	26.35
3L3L	Glu102	OE1	33.25
		OE2	37.60
3GN7	Glu102	OE1	29.15
		OE2	32.50
3GJW	Asp105	OD1	84.65
		OD2	19.65
	Asp109	OD1	30.95
		OD2	23.05

Table 2

Statistical data of structure-based pharmacophore models.

Pharmacophore Model	TA	TI	TP	TN	FP	FN	Se	Sp	Acc
SB_Pharma1	742	30403	346	27926	2477	396	0.466	0.918	0.908
SB_Pharma2	742	30403	364	26037	4366	378	0.491	0.856	0.848
SB_Pharma3	742	30403	275	25873	4530	467	0.371	0.851	0.840
SB_Pharma4	742	30403	225	28073	2330	517	0.303	0.923	0.908

TA: Total number of actives; TI: Total number of inactives; TP: True positives; TN: True negatives; FP: False positives; FN: False negatives; Se: Sensitivity; Sp: Specificity; Acc: Accuracy.