Revealing the selective mechanisms of inhibitors to PARP-1 and PARP-2 via multiple computational methods

Background Research has shown that Poly-ADP-ribose polymerases 1 (PARP-1) is a potential therapeutic target in the clinical treatment of breast cancer. An increasing number of studies have focused on the development of highly selective inhibitors that target PARP-1 over PARP-2, its closest isoform, to mitigate potential side effects. However, due to the highly conserved and similar binding sites of PARP-1 and PARP-2, there is a huge challenge for the discovery and design of PARP-1 inhibitors. Recently, it was reported that a potent PARP-1 inhibitor named NMS-P118 exhibited greater selectivity to PARP-1 over PARP-2 compared with a previously reported drug (Niraparib). However, the mechanisms underlying the effect of this inhibitor remains unclear. Methods In the present study, classical molecular dynamics (MD) simulations and accelerated molecular dynamics (aMD) simulations combined with structural and energetic analysis were used to investigate the structural dynamics and selective mechanisms of PARP-1 and PARP-2 that are bound to NMS-P118 and Niraparib with distinct selectivity. Results The results from classical MD simulations indicated that the selectivity of inhibitors may be controlled by electrostatic interactions, which were mainly due to the residues of Gln-322, Ser-328, Glu-335, and Tyr-455 in helix αF. The energetic differences were corroborated by the results from aMD simulations. Conclusion This study provides new insights about how inhibitors specifically bind to PARP-1 over PARP-2, which may help facilitate the design of highly selective PARP-1 inhibitors in the future.


