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

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⁻¹ Å⁻² 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⁻¹ Å⁻² 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 $\Delta V\left(r\right)$ of the system is calculated according to the following Eq. (2):

(1)${\displaystyle \Delta V\left(r\right)=0\Delta V\left(r\right)\ge E}$ (2)${\displaystyle \Delta V\left(r\right)=\frac{{\left[E-\Delta V\left(r\right)\right]}^2}{\alpha +\left[E-\Delta V\left(r\right)\right]}\Delta V\left(r\right)<E}$

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 α. (3)${\displaystyle E_{total}=V_{total\text{\_}avg}+n\times {\mathrm{\alpha }}_{totalral}}$ (4)${\displaystyle {\mathrm{\alpha }}_{totalral}=0.2\times N_{atoms}\left(N=1,2,3\dots \right)}$ (5)${\displaystyle E_{dihedral}=V_{dihedral\text{\_}avg}+\left(3.5\times N_{res}\right)}$ (6)${\displaystyle {\mathrm{\alpha }}_{dihedral}=3.5\times \frac{N_{res}}{5}}$

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. (7)${\displaystyle C_{ij}=\frac{〈\Delta r_i\Delta r_j〉}{\sqrt{〈\Delta r_i^2\Delta r_j^2〉}}}$

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: (8)${\displaystyle \Delta G_{bind}=\Delta {\mathrm{G}}_{\mathrm{R}+\mathrm{L}}-\left(\Delta {\mathrm{G}}_{\mathrm{R}}+\Delta {\mathrm{G}}_{\mathrm{L}}\right)=\Delta {\mathrm{E}}_{\mathrm{MM}}+\Delta {\mathrm{G}}_{\mathrm{sol}}-\mathrm{T}\Delta \mathrm{S}}$ (9)${\displaystyle \Delta {\mathrm{E}}_{\mathrm{MM}}=\Delta {\mathrm{E}}_{\mathrm{int}}+\Delta {\mathrm{E}}_{\mathrm{vdW}}+\Delta {\mathrm{E}}_{\text{elec}}}$ (10)${\displaystyle \Delta {\mathrm{G}}_{\mathrm{sol}}=\Delta {\mathrm{G}}_{\mathrm{GB}}+\Delta {\mathrm{G}}_{\mathrm{SA}}}$

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⁻¹ Å⁻². 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 of PARP-2 (especially the αF and loop 347–356), which might be the dominating driving force for the redistributions of free energies.

Binding free energy calculations using the MM/GBSA method based on classical MD simulation trajectories

MM/GBSA binding free energy calculations based on classical MD simulation trajectories were applied for the assessment of the energy properties of NMS-P118 and Niraparib when they were bound with PARP-1 or PARP-2. As listed in Table 1, the predicted binding free energies (ΔG_binding) for PARP-1/NMS-P118, PARP-2/NMS-P118, PARP-1/Niraparib and PARP-2/Niraparib were −46.06 ± 0.15, −36.40 ± 0.14, −44.35 ± 0.22, and −43.89 ± 0.24 kcal/mol, respectively. It is clear that the ΔG_binding was highly correlated with the experimental data reported, and that for each system, different energy terms contribute differentially to the ΔG_binding. In this study, only the polar contributions (ΔE_elec + ΔG_GB) for the most vital factor were discussed. The polar contributions for the PARP-1/NMS-P118 and PARP-2/NMS-P118 were significantly different at 10.65 ± 0.34 and 19.29 ± 0.16 kcal/mol, respectively. In comparison, those for PARP-1/Niraparib and PARP-2/Niraparib were quite similar at 9.72 ± 0.29 and 9.84 ± 0.17 kcal/mol. The non-polar contributions (ΔE_vdW + ΔG_SA) for PARP-1/NMS-P118 and PARP-2/NMS-P118 (−56.70 ± 0.12 and −55.69 ± 0.19 kcal/mol) were almost identical with those for PARP-1/Niraparib and PARP-2/Niraparib (−53.92 ± 0.26 and −53.72 ± 0.41 kcal/mol). Taken together, these results demonstrate that the polar contribution has a significant impact on drug selectivity to PARP-1 and PARP-2.

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 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, Gln-322, Glu-335, and Tyr-455 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 distributions of PARP-1/NMS-P118 were remarkably different from those of PARP-2/NMS-P118 (Figs. 10A–10B). Meanwhile, those for PARP-1/Niraparib and PARP-2/Niraparib were also different (Figs. 10C–10D). Compared to the PARP-1/NMS-P118 complex, the PARP-2/NMS-P118 complex exhibited more structural clusters and a wider range of conformational distributions (Figs. 10A–10B). However, the conformational distributions for Niraparib bound to PARP-1 and PARP-2 had a somewhat similar range (Figs. 10C–10D), suggesting that the selective PARP-1 inhibitors could stabilize the protein conformation of PARP-1. These results were in agreement with the RMSDs, DCC and PCA analyses from the classical MD simulations.

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.