In silico investigation of PARP-1 catalytic domains in holo and apo states for the design of high-affinity PARP-1 inhibitors

Abstract

The rational design of high-affinity inhibitors of poly-ADP-ribose polymerase-1 (PARP-1) is at the heart of modern anti-cancer drug design. While relevance of enzyme to DNA repair processes in cellular environment is firmly established, the structural and functional understanding of the main determinants for high-affinity ligands controlling PARP-1 activity is still lacking. The conserved active site of PARP-1 represents an ideal target for inhibitors and may offer a novel target at the treatment of breast cancer. To fill the gap in the structural knowledge, we report on the combination of molecular dynamics (MD) simulations, principal component analysis (PCA), and conformational analysis that analyzes in great details novel binding mode for a number of inhibitors at the PARP-1. While optimization of the binding affinity for original target is an important goal in the drug design, many of the promising molecules for treatment of the breast cancer are plagued by significant cardiotoxicity. One of the most common side-effects reported for a number of polymerase inhibitors is its off-target interactions with cardiac ion channels and hERG1 channel, in particular. Thus, selected candidate PARP-1 inhibitors were also screened in silico at the central cavities of hERG1 potassium ion channel.

• MD simulations
• molecular modeling
• PARP-1

Introduction

Poly [ADP-ribose] polymerase 1 (PARP-1) is an archetypical enzyme from the PARP group, which contains 18 identified proteins in humans1. PARP-1 catalyzes the polymerization of ADP-ribose units from NAD+ co-enzymes (as donor group) on target proteins2. An entire family shares well-conserved catalytic domain with “PARP signature” motif. PARP-1 structural and functional organization is also well preserved with four major domains: (i) zinc-finger DNA-binding domain (DBD), (ii) caspase-cleaved domain, (iii) auto-modification domain, and (iv) a C-terminal catalytic domain containing an active site for the enzyme3^,4. Because of its significance in control of DNA repair, PARP family has been considered as an important target for anti-cancer therapeutics5. It was suggested that a number of PARP-1 inhibitors could be effective against tumors with inherent DNA repair defects, such as common breast tumors with mutations in the DNA repair proteins BRCA1 or BRCA2. PARP inhibitors also showed considerable promises as part of therapeutic intervention in several other diseases, such as stroke, cardiac ischemia, cancer, inflammation, and diabetes6–9. PARP inhibitors can be used as anti-cancer agents if a genetic treatment is considered; indeed they are already started to use for the treatment of triple-negative breast cancer10. PARP catalyzes the cleavage of NAD+ into ADP and ADP-ribose. Expectably, most of the developed inhibitors resemble the nicotinamide moiety of NAD+ to achieve high-affinity binding to the corresponding catalytic pocket in PARP-110. The available X-ray data involving the catalytic domain of PARP-1 with novel bound inhibitors at the binding site leads virtual screening studies on PARP-1 inhibitors11.

Here molecular dynamics (MD) simulations are performed for better understanding of the inhibitor interactions with PARP-1 enzyme. All-atom MD simulations were then used for elucidation of inhibitor and enzyme dynamics. In the present work, structural and dynamical profiles of both target and inhibitors are thoroughly analyzed with derived MD trajectories. The functionally important transitions and structural features were defined with principle component analysis (PCA) and clustering analysis of MD trajectories. The data from MD simulations lead us an opportunity to explore the intermolecular ligand/protein interactions at atomic scale. In particular, MD simulations enabled comparison of the structural and dynamical re-arrangements between holo and apo states of PARP-1. A set of known drugs on the market is then tested to map their multi-functional potentials as potent PARP-1 inhibitors. The well-known feature of many polymerase inhibitors is their high cardiotoxicity due to blockade of hERG1 channel. It is well established that human ether-à-go-go related gene 1 (hERG1) K channel is overexpressed in many tumors and, therefore, pro-arrhythmic potential of all novel inhibitors is an important concern in structure driven drug development. To circumvent potential risks associated with drug-induced QT prolongation, we report on an in silico cardiotoxicity screening against previously developed hERG1 models12–19 by our group.