INTRODUCTION
As one of the most commonly diagnosed malignancies in women, breast cancer accounts for about a quarter of all female cancer cases (Siegel, Miller & Jemal, 2018). In the past few years, the incidence of breast cancer has continued to rise, with more than 1 million new cases worldwide each year (Siegel, Miller & Jemal, 2018). It is estimated that a total of 5-10% of all breast cancer cases are genetically susceptible to the disease, with multiple breast cancer susceptibility genes having been proposed, including breast cancer 1 (BRCA1) and breast cancer 2 (BRCA2), two of the major genes (Begg et al., 2008). Currently, sequencing of these two genes is considered as the optimal approach to determining the mutation status in breast cancer patients (Pfeffer, Ho & Singh, 2017). Previous studies have shown that homologous recombination-deficient tumor cells resulting from BRCA1 or BRCA2 gene mutations are hypersensitive to the inhibitory effect of poly-ADP ribose polymerase-1 (PARP-1) (Lin & Kraus, 2017;McCann, 2019). One mechanism that has been proposed as a tentative explanation is that inhibition of PARP-1 blocks DNA single strand break repair and leads to the formation of unrepaired double strand breaks at the replication fork (Kim et al., 2020;Min & Im, 2020;Wang, Luo & Wang, 2019b). Aside from DNA damage repair, PARP-1 is also involved in a wide variety of cellular processes, such as cell proliferation and cell death. The implication of PARP-1 in these processes is a result of the diverse substrates in PARP-1, such as nuclear proteins, which are involved in apoptotic cell death, cell cycle regulation, chromatin decondensation, inflammation, and transcriptional regulation (Jubin et al., 2016). Due to these functions, PARP-1 inhibitors have been developed as the first class of cancer therapeutics in clinical trials (Min & Im, 2020).
Currently, four PARP-1 drugs have been approved by the US Food and Drug Administration (FDA), including olaparib, niraparib, talazoparib, and rucaparib. In addition to these commercially available PARP-1 drugs, numerous PARP-1 inhibitors have also entered different phases of clinical research, targeting multiple types of tumor either collectively or as single agents (Min & Im, 2020). Most of the marketed PARP-1 drugs and inhibitors exhibit poor selectivity when targeting PAPR-1. For instance, olaparib, the first PAPR-1 drug approved by FDA, also interacts with PAPR-1 close homologues PARP-2 and PARP-3 (Gunderson & Moore, 2015;Min & Im, 2020). Similarly, Niraparib, Talazoparib and rucaparib were also found to exhibit no selectivity between PARP-1 and PARP-2 (Min & Im, 2020). There is ample evidence in previous studies that inhibition of PARP-2 could produce potentially undesirable side effects (Farres et al., 2013;Navarro et al., 2017). For instance, Farres et al. (2013) reported that loss of PARP-2 leads to a shortened red blood cell lifespan and impaired differentiation of erythroid progenitor cells, thereby causing chronic anemia. In light of this situation, a great deal of efforts has been put into the design and development of potent PARP-1inhibitors with high selectivity, especially between PARP-1 and PARP-2 (Eltze et al., 2008;Fatima et al., 2014;Papeo et al., 2015). However, given high sequence similarity (84% identity and 90% similarity) and conserved catalytic domain (Figs. 1A-1C), increasing the selectivity of inhibitors remains a huge challenge (Yelamos et al., 2011).
Up to now, some potent inhibitors with high selectivity to PARP-1 over PARP-2 have been discovered, such as WD2000-012547, BYK204165 and NMS-P118 (Eltze et al., 2008;Fatima et al., 2014;Papeo et al., 2015). Despite these promising results, little computational research has been conducted to elucidate the selective mechanisms underlying the effect of these inhibitors. To this end, two representative inhibitors (Niraparib, NMS-P118) with divergent selectivity to PARP-1 and PARP-2 were utilized in this study to demonstrate such mechanisms (Ison et al., 2018;Papeo et al., 2015). Niraparib (also formerly known as MK-4827) is a novel, highly selective, and orally available PARP-1 and PARP-2 small molecule drug developed by Tesaro, approved by FDA in 2017 to treat ovarian cancer  Ison et al., 2018). Recent clinical studies show that the combined use of niraparib and pembrolizumab can have promising antitumor effects on patients with advanced or metastatic triple-negative breast cancer (Lin, Yang & Zhao, 2019). NMS-P118, another inhibitor originally designed by Papeo et al., has been shown to exhibit ∼150-fold selectivity to PAPR-1 over PARP2 . In addition, the researchers also reported the co-crystal structure of PARP1 and PARP-2 bound to NMS-P118, and proposed that helix αF in the two proteins may be responsible for the drug selectivity, as the induced binding pocket is larger in PARP-1 than in PARP-2 (Papeo et al., 2015). However, this explanation seems far from satisfactory. Here, classical molecular dynamics (cMD) simulations and accelerated molecular dynamics (aMD) simulations were employed to clarify the selective mechanisms between PARP-1 and PARP-2 using two representative inhibitors (Niraparib, NMS-P118). A combination of classical MD simulations, the root-mean-square deviations (RMSD), principal component analysis (PCA), dynamical cross-correlation (DCC) analysis, and root-mean-square fluctuations (RMSF) was applied to investigate the effects of the inhibitors on the protein flexibility and dynamic behavior of key parts of PARP-1 and PARP-2. Afterwards, binding free energy calculations and per-residue free energy decompositions based on the molecular mechanics/generalized Born solvent area (MM/GBSA) method were performed to highlight key residues related to selectivity. Next, aMD simulations combined with RMSD, PCA, DCC and free energy landscape (FEL) analyses were carried out to examine in detail the local energy minima and the conformational space that were not illuminated in the classical MD simulations. Overall, these results can be effective to deepen our understanding of the selective mechanisms between PARP-1 and PARP-2, and may help facilitate the design of novel inhibitors to improve drug selectivity.

Preparation of the initial systems
The three-dimensional structures of human PARP-1 bound to NMS-P118 (PDB ID: 5A00), Niraparib (PDB ID: 4R6E), and PARP-2 bound to NMS-P118 (PDB ID: 4ZZY) were obtained from the Protein Data Bank (Papeo et al., 2015;Thorsell et al., 2017). The Loops/Refine Structure module of UCSF Chimera program was employed to model the missing side-chains and loop structures (Pettersen et al., 2004). The PDB2PQR Server was employed to estimate the protonation states of ionizable side chains (Dolinsky et al., 2004). The initial coordinates of PARP-2 bound to Niraparib were constructed using the AutoDock program (Morris et al., 2009). The grid size of cubic box, which was centered on the binding pocket, was set to 60 ×60 ×60 xyz points with a grid spacing of 0.375 Å. The AutoDockTools program was employed to assign the Gasteiger partial charges to PARP-2 and Niraparib. The affinity maps of grids were estimated using AutoGrid program. The docking protocol involved the generation of 200 conformations. The maximum number of energy evaluations and iterations were set to 25,000,000 and 3,000, respectively. Other parameters were set to default. The top-ranked structure was used for the subsequent molecular dynamics (MD) simulation analyses.

