1. Introduction

Polyadenosine diphosphate-ribose polymerase-1 (PARP1) was found by Pierre Chambon about six decades ago (Singh et al., 2014). This enzyme was extracted and purified from the liver nuclei of chicken (Jeggo, 1998). It uses NAD+ and cleaves/transfers the ADP-ribose to acceptor proteins to produce long poly ADP-ribose (PAR) chains (Ekblad, Camaioni, Schuler, & Macchiarulo, 2013; Hottiger, Hassa, Luscher, Schuler, & KochNolte, 2010; Steffen, Brody, Armen, & Pascal, 2013). PARP1 is encoded by the ADP ribosyl transferase-1 gene (MegninChanet, Bollet, & Hall, 2010) and plays its part by stopping the death of cells. PARP1 enzyme plays important roles in various biological functions such as in cell development, inflammation, during stress responses, cancer, etc. (Singhetal., 2014). Its key role, target protein parylation (PAR formation) is commonly involved in many diverse biological functions and processes such as transcriptional regulation, genomic stability preservation and cell death (Liscio, Camaioni, Carotti, Pellicciari, & Macchiarulo, 2013; Luo & Kraus, 2012; Schreiber, Dantzer, Ame, & de Murcia, 2006).

In the modern drug design approaches, molecular targeting must be focused. The PARP1 is one of the targets used in the design of anticancer drugs. It is active in many cellular processes, from the repair of DNA to the destruction of cells (Kraus & Hottiger, 2013). PARP1 binds the damage site when DNA is damaged and catalyzes PAR synthesis recruiting host DNA repair proteins. PARP1 AutoPARylation causes electrostatic destabilization and dissociation from the site of DNA damage which enables protein repair to be found to the lesion (Ogata, Ueda, Kawaichi, & Hayaishi, 1981; Satoh & Lindahl, 1992). Olaparib (Lynparza) and Rucaparib (Rubraca) which are FDAapproved PARP inhibitors for the treatment of BRCA-mutated ovarian cancer (Lord & Ashworth, 2017) and Niraparib (Zejula) for the treatment of recurrent ovarian cancer with or without BRCA mutation (Scott, 2017) are significant milestones for the discovery of novel compounds targeting this enzyme.

The goal in anticancer drug design targeting this enzyme is the inhibition of excessive activity of PARP1 function. PARP-1 inhibition enhances the radiation and chemotherapeutics influence (Michels, Vitale, Saparbaev, Castedo, & Kroemer, 2014; Plummer, 2014). PARP1 inhibitors, particularly in BRCA-deficient cells, have anticancer activity (Curtin & Szabo, 2013). Since PARP1 is a main enzyme which controls such carcinogenic modifications in the body, it is considered as an important molecular target for antitumor agents hence PARP1 inhibitors are used in both academy and in industry as crucial anticancer hits promising anticancer drugs. In addition, PARP-1 inhibitors may be of significant benefit to certain disorders such as diabetes mellitus, inflammation, aging, cardiovascular disorders and neurodegenerative disorders (i.e. Parkinson’s and Alzheimer’s diseases) (Burkle, Diefenbach, Brabeck, & Beneke, 2005; Graziani, Battaini, & Zhang, 2005; Sodhi, Singh, & Jaggi, 2010; Szabo, 2005).

In recent years, the fact that new lead molecules identified from small molecule data banks which give promising results in different biological problems has increased the interest in this field (Ekhteiari Salmas, Senturk, Yurtsever, & Durdagi, 2016; Ekhteiari Salmas, Unlu, Yurtsever, Noskov, & Durdagi, 2016; Ekhteiari Salmas et al., 2017; Kanan et al., 2019; Mirza, Salmas, Fatmi, & Durdagi, 2016; Sahin, Zengin Kurt, Sonmez, & Durdagi, 2019). Piperazine has received considerable interest in the development of many effective biologically active compounds, including anti-proliferative agents (Shaquiquzzaman et al., 2015; Zeng et al., 2009). FDA approved piperazine-derivative anticancer drug imatinib is currently on the market and it is used to treat many cancers, including myelogenous leukemia (O’Brien et al., 2003). FDA has recently approved another piperazine including anticancer drug Palbociclib for the treatment of metastatic breast cancer (Morikawa & Henry, 2015). In addition, for the treatment of chronic myeloid leukemia, Ponatinib that includes piperazine ring has been also produced (Pavlovsky, Chan, Talati, & Pinilla-Ibarz, 2019). Moreover, one of the currently used FDA approved PARP1 inhibitors Olaparib (Johannes et al., 2015) also includes piperazine ring.