Methodology

Input receptor structure

Crystal structure of PARP-1 (Protein Data Bank (PDB) code: 4HHY)11 which involves a bound-novel inhibitor inside the binding site was taken from the X-ray structure (Figure 1a). Hydrogen atoms were added to the receptor by employing the Protein Preparation Wizard20 integrated in the Maestro molecular modeling suite of Schrodinger software package (Tripos, St. Louis, MO). The PROPKA21 program was used to assign the protonation states of amino acids as in physiological medium, i.e. the protonation and tautomeric states of Asp, Glu, Arg, Lys, and His residues were assigned at the pH of 7.4. The hydrogen atoms were then relaxed with a restrained minimization procedure.

Figure 1. (a) Chemical structure of PARP-1 inhibitor used in MD simulations. (b) The equilibrated system includes PARP-1 protein (brown ribbons), inhibitor (CPK model), ions (green and yellow spheres), and water molecules (white surface).

Simulation protocol

The protein–ligand complex was simulated for 25 ns in 0.15 M KCl solution and in the presence of explicit water molecules. The TIP3P water model was employed. The aqueous medium was modeled by Visual Molecular Dynamics (VMD)22 program. MD simulations were carried out using all-atom CHARMM36 force-field parameters with CMAP corrections23 implemented in NAMD 2.9 program package (Theoretical Biophysics Group, Baltimore, MD)24. All bonds involving hydrogen atoms were restrained using SETTLE25 algorithm. The Langevin piston method with Langevin thermostat was used to maintain constant temperature (310 K) and constant pressure (1 atm) during simulations. Periodic boundary conditions (PBC) and particle mesh Ewald (PME) methods were used to ensure the correct treatment of long-range electrostatic interactions. The non-bonded cut-off was set as 10 Å for inter-atomic Lennard–Jones (LJ) and Coulombic interactions. Water molecules around protein were equilibrated for 1 ns by keeping the protein atoms at their fixed positions then the constraints on the atomic positions of protein atoms were reduced gradually. The MD trajectories were analyzed utilizing VMD, PyMOL26, Bio3D27, and R-based scripts. The snapshot picture of the simulation box at the end of the equilibration is given in Figure 1(b).

Principal component analysis

Principal component analysis (PCA)28 is a statistical approach to identify stable basins in high-dimensional data allowing for reducing dimensionality of the process. In order to implement PCA in MD simulations, 2000 trajectory frames derived throughout the MD simulations are incorporated into the Bio3D scripts. A binary MD trajectory file includes coordination data of protein atoms (Cα atoms used in this study). The coordinates of atoms were written in the column of matrix and they are aligned with respect to the initial conformer to remove the effect of diffusion of atoms during MD simulations in PCA.

Mathematical calculations of PCs from MD trajectories include two consequence steps:

• (i) generation of the covariance matrix (C), (1)

• (ii) the next stage is diagonalization of the 3N × 3N covariance matrix C that can be applied by eigenvalue decomposition (EVD) as (2)

where Λ indicates the diagonal matrix involving the eigenvalues and V indicates a matrix involving eigenvectors.

In order to identify the principal dynamics of the MD simulations and separate it from the unrelated noise, we applied a PCA for the MD frames of the two systems by diagonalizing the positional fluctuations of covariance matrix (S) of the Cα atoms relative to average structure of a given conformation. We also performed dynamic cross-correlation maps (DCC) to show dynamic cross-correlated motions of Cα atoms in the trajectory re-oriented relative to an average structure.

MD simulations of biomolecules result in a collection of noisy states, which poses a considerable challenge in finding principles governing ligand binding to the receptor. A PCA can reduce complexity of the protein dynamics by mapping correlations existing in the trajectory.