Classical MD simulation
The prepared crystal structures of PARP-1 bound to NMS-P118 and Niraparib, PARP-2 bound to NMS-P118, and modeled complex of PARP-2 bound to Niraparib were applied to determine the dynamic structural behavior via Assisted Model Building with Energy Refinement 18 (Amber 18) program. The restrained electrostatic potential (RESP) fitting technique was employed to estimate the partial atomic charges of NMS-P118 and Niraparib (Wang, Cieplak & Kollman, 2000). The parameters of protein and ligand were derived from the ff14SB force field the General Amber Force Field 2 (GAFF2) in Amber 18 (Maier et al., 2015;Vassetti, Pagliai & Procacci, 2019). Each of the prepared complexes was solvated in a cubic box containing TIP3P water molecules, with a minimum distance of 15 Åfrom any edge of the box to any complex atom. Counter ions of an appropriate quantity were added to the system to preserve overall charge neutrality.
Prior to the classical MD simulation, two-step minimizations, heating and equilibration were performed. At first, two-step minimizations were undertaken to eliminate bad contacts between the solvent molecules and the complexes. To reduce the counterions and water molecules to a minimum, a harmonic constraint of 20 kcal mol −1 Å −2 was first imposed on the four complexes. Then, restriction was eliminated in order for all atoms to move freely. During each stage, the steepest descent minimization of 7,000 steps was performed, followed by conjugate gradient minimization of 7,000 steps. Thereafter, Langevin thermostat with a position restraint of 20 kcal mol −1 Å −2 was applied to gradually heat up each complex from 0 K to 300 K over 300 ps. Then, each complex was equilibrated at 300 K with 1000 ps simulation time in the isothermal isobaric (NPT) ensemble. Finally, 800 ns production classical MD simulation was carried out for each complex in the NPT ensemble with a time step of 2fs. During the simulations, the Langevin temperature scalings and Berendsen barosta were utilized to maintain the temperature and pressure, respectively (Izaguirre et al., 2001;Praprotnik, Delle Site & Kremer, 2005). The Particle mesh Ewald (PME) method was employed to estimate the long-range electrostatic interactions, with the cutoff parameter of nonbonded interaction set to 10 Å (Essmann et al., 1995). SHAKE method was applied to constrain all covalent bonds connecting hydrogen atoms (Krautler, VanGunsteren & Hunenberger, 2001). The coordinates for each complex were saved at an interval of 10 ps for subsequent analysis. The RMSDs and RMSF of the trajectories were calculated using CPPTRAJ module in Amber 18 program.

aMD simulations
The aMD is an enhanced sampling technique that alters the energy landscape through the addition of a boost potential V(r) to the original potential energy surface. The V(r) stands for either the total potential energy (E total ) or the dihedral energy (E dihedral ) of a system (Hamelberg, Mongan & McCammon, 2004). When V(r) is equal or greater than a predefined reference energy E, none of addition energy will be added (Eq. (1)). On the contrary, the modified V (r) of the system is calculated according to the following Eq. (2): where α is the acceleration factor that governs the size of the boost. The α is calculated with Eqs.
(3) and (4), and the E is calculated according to Eqs. (5) and (6). The boost parameters E and α for the total boost (E total and α total ) and dihedral boost (E dihedral and α dihedral; ) are based on the corresponding average E total (V total_avg ) and average E dihedral (V dihedral_avg ), which are calculated from the classical MD simulations prior to the aMD simulations. N atoms and N res represent the number of atoms and residues in the system, respectively. In Eq.
(3), the n is an integer defined as the magnitude of the threshold, which is a multiple of the α.
Herein, the equilibrated structures extracted from classical MD simulations were selected as the initial structures for the aMD simulations. The Dual-boost approach was applied by adding V(r) to both the E total and E dihedral of the system. The V(r) was obtained based on the N atoms , N res , V total_avg , and V dihedral_avg from the first 40 ns of the classical MD simulations. Then, 800 ns aMD simulations were employed using the dual-boost approach. During aMD simulations, the PME method was employed to estimate the long-range electrostatic interactions, with the cutoff parameter of nonbonded interaction set to 10 Å (Essmann et al., 1995). SHAKE method was applied to constrain all covalent bonds connecting hydrogen atoms (Krautler, VanGunsteren & Hunenberger, 2001). The Langevin temperature scalings was used to handle the temperature (Izaguirre et al., 2001). The coordinates for each complex were saved at an interval of 10 ps and the trajectories were calculated using the CPPTRAJ module in Amber 18 program.