Text mining helps to find molecules of interest by screening a large database with keywords quickly (Krallinger & Valencia, 2005). Thus, in the current study, a text mining study was considered for identifying the piperazine derivatives from a large database. For this aim, 211960 molecules obtained from SpecsSC database have been prepared using MarvinSketch program (MarvinSketch 18.30.0, Chemaxon) in IUPAC text file format. To find the compounds that include “piperazine” phrase, using an in-house Python-based text mining script, this text file was screened and 4674 molecules were obtained. Filtered piperazine-based compounds were then used in our rigorous virtual screening algorithm to find new hits against PARP1.

2. Materials and methods
2.1. Ligand preparation

The 2D structure of selected 4674 ligands were converted to 3D structure using the Maestro molecular modeling suite’s LigPrep module (LigPrep, Schrodinger, LLC, New York, NY, 2015) with the OPLS-2005 forcefield (Bankset al., 2005). The correct protonation states of the ligand in physiological environments is a problem to be noticed. The Epik module (Shelley et al., 2007) of Maestro molecular modeling suite has been used at neutral pH for the correct estimation of protonation states of the filtered compounds from text-mining studies. It must be noted that potential steroisomers and tautomers of identified compounds were also generated.

2.2. Protein preparation

The protein data bank (PDB) (http://www.rcsb.org) provides 3D structures of different nucleic acids, proteins, protein fragments and ligand complexes of proteins. Since 5DS3 PDBcoded structure includes one of the FDA-approved drug (Olaparib, which includes piperazine group), we used this protein structure as target. However, since this target includes only residues between 788 and 1012, we decided to use another PARP1 target (PDB ID: 4HHY) (Souers et al., 2013), which includes residues between 660 and 1011. Thus, we superimposed these two target structures and aligned the Olaparib and crystallized water molecules around the Olaparib at the 5DS3 PDB-coded structure to the binding pocket of 4HHY. Then, Olaparib-bound 4HHY PDB-coded structure was used as target file. The entire structure was pre-processed at the Protein Preparation module of Maestro with the default settings. Prime module of Maestro fixes the missing side chains, atoms and loops on the backbone (Jacobson et al., 2004). The impact of water molecules at the catalytic domain must be taken into account, thus water molecules within a sphere of 5.0 Å radius from center of mass of ligand were used in this study. The PROPKA and OPLS-2005 force fields were used for protonation states, structural optimization and minimization (Protein Preparation, Version 2.5, Schrodinger, LLC, New York, 2015), respectively.

2.3. Molecular docking simulations

Molecular docking studies examine associations and candidates in protein-protein and ligand-protein complexes by a predicted affinity scoring method (,S()led,z() & Caflisch, 2018).Docking processes predict ligands to bind in the binding pocket of target protein with the most appropriate conformation using different scoring and conformational search algorithms. The Glide was used to perform the molecular docking between the selected ligands and PARP1 target to achieve the Salivary microbiome poses with high docking scores. In Glide, (Raha & Merz, 2004), the conformations obtained during the docking are scored using the Glide scoring function. A total of 10 bioactive conformers were produced for each ligand and residues were positioned in the binding pocket inside 10Å of co-crystallized ligand. All ligands were docked in the PARP1 binding pocket using the grid-based docking program Glide Standard Precision (SP) of the Maestro molecular modeling method (Friesner et al., 2004; Jacobson et al., 2004) to produce the top-docking poses.

2.4. Molecular dynamics (MD) simulations

To decide the correct binding matches, studied complexes at the docking studies require MD simulations (Hou, Wang, Li, & Wang, 2011). For MD simulations, the top-docking complexes were used with the Desmond v 4.9. MD simulations are conducted to determine structural and dynamics changes throughout the protein-ligand interactions to investigate the conformational stability of identified hit compounds at the PARP1 binding site (Desmond, Version 4.9). The complex structures are solvated with the orthorhombic simple point charge (SPC) water models (Berendsen, Postma, van Gunsteren, & Hermans, 1981). Counter ions were used to neutralize the systems (0.15M NaCI solution). The system was set as Lennard-Jones interactions cut-off of 10Å on periodic boundary conditions (Friesner et al., 2004). In the integration steps, 2.0 fs timestep was used. Nose-Hoover thermostat (Hoover, 1985) and Martyna-Tobias Klein protocols (Ma & Tuckerman, 2010) were used to control the temperature and pressure of the systems at 310K and 1.01325bar, respectively. Before the production run, the system is minimized and equilibrated. 100-ns MD simulations for each system were then performed.

Figure 1. Applied virtual screening workflow.

2.5. Molecular mechanics/generalized born surface area (MM/GBSA)

MM/GBSA was used for protein-ligand complexes to measure the free energy of binding. Using the Maestro Prime module (Jacobson et al., 2004), MM/GBSA calculations were carried out for the studied systems. The trajectories of proteinligand complexes were used from the MD simulations of each selected hit compound at every 10 ps. (Desmond, Version 4.9) . VSGB solvation system (Li et al., 2011), which is a practical solvation parameterization, and OPLS-2005 forcefield (Banks et al., 2005) have been used for the flexibility of proteins. In the binding free energy predictions of selected hits by MM/GBSA calculations, 100 trajectories during the whole production MD simulations (100-ns) were used and average values were calculated.

3. Results and discussion

In the current study, we aimed to identify novel piperazinebased PARP1 inhibitors using text-mining and integrated molecular modeling approaches. The entire procedure used in the current study has been shown in Figure 1.Using the LigPrep module of the Schrodinger Maestro Molecular Modeling Suite, the two-dimensional molecular structure of obtained 4674 piperazine-based small molecules have been converted into energy-optimized 3D structures. The molecular docking technique was performed by Glide/SP to obtain the hit compounds that give high docking scores at the binding pocket of the PARP1 enzyme. The grid-based docking was initially conducted on the PARP1 binding site for 4674 piperazine derivatives obtained by text-mining. Based on docking scores within 4674 compounds, we selected top-10 compounds that have high docking scores as hit molecules and MD simulations were performed for these complexes. In order to compare the binding free energies, one of the known FDA approved drugs (Olaparib) was also considered in MD simulations.Complex structures are placed in standard orthorhombic boxes including the TIP3P water models and the box size buffer length range was set to 10Å in all three directions. Na+/Cl-ions were used for neutralization process, thus 0.15M NaCl solution was used to imitate the body environment. Simulations were carried out at constant temperature and pressure (NPT ensemble).100ns MD simulations were performed for each of 10 hit molecules at the binding pocket of the PARP1 as well as FDA approved PARP1 inhibitor Olaparib and the corresponding coordinates throughout the MD simulations were collected separately in the trajectory file. The obtained results of ligand docking and post-processing binding free energy analysis from MD simulations trajectories were summarized at Table 1.

Root-mean square deviations (RMSDs) of Ca atoms of studied systems were plotted at Figure S1 (Supplementary material). It can be seen that all of the studied systems relaxed and do not represent significant fluctuations at the RMSD values after 20-ns of MD simulations. Average RMSD values of all of the studied systems were less than 3.0Å.RMSDs of studied ligands were also investigated throughout the MD simulations. The LigFitLig RMSD plot is shown in Figure S2 (Supplementary material). This plot represents the RMSD of studied compound based on its initial conformation during the simulation. Hence, the RMSD value represents the internal fluctuations of nonhydrogen atoms of the ligand.LigFitProt RMSD plot shows RMSD of ligands when the protein-ligand complex is initially aligned as a reference state with the backbone atoms of the protein. Then, the RMSD of the heavy atoms (i.e., non-hydrogens) of the compounds were calculated (i.e., structural stabilities of the compounds at the binding site during the simulations). . The LigFitProt RMSD plot is depicted in Figure S3 (Supplementary material). Corresponding average RMSD values are represented in Table 2. Most of the analyzed compounds have slight rotational movements throughout the MD simulations in the binding site.

Root-mean square fluctuation (RMSF) values were also calculated during simulations to examine the impact of defined hits on the mobility of target protein’s backbone atoms.RMSF of polypeptide backbone atoms of each amino acid residue were produced in the complex analyzes to identify the target structure’s fluctuation regions. High RMSF values in the RMSF plot suggest highly mobile areas during MD simulations. The RMSF graph of 10 hit molecules and FDA approved PARP1 inhibitor Olaparib is shown in Figure 2. Peaks on this graph indicate the most fluctuating residues in the protein during the simulations. The tails (N-and C-terminals) are usually found to be more fluctuating than any other protein component. While secondary protein structures show more rigid profiles throughout the simulations, loop domains fluctuate more. . Figure 2 shows the effect of all ligands on the binding pocket as well as on the entire protein structure. Figures 3 and 4 display associations that arise in the chosen path over 15% of the simulation time (0.00– 100.00ns). For selected hit molecule-1388 which have similar MM/GBSA score with Olaparib, crucial interactions are constructed via hydrogen bonds with Gln759 and Met890 and π-π stacking this website interactions is formed with Tyr889. It must be noted that Gln759 was observed as one of the common crucial residues both in Olaparib and molecule-1388.Figures 5 and 6 highlight a timeline representation of interactions and contacts. The top panel indicates the total number of specific contacts between ligand and protein while the bottom panel represents which amino acid residues at the target protein interacts with the studied ligand during MD simulations.Many residues, according to the scale to the right of the map, allow more than one direct interaction with the ligand, which is defined by a darker orange color.Protein-ligand interactions graph of Mol-1388 shows that following residues have stable nonbonded chemical contacts throughout the simulations: Ala762, Gln759, Tyr889 and Met890. While protein-ligand interactions graph of Mol-1388 shows that Ala762, Gln759, Tyr889 and Met890 form stable nonbonded chemical interactions (Figure 5), Olaparib adopts corresponding interactions with following residues: Gly863, Arg878, Tyr896, Ser904 and Tyr907 (Figure 6).Moreover, together with Olaparib, we also investigated other three known PARP1 inhibitors with their co-crystallized PARP1 complexes using same MD simulations parameters with hit compounds. MM/GBSA results show that while Olaparib has slightly better average interaction energy score than the identified hit compounds, some of the selected hits (i.e., Mol-1130, Mol-1388, Mol-3915) have similar or better average interaction energy values than other three FDA approved drugs.

Figure 2. RMSF graph for selected 10 piperazine-derivative compounds and FDA approved drug (Olaparib) throughout the MD simulations. X-axis as the figure shows residue numbers at the target protein.

Figure 3. 2D Protein-ligand interactions graph for the selected hit compound Mol-1388.

Figure 4. 2D Protein-ligand interactions graph for the FDA-approved PARP1 inhibitor Olaparib.

Figure 5. Protein-ligand interactions graph of Mol-1388 throughout the simulations.

Figure 6. Protein-ligand interactions graph of Olaparib throughout the simulations.

In order to check the selectivity profiles of our hit compounds, we also performed molecular docking of top-hit compound (Mol-1388) and it is compared with Olaparib at the binding pockets of PARP2 (PDB ID: 4TJV). Results show that docking scores of Mol-1388 and Olaparib are -6.11kcal/ mol and -14.92kcal/mol. Thus, the identified piperazine derivative Mol-1388 is found more selective than Olaparib based on docking scores.Olaparib is an FDA-approved targeted therapy for cancer and it is the FDA approved inhibitor of PARP1. This compound also includes piperazine ring in its structure. Piperazine has received considerable interest in the development of many effective biologically active compounds. 4674 compounds that include piperazine ring from Specs SC database (among more than 210.000 compounds), were screened at the binding pocket of PARP1 target (Figure S4, Supplementary material). Thus, we have studied with piperazine derivatives and considered Olaparib as positive control. In order to compare the fingerprint of identified hit molecule and olaparib at the binding pocket of the PARP1, these two structures were superimposed at the binding site of the enzyme (Figure S5, Supplementary material). Main common residues that are important in both of Olaparib and hit molecule Mol-1388 were Tyr896 and Tyr907.MetaCore/MetaDrug (MC/MD) platform of Clarivate Analytics includes binary QSAR models for activity predictions, pharmacokinetic (ADME) and toxicity estimation for small compounds . MC/MD contains individually designed disease-QSAR models for the 25 common diseases, including cancer, and estimates the potential therapeutic activity of the screened small molecule with normalized values between 0 and 1 based on structural similarity. Estimated QSAR values greater than 0.5 indicate potential therapeutic activity. Here, we used this platform to check the therapeutic activity prediction of hit compound mol-1388 using “cancer-QSAR” model. Cancer-QSAR model of MC/MD is created with 886 training set compounds and 167 test set compounds. The model’s sensitivity, specificity, accuracy, and Matthews Correlation Coefficient (MCC) values were found as 0.89, 0.83, 0.86 and 0.72, respectively. The obtained statistical results show high accuracy and high predictive power of the model. Model predicts therapeutic activity potentials of screened compounds with normalized values between 0 (inactive) and 1 (active). Cancer-QSAR model in MC/MD predicted the therapeutic activity potential of Mol-1388 as 0.60. This compound was also screened in 26-different toxicity models and it is found that compound does not show any toxicity within these 26 toxicity models.

4. Conclusions

In this work, we performed advanced text-mining and combined molecular modeling approaches to identify the piperazine-based novel PARP1 small molecule inhibitors. For this purpose, 4674 compounds that include piperazine were selected from Specs-SC small molecule database. These obtained compounds were screened against the PARP1 target using molecular docking and the top-docking poses of selected 10 hits were then used in long (100ns) MD simulations. For all ligands, MM/GBSA binding free energy calculations were carried out. The same procedure has also been applied for the FDA approved PARP1 inhibitors. Since many FDA approved anti-cancer drugs including PARP1 inhibitor Olaparib includes piperazine ring, in this study we used the advantage of text mining approaches which can decrease the total number of compounds in short time periods from Disseminated infection small molecule large databases to reasonable numbers of compounds for further analyses. Eventually, throughout this strict combination screening we investigated hit compounds that may be used to inhibit PARP1 activity. This new scaffold may open up new avenues for the designing of small inhibitors against PARP1.