Per-residue interaction energy calculations, and hydrogen bonding network

Per-residue interaction energy calculations, contact maps, and hydrogen bonding analyses were performed with the available plugin in VMD. The per-residue contributions to protein–ligand interaction energy were evaluated using a cut-off distance of 12 Å. The following criteria were used for the evaluation of hydrogen bonding: the angle cut-off values were <20° and the heavy atom–light atom distance cut-off was 3.0 Å.

Molecular docking simulations

A set of approved drugs was retrieved from Drug Bank database29. The lowest energy conformers of the ligands were searched using LigPrep module30 of Maestro Suite31. The ionization states and the partial charges of the ligand atoms were assigned according to the pH of 7.4. Molecular docking studies were performed with the Glide docking32 program. Glide has three scoring functions: extra-precision (XP), standard precision (SP), and high-throughput virtual screening (HTVS). In this study, initially the SP method was performed for large data set of inhibitors and top-docking poses were selected and re-docked using the XP method. The volume of the inner and outer active binding-site of PARP1 target was defined as (14 × 14 × 14) ų and (30 × 30 × 30) ų, respectively. The active site was defined by the centroid rendered from inhibitor observed in crystallographic structure. For top-10 docking poses generated by Glide/XP, a multi-stage induced fit docking (IFD) method was also performed33. For IFD, the van der Waals radius rendering for the complex was set by a factor of 0.5 and ligands were docked with the default Glide/SP. Then, Prime34 protocol was used to refine positioning of the protein side chains. Finally, the ligand poses were scored and top-poses were re-docked using standard Glide/XP without any scaling applied to van der Waals radius.

Results and discussion

In this study, a combined approach involving docking simulations, MD simulations, and PCA was applied to identify the key structural features governing effective inhibition of PARP-1. MD simulations performed on two different conformational states of PARP-1, named as apo and holo, were also used to compare state-dependent inhibitory mechanism. At the end of 25 ns MD simulations, the total energy, temperature, and volume fluctuations of the system reached to a plateau indicating the well-equilibrated system. Figure 2 compares the equilibrium structure obtained from MD simulations with the X-ray structure of PARP-1. The bound inhibitor displays very little relative displacement in the course of MD simulations relative to its position in the crystal structure (Figure 2).

Figure 2. (a) The crystal structure of the PARP-1 complex (green) superimposed with the representative structure (red). (b) Equilibrated inhibitor pose and its interactions with active site residues.

Ligand interaction diagram was used to show 2D receptor–ligand interactions in a representative MD trajectory (Figure 3). Five stable hydrogen-bonding interactions which were reported to be essential for drug stabilization in PARP-1-ligand complexes were reproduced in our simulations. These H-bonds are (i) between the inhibitor carbonyl's oxygen and the OH group of Ser904 amino acid; (ii) between the nitrogen atom in the ring of the ligand and backbone oxygen of the Gly863; (iii) between the amide backbone of Gly863 and carbonyl oxygen of the ligand; (iv) between the amide backbone of Ile879 and carbonyl oxygen of the ligand; (v) between carbonyl oxygen of the ligand and backbone OH group of the Arg878. Additionally, the number of hydrophobic interactions between the inhibitor and protein residues (i.e. Phe897, Trp861, Tyr907, Ala880, Ile872, and Ile895) are available.

Figure 3. Mainly observed conformations of inhibitor at the binding pocket throughout 25 ns MD simulations. 2D ligand interaction diagram for the representative complex structure has also been shown in the figure.

Clustering analysis

Clustering analysis in PARP-1–ligand complex allowed us to identify the stable binding poses of the ligand. The inhibitor in the active site cavity has a dominant conformation which comprises the 99% of frames in analyzed trajectories as shown in Figure 3.

System stability analysis