Principal component analysis (PCA)
PCA was performed on the trajectories from both classical MD and aMD simulations via the CPPTRAJ module in Amber 18 program. All snapshots of each trajectory were aligned to eliminate the translational and rotational motions of all protein C α atoms. Then, a covariance matrix (3N × 3N) was generated from the Cartesian coordinates. A set of eigenvectors and eigenvalues were generated by diagonalizing the matrix. The top two eigenvalues (principal component 1 and principal component 2, PC1 and PC2) were used for subsequent analysis.

DCC analysis
The Bio3D package of R was employed to conduct DCC analysis of the backbone atoms (C α ) (Skjaerven et al., 2014). The cross-correlation matrix (C ij ) between residues i and j was generated based on the trajectories from both classical MD simulations and aMD simulations, with 5,000 snapshots for each complex. The C ij is calculated based on the following Eq. (3). The r i and r j represent the shift from the mean position of the ith or jth of protein C α atom. angle bracket represents an average of these two values over the sampled period.

Binding free energy calculations based on classical MD simulations
The MM/GBSA method is often used to perform classical MD simulations combined with free-energy calculations, which can serve as an effective tool for quantitative prediction of protein-ligand binding energies (Wang et al., 2019a). The binding free energies ( G bind ) are calculated on the basis of the MM/GBSA method according to the following equations: In Eq. (8), G R+L , G R , and G L stand for the free energies of receptor-ligand complex, receptor and ligand, respectively. The sum of molecular mechanics interaction energy ( E MM ), solvation energy ( G sol ) and the change of the conformational entropy at temperature T (-T S) are equal to the sum of G R+L , G R , and G L . In Eq. (9), E MM is given as the sum of intermolecular interaction energy ( E int ), van der Waals energy ( E vdW ), and electrostatic energy ( E elec ). The solvation free energy can be divided into polar ( G GB ) and nonpolar ( G SA ) and contributions (Eq. (10)). Here, for each complex, 1,000 structures extracted from the classical MD simulations with a simulation time between 600 and 800 ns were utilized to conduct binding free energy calculations and per-residue energy decomposition. The E int was canceled as the single trajectory strategy was executed. Based on parameters by Onufriev et al., the G GB was determined using a modified GB model (GB OBC1 ) (Onufriev, Bashford & David, 2000). The G SA was determined using the solvent accessible surface area (SASA) model ( G SA = σ *SASA). The parameter σ was set to 0.0072 kcal mol −1 Å −2 . The −T S was excluded from consideration because of relatively low prediction accuracy and high computational demand (Hou et al., 2011).

Free energy landscape (FEL) calculation based on aMD
The cumulant expansion to the second order method was employed to determine the FEL for each simulated complex from aMD simulations, because this method offers a good approximation for calculating the reweighting factor. In this study, the V(r) combined with PC1 and PC2 from PCA of the aMD simulation trajectories were utilized to recover the FEL (Miao et al., 2014b;Roe & Cheatham 3rd, 2013).