The RMSD values for Cα atoms can serve as a measure for relative flexibility of the binding pockets in apo and holo states of PARP-1 (Figure 4a). To illuminate the effect of a given conformation on the stability of binding pocket, we calculated the RMSD values with respect to Cα atoms of the residues comprising binding site and participating in hydrogen bonding (Figure 4b) using the protein X-ray diffraction structure. Two independent simulations with and without bound inhibitors, i.e. at the holo and apo states, were carried out to determine the effect of inhibitor on the protein dynamics. The RMSD values for Cα atoms of residues forming binding pocket showed a slightly higher RMSD in the apo state than in the holo state at the equilibration stage of the MD simulations. After the equilibration, the PARP-1 structure at both apo and holo forms did not display any significant conformational transitions (Figure 4). The RMSD values for all Cα atoms of the PARP-1 protein at apo and holo states fluctuated around a mean value (1.83 Å and 1.65 Å for holo and apo forms, respectively) and remained almost the same during the equilibration (Figure 4).

Figure 4. (a) RMSD versus time plot for all Cα atoms at the target; (b) RMSD versus time plot for Cα atoms of residues that are participating in hydrogen bonding with the inhibitor.

Figure 5 illustrates per-residue RMSF values for Cα atoms of the protein. The Cα atoms of the residues which involve in loops have higher flexibility than the ones in the α-helices and β-sheets at both apo and holo forms. The incorporation of the inhibitor induces some conformational changes as reflected in the regions formed by residues from 696 to 713 and also from 881 to 888. The RMSF plot confirms the experimentally observed H-bond interaction between the amino acids, Ser904, Gly863, and inhibitor at the active binding site at both apo and holo proteins. As expected, the apo form is more flexible than the holo form, inhibitor-bearing structure (Figure 5).

Figure 5. RMSF of Cα atoms of all residues at apo and holo states of PARP-1 protein during the MD simulations. The arrows indicate the residues participating in hydrogen bonding in the holo form.

Per-residue interaction energy analysis

The quantitative analysis of the interactions between the inhibitor and the amino acids of the residues in the binding pocket was done by calculating the interaction energies as an average of 25 ns MD simulations time. The interacting groups are categorized according to their position with respect to the inhibitor molecule. For better visualization, the structure of the inhibitor is divided into three segments and the interactions of the each segment with the specific amino acids, and their binding energies are shown in Figure 6. The total non-bonded interaction energies including the H-bond of the first, second, and third segments are found to be −24, −14, and −16 kcal/mole, respectively. The first segment interacts via strong hydrophobic and weak π–π stacking interactions with Tyr907 and Tyr896, via H-bonding interactions with Ser904 and Gly863. The chemical structure of the first segment resembles the active pocket of nicotinamide adenine dinucleotide (NAD⁺) for which the strong H-bond interaction is reported35^,36. The non-bonded interactions are mostly van der Waals type and relatively stronger interactions involve His862, Tyr907, Tyr896, and Arg878 amino acids with the energies varying between −4 and −9 kcal/mole. This computational result is in a good agreement with the experimental finding which pointed out that His862 and Glu988 have significant contributions to the interactions at the active site35^,36. Among the PARP-1 active site amino acids, Glu988 is important for NAD⁺ catalysis37^,38, Lys903 is critical for the DNA polymerization process39. The second part of the inhibitor molecule forms strong hydrophobic interaction with His862 of beta sheet-like Glu988, it also participates in NAD⁺ binding40. The third segment of the inhibitor interacts with Ala880, Leu769, Ile879, Pro881, and Arg878 which plays a critical role in stabilizing the adenosine portion of the substrate NAD⁺41.

Figure 6. The average ligand–target interaction energies between the given residues and the inhibitor.

H-bond formation ratios during the MD simulations for the key-amino acids have been represented in Figure S1 in the Supplementary materials. The percentage of hydrogen bond occupancies of the key-amino acids during 25 ns MD simulations has been found as Gly863 (52%), Ser904 (33%), Arg878 (40%), and Ile879 (30%) (Figure S1, Supplementary materials).

PCA

Eigenvalues defining collective modes of fluctuations in apo and holo states of PARP-1 protein are shown in Figure 7. Residue cross-correlation calculations allowed us to examine the effect of ligand on the conformational dynamics of the protein. It is clear from these calculations that inhibitor binding has an impact on a low-frequency harmonic mode propagating along an entire PARP-1 structure. Figure 8 illustrates positive (moving of amino acid residues in the same directions) and negative (moving of amino acid residues in opposite directions with respect to each others) correlations for studied systems. Positive and negative correlations between amino acid residues conclusively determine the conformational transition pathway during inhibitory treatments of protein. Higher negative correlations were significantly established between PARP-1 amino acids in the presence of inhibitor, the corresponded domains and residues have been indicated in correlation profile.

Figure 7. Comparison of PCA profiles of apo and holo proteins.

Figure 8. Cα-residue cross-correlation profiles for apo and holo proteins. Blue and red parts indicate the positive and negative correlations, respectively.

Protein conformation clustering

To elucidate the effects of drug binding on conformational dynamics of the protein, we employed RMS-based clustering analysis using Cα positions42. Reduction of the protein dynamics to Cα dynamics allows us to decrease the dimensionality of the system of interest by eliminating fast side-chain dynamics. Figure 7 shows the PC profile and displays the relation of variance versus applied eigenvalue mode at apo and holo states. The first four modes, i.e. PC1, PC2, PC3, and PC4, combine 43.0% and 55.1% of the variance for apo and holo forms, respectively. Unlike to the apo state which includes low occupancy for the PC1 (17.5%), the holo form of the first PC indicates a high-degree linear conformational transition with respect to the inhibitor (37.3% of the whole conformational transition). Highly fluctuated regions during MD simulations are analyzed using PCA and represented in Figure S2 in the Supplementary materials.

Molecular docking

Database mining allows rapid and exhaustive search of potential candidates for further biological testing. To improve its performance on the rational design of drug candidates, it can be matched to an array of De Novo drug design algorithms43. To initiate database development, we focus on an array of small ligands with known PARP-1 inhibition properties to create our initial hypothesis. The available drug-bank databases are analyzed and ∼9000 small molecules are considered as potential inhibitor molecules of PARP-1 target. Their properties such as molecular weight, ALogP, hydrogen bond numbers as acceptors and donors, molar refractivity, number of rotatable bonds, and polar surface area were used to filter out the ligands with a low degree of drug-likeness44. The remaining compounds are then subjected to Glide/SP docking protocol. Docking scores of compounds which are better than −8 kcal/mole were selected as potential PARP-1 inhibitors. The Glide/XP and Glide/IFD protocols were then employed for the selected ligands as given in Figure S3 in the Supplementary materials.

In addition to test the validation of used docking protocols in examining the PARP-1 inhibitory activity of compounds, several known active and inactive PARP-1 inhibitors were selected and docked into the binding pocket of the target structure (Table 1). Each compound was identified as its ZINC database ID, experimental binding affinities (IC₅₀)45–52, and estimated docking scores were reported for individual compounds. While docking results showed higher (absolute) docking scores for PARP-1 active inhibitors, inactive PARP-1 molecules showed lower docking scores (Table 1). These results significantly suggest the usefulness of low time-consuming in silico methods in discovering of novel potent inhibitors.

Table 1. Experimental binding affinities and estimated docking scores (Glide/XP) of active and inactive PARP-1 inhibitors at interaction with PARP-1 target. Each compound is identified by its ZINC database ID.