Evaluation of the stability of simulated complexes from classical MD simulations
To validate the docking results of the modeled complex of PARP-2 bound to Niraparib, structural alignments of the PARP-2/NMS-P118 and modeled PARP-2/Niraparib were performed with the corresponding crystal structures of PARP-1/NMS-P118 and PARP-1/Niraparib. As shown in Fig. 2A, alignment of the crystal structures between PARP-1/NMS-P118 (PDB code: 5A00) and PARP-2/NMS-P118 (PDB code: 4ZZY) exhibited relatively high similarity, with a RMSD of 0.704 Åfor heavy atoms. Similarly, the alignment of the modeled complex of PARP-2/Niraparib with the crystal structure of PARP-1/Niraparib (PDB code: 4R6E) mostly showed similarities with only minor differences, with a RMSD of 0.767 Å for heavy atoms (Fig. 2B). These results suggest that the predicted model is sufficient to study the dynamic features through further MD simulations. Firstly, 800 ns classical MD simulations for the modeled structure and the three crystal structures were performed. As a prerequisite for all further analyses, the dynamic stability of the simulated complexes was monitored by studying the RMSDs for all protein backbones (C α ) atoms and all ligand heavy atoms for each complex with the starting structure as a function of simulation time. Theoretically, smaller fluctuations of RMSDs indicate greater stability of the complex. As shown in Fig. 3, the time evolution of the RMSD values of C α atoms and ligand heavy atoms in each complex tend to converge after ∼100-300 ns simulations. During the simulation of the last 400 ns, the RMSD curves of both PARP-1 C α atoms displayed minor fluctuations (<1 Å), indicating that NMS-P118 constrained the protein structural flexibility of PARP-1. In contrast, those for both PARP-2 and NMS-P118 exhibited greater fluctuations compared to the complex of PARP-2/NMS-P118 (Fig. 3B). This attested to the highly unstable nature of PARP-2 when it was bound to selective PARP-1 inhibitor NMS-P118. In comparison, the RMSD curves for Niraparib in both PARP-1 and PARP-2 oscillated with minute fluctuations (<1 Å) during the last 400 ns simulation, indicating that the non-selective drug Niraparib constrained the structural flexibility of both PARP-1 and PARP-2. Based on these results, the structural and energetic properties for each complex were further analyzed in detail.

Dynamic features of each complex from classical MD simulations
The different flexibility for NMS-P118 and Niraparib in different proteins may result in different protein conformational transitions. To characterize the protein conformational transitions over time, PCA was performed to identify the trend in large-scale collective motions. Generally, the eigenvectors with the highest eigenvalues capture the majority of variance in the original protein conformational distributions (Mashiko, 2018;Sittel, Filk & Stock, 2017). Herein, the protein conformational ensembles were investigated by projecting the first two principal components (PC1 and PC2) from classical MD simulations onto a two-dimensional space. Generally, when PC1 and PC2 are plotted against each other, structures with high similarity are clustered together. As a result, a cluster represents a different state of protein conformation. As shown in Fig. 4, the protein conformational distributions for each complex were dynamic during 800 ns classical MD simulations and eventually reached overall stability. It is apparent that the conformational distributions of PARP-1/NMS-P118 were remarkably different from those of PARP-2/NMS-P118, while those for PARP-1/Niraparib and PARP-2/Niraparib were also different. Analysis of the protein conformational distributions clearly showed that PARP-2/NMS-P118 samples a wider conformational space compared to PARP-1/NMS-P118. However, those for PARP-1/Niraparib and PARP-2/Niraparib shared a certain degree of similarity. These results demonstrate that the selective PARP-1 inhibitor NMS-P118 bound to PARP-1 had a different protein conformational flexibility compared to PARP-1, not the non-selective drug Niraparib. To further investigate the correlation between the motion of the residues in the proteins, DCC analysis was performed. As plotted in Fig. 5, the direction of correlation is represented by a color gradient ranging from blue (negative correlation) to red (positive correlation). The correlation coefficient (−1 to +1), corresponding to three different colors: dark blue (−0.25 to −1) represents anti-correlation; dark red (0.25 to 1) represents positive correlation; blue represents anti-correlation (−0.25 to −1); and light red or light blue (−0.25 to +0.25) represents weak or no-correlation. It can be observed that both the red and blue regions in PARP-2/NMS-P118 were larger and more intense than those in PARP-1/NMS-P118, implying elevated correlation or anti-correlation motions in PARP-2/NMS-P118 (Figs. 5A and 5B). In comparison, color regions for PARP-1 and PARP-2 bound to the non-selective drug Niraparib were quite similar (Figs. 5C and 5D). These results indicate that the relative motions of different protein subdomains may correlate with different protein flexibility, which was responsible for drug selectivity.
To further highlight the key sub-domains, the C α of RMSF analyses were conducted. Low RMSF values of residues represented less flexibility, whereas high RMSF values indicated greater fluctuations in relation to their average position during simulation. As shown in Fig. 6, the RMSFs of PARP-1/NMS-P118 showed a high degree of similarity with those of PARP-1/Niraparib. However, the helix αF connected with loop 347-356 in PARP-2 exhibited amplified fluctuations when it binds to the selective PARP-1 inhibitor NMS-P118 compared with the non-selective drug Niraparib. The above results indicate that the binding of selective PARP-1 inhibitor NMS-P118 leads to increased fluctuations
To gain further insights into the vital residues in drug selectivity, per-residue decomposition based on the MM/GBSA method was employed to assess residue contributions to the binding of the protein-ligand complexes. Per-residue energy differences between Niraparib and NMS-P118 systems ( G = G Niraparib -G NMS-P118 ) were plotted to identify the key residues. Negative values represent the residues of Niraparib formed stronger interactions with the protein than those of NMS-P118, whereas positive values indicated quite the opposite, namely, that the residues of Niraparib formed weaker interactions with the protein than those of NMS-P118. As shown in Fig. 7A, the differences Table 1 Binding free energies of NMS-P118 and Niraparib in PARP-1 and PARP-2 (kcal/mol). Notes. E vdW , Van der Waals energy; E ele , Electrostatic energy; G GB , Electrostatic contribution to solvation; G SA , Non-polar contribution to solvation; E nonpolar , Non-polar interaction; E polar , polar interaction; G bind , Binding free energy.

NMS-P118 Niraparib
of G between PARP-1/NMS-P118 and PARP-1/Niraparib were quite small with G less than 0.5 kcal/mol. Alignment of the representative structures of PARP-1/NMS-P118 and PARP-1/Niraparib were highly similar (Fig. 7B). However, the residues of Ser-328,  formed significantly stronger interactions with Niraparib than with NMS-P118 in PARP-2 (Fig. 7C). Notably, the key residues of Ser-328, Gln-322 and Glu-335 were located in the helix αF, which exhibited amplified fluctuations in RMSF analysis. This may be due to the fact that the distinctive flexibility of the helix αF in PARP-2 bound to selective PARP-1 inhibitors induced energetic redistributions.

Evaluation of the stability of the simulated complexes from aMD simulations
To explore the conformational behaviors in greater detail, aMD simulations were first employed. Thereafter, the RMSDs of C α atoms and ligand heavy atoms were monitored. As shown in Fig. 8, the RMSD values of C α atoms and ligand heavy atoms for each complex reached equilibrium after 120-200 ns aMD simulations, suggesting that simulated complexes became dynamically stable through 800 ns aMD simulations. Interestingly, the fluctuations of selective PARP-1 inhibitor NMS-P118 bound to PARP-1 were much smaller than when bound to PARP-2. In contrast, the fluctuations of non-selective drug Niraparib were similar in both PARP-1 and PARP-2. A comparison of these results revealed that the selective PARP-1 inhibitor allowed for larger protein conformational changes of PARP-2. These findings were corroborated by the DCC analysis, which further revealed that the red and blue regions in PARP-2/NMS-P118 were both larger and more intense than those in PARP-1/NMS-P118, while those for PARP-1 and PARP-2 bound to Niraparib were quite similar (Fig. 9).
Based on the above findings, PCA was used to identify the various protein conformations obtained during the aMD simulations. As shown in Fig. 10, the protein conformations of for each complex was characterized by dynamic fluctuations during 800 ns aMD simulations. Similar with the results from classical MD simulaitons, the conformational The FEL was employed to further demonstrate the relationship between the changes of conformation and energy (Fig. 11). Generally, more energy wells (dark blue regions) represent greater conformational changes of the protein during aMD simulation (Han et al., 2019;Miao, Nichols & McCammon, 2014a). As shown in Fig. 11A, only a deep energy well for NMS-P118 bound to PARP-1 was observed throughout the whole 800 ns aMD simulation. In contrast, two major deep energy wells with a much wider range of distributions were observed for the complex of PARP-1/NMS-P118 (Fig. 11B). Nevertheless, those for Niraparib bound to PARP-1 and PARP-2 showed a similar range of distributions and both were confined to a single deep energy. These results highlight the unstable nature of PARP-2 bound to selective PARP-1 inhibitor NMS-P118.