Compound (inactive)	Docking score (kcal/mol)	Binding affinity (experimental IC50 (nM))	Compound (active)	Docking score	Binding affinity (experimental IC50 (nM))
ZINC02041009	−4.78	8500	ZINC13860431	−11.34	100
ZINC01590754	−5.00	9100	ZINC53297239	−11.05	10
ZINC17970819	−5.17	9500	ZINC03821234	−8.02	19
ZINC00150158	−5.20	8100	ZINC28467935	−7.87	1
ZINC00152996	−5.21	5400	ZINC43206370	−6.86	3
ZINC01512423	−5.26	2010	ZINC34031449	−6.79	3
ZINC08649468	−5.42	6100	ZINC13536325	−6.74	20

MD simulations53 and SAR studies54 on PARP-1 inhibitor complexes revealed that the potent PARP-1 inhibition depends on the position of the carboxamide oxygen. This functional group is stabilized by hydrogen bonds that it forms with Ser904 and Gly863 in the PARP-1 NAD+ binding site. According to our docking results, the residues in the ligand binding-site included Ser904, Gly863, and also Arg865, Glu988, Tyr907, and Ser864.

Selected PARP-1 potent inhibitors were also tested in silico against their cardiotoxicity risks. They were docked into the central cavities of previously generated and validated human ether-a go-go related gene (hERG1) K-channel models by our group12–19. Table 2 shows the Glide/XP docking scores of these drugs in hERG1. Figure 9 shows top-docking pose of carbergoline which has the lowest docking score at the pore domains of hERG1 K channel. Our results are in agreement with the experimental studies which indicate that cabergoline is a prolactinemic drug with low toxicity. It is an anticancer drug and controls the prolactin release in metastatic breast cancer55^,56. The second identified anti-cancer drug is raloxifene. According to in vitro studies, raloxifene decreases and prevents breast cancer disease in women57. Idarubicin drug is characterized with a high affinity against PARP-1 enzyme and low hERG blocking properties. Idarubicin is a member of a generic group of anticancer drugs named as anthracycline antibiotics58.

Figure 9. 2D ligand interaction diagram of cabergoline with the hERG pore domain and its 3D top-docking pose of at hERG1 K channel.

Table 2. Docking scores of selected 12 potent drugs against hERG1 channel pore domain (kcal/mol).

Generic name	Glide/XP	Generic name	Glide/XP
S-Adenosylmethionine	−6.11	Arbutamine	−8.78
Masoprocol	−8.50	Lisdexamfetamine	−7.58
Idarubicin	−5.98	Cabergoline	−4.63
Raloxifene	−5.44	Bimatoprost	−7.69
Propafenone	−7.41	Sunitinib	−6.15
Nebivolol	−8.00	Hydroxychloroquine	−6.13

The selected molecules which show high docking scores at the active site of the PARP-1 target are re-screened at the hERG1 pore domain for checking in silico their cardiotoxicity risks. Molecules which show low predicted hERG1-blocking affinities at the central cavities of the channel can be considered as safer drugs.

Conclusions

In this study, we have combined molecular docking, MD simulations, and binding energy analysis on both holo and apo forms of PARP-1 target to obtain some insights into the conformations and interactions for binding of PARP-1 inhibitors. Our work indicated that the residues Ser904, Gly863, and Arg878 forming hydrogen bond with the ligand are highly responsible for the stability of the binding pocket. Especially, Ser904 and Gly863 are primarily corresponded for establishing the inhibitory potency of PARP-1 inhibitors. PC analysis merged with clustering to elucidate conformational transitions in the presence of PARP-1 inhibitor. PCA analysis showed that ligand binding to PARP-1 hinders protein dynamics compared with the ligand-free state. We also screened in silico a set of known compounds as PARP-1 inhibitor candidates both in principal target and against common off-target (central cavities of hERG1 K channel).

Supplementary material available online

Supplementary Figures S1-S3.

Acknowledgements

Part of computations for the work described in this paper was supported by Turkish Scientific and Technical Research Council (TUBITAK) ULAKBIM High Performance Computing Center. S.D. acknowledges support from Bilim Akademisi – The Science Academy, Turkey, under the BAGEP program.

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of this article.