DISCUSSION
PARP-1 inhibitors have been widely studied as potential cancer therapeutics for breast and ovarian cancers (Min & Im, 2020). The reported clinical candidates and preclinical PARP-1 inhibitors were designed with the purpose of imitating the nicotinamide portion of nicotinamide adenine dinucleotide (NAD + ), with which they compete for the corresponding PARP-1 binding site. A number of recent studies have indicated that due to their high sequence and structural similarity (Figs. 1A-1C), non-selective PARPs drugs may result in potentially undesirable side effects, especially between PARP-1 and PARP-2 (Eltze et al., 2008;Fatima et al., 2014;Papeo et al., 2015). In this aspect, the ideal drug candidate should be a highly selective PARP-1 inhibitor with greater subtype specificity. Until now, only a few highly selective PARP-1 inhibitors have been reported. However, even fewer computational modeling researches have been conducted to elucidate their underlying selective mechanisms. In the present study, a comprehensive molecular computational method was employed to demonstrate the selective mechanisms via two representative inhibitors (NMS-P118, Niraparib) with different selectivity to PARP-1 and PARP-2 (Figs. 1D-1E). Initially, molecular docking was applied to predict the complex of PARP-2/Niraparib. To examine the prediction accuracy, the predicted and the crystal structures were aligned. The RMSDs between the predict pose (PARP-2/Niraparib) and the crystal structure (PARP-1/Niraparib) were highly similar, a finding that was consistent with the results observed by alignment of crystal structures of NMS-P118 bound to PARP-1 and PARP-2 (Fig. 2). Based on the above findings, these structures were employed to further explore the dynamic behavior via classical MD simulations and aMD simulations. Classical MD simulations in conjunction with the RMSD, PCA, and DCC analyses provided compelling evidence that the conformation fluctuations were different for PARP-2 bound to selective PARP-1 inhibitors and PARP-2 bound to non-selective inhibitor . Further RMSF analysis revealed that the major structural variations were the conformational changes of the helix αF, which may account for drug selectivity (Fig. 6). According to the binding free energy calculations, the polar contributions had an obvious impact on drug selectivity (Table 1). Per-residue decomposition analysis further revealed that drug selectivity was primarily controlled by the residues of Ser-328, Gln-322, Glu-335, and Tyr-455, most of which are located in the helix αF (Fig. 7). As there are possible energy barriers between various meta-stable states, the classical MD simulations may still not be sufficient to sample the possible conformations (Hamelberg, Mongan & McCammon, 2004;Miao, Nichols & McCammon, 2014a). Therefore, an enhanced sampling technique that can sample the protein conformation at various metastable states is still required. Most of enhanced sampling techniques often require expert knowledge of the studied complexes, as these techniques require reaction coordinates. This limitation is overcome by the aMD simulations, which can explore the conformational behavior of the biomacromolecule without this requirement (Hamelberg, Mongan & McCammon, 2004;Miao, Nichols & McCammon, 2014a). Therefore, within the framework of this study, aMD simulation was employed to further sample the possible conformational ensembles. The results of RMSD, PCA, DCC analyses from aMD simulations indicated that the selective PAPR-1 inhibitor NMS-P118 may significantly disrupt the stability of the PARP-2 protein, not the non-selective drug Niraparib, a finding that was consistent with the results from classical MD simulations (Figs. 8-10). The FEL analysis further demonstrated the unstable nature of PAPR-2 bound to the selective PAPR-1 inhibitors NMS-P118 (Fig. 11). The preferential binding of NMS-P118 to PARP-1 takes precedence over PARP-2, resulting in a better inhibitory effect of PARP-1 than PARP-2, and greater PARP-1 selectivity. In summary, more selective PAPR-1 inhibitors may be required to evaluate the above findings via molecular modeling, which might help facilitate the rational design of high selective PAPR-1 inhibitors.