Computational Tool to Design Small Synthetic Inhibitors Selective for XIAP-BIR3 Domain

X-linked inhibitor of apoptosis protein (XIAP) exercises its biological function by locking up and inhibiting essential caspase-3, -7 and -9 toward apoptosis execution. It is overexpressed in multiple human cancers, and it plays an important role in cancer cells’ death skipping. Inhibition of XIAP-BIR3 domain and caspase-9 interaction was raised as a promising strategy to restore apoptosis in malignancy treatment. However, XIAP-BIR3 antagonists also inhibit cIAP1-2 BIR3 domains, leading to serious side effects. In this study, we worked on a theoretical model that allowed us to design and optimize selective synthetic XIAP-BIR3 antagonists. Firstly, we assessed various MM-PBSA strategies to predict the XIAP-BIR3 binding affinities of synthetic ligands. Molecular dynamics simulations using hydrogen mass repartition as an additional parametrization with and without entropic term computed by the interaction entropy approach produced the best correlations. These simulations were then exploited to generate 3D pharmacophores. Following an optimization with a training dataset, five features were enough to model XIAP-BIR3 synthetic ligands binding to two hydrogen bond donors, one hydrogen bond acceptor and two hydrophobic groups. The correlation between pharmacophoric features and computed MM-PBSA free energy revealed nine residues as crucial for synthetic ligand binding: Thr308, Glu314, Trp323, Leu307, Asp309, Trp310, Gly306, Gln319 and Lys297. Ultimately, and three of them seemed interesting to use to improve XIAP-BR3 versus cIAP-BIR3 selectivity: Lys297, Thr308 and Asp309.


Introduction
Apoptosis escape is one of the major causes of cancer development and progression. It also contributes to chemoresistance [1]. Restoring apoptosis in cancer cells is therefore a promising strategy for the new anticancer therapy development. X-linked inhibitor of apoptosis protein (XIAP), also known as "inhibitor of apoptosis protein 3" (IAP3) or "baculoviral IAP repeat-containing protein 4" (BIRC4), is a protein inhibiting cell apoptosis. To stop apoptotic cell death, XIAP binds to caspase-9, caspase-3 and caspase-7, three enzymes that are essential for apoptosis initiation and execution [2,3]. Binding to caspases allows XIAP to inhibit their activation to block intrinsic, as well as extrinsic, signaling apoptosis pathways. XIAP is overexpressed in many tumors [4], and its anti-apoptotic effect contributes to the cancer cell's escape from apoptosis. Moreover, it plays an important role in chemoresistance and has therefore become an important target for the malignancy treatment [4].
XIAP belongs to the inhibitor of apoptosis proteins (IAPs) family, which consists of eight different proteins. They all share the zinc-binding baculovirus-IAP-repeat (BIR) domain comprising around 70 amino acids. Among IAP family members, XIAP is the only antagonist binding affinities, using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach. MM-PBSA is widely employed to predict ligand-receptor free energy thanks to its efficiency and reasonable accuracy. Then, we exploited the molecular dynamics simulations that were producing the best agreement with experimental binding data in order to build a 3D pharmacophore for a synthetic XIAP-BIR3 antagonist. The pharmacophore performance was assessed and optimized with a testing chemolibrary. The MM-PBSA data analysis correlated with the held features of our final optimized pharmacophore. This also allowed us to highlight the residues to target and increase XIAP-BIR3 vs. cIAP1/2-BIR3 selectivity.

MM-PBSA Binding Free-Energy Prediction
To reach our goal, we started by setting up a method for rapid and reliable predictions of synthetic ligand-binding affinities to the XIAP-BIR3 domain. To assess the predictive performance, we selected four XIAP-BIR3 synthetic ligands for which the X-ray structures of their complexes with XIAP-BIR3 (residues 248-352) were known, and the binding constants were measured experimentally. At the time of our work, 36 3D structures with XIAP-BIR3 were available in the PDB database [19]; 4 NMR solution structures and 32 structures were solved using X-ray diffraction. Among them, only 16 were co-crystallized with a non-peptide or non-peptidomimetic ligand. For our study, we selected four complexes from sixteen available ones with the following PDB IDs: 5C7C [20], 5M6M [21], 5OQW [22] and 5M6L [21]. The choice of these complexes was based on the ligand structural diversity/similarity, as well as on the difference/similarity in their binding activities to the XIAP-BIR3 domain determined by fluorescence polarization assay (see Figure 1). simulation strategies. Our strategies included a ranking of relative XIAP-BI3 antagonist binding affinities, using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach. MM-PBSA is widely employed to predict ligand-receptor free energy thanks to its efficiency and reasonable accuracy. Then, we exploited the molecular dynamics simulations that were producing the best agreement with experimental binding data in order to build a 3D pharmacophore for a synthetic XIAP-BIR3 antagonist. The pharmacophore performance was assessed and optimized with a testing chemolibrary. The MM-PBSA data analysis correlated with the held features of our final optimized pharmacophore. This also allowed us to highlight the residues to target and increase XIAP-BIR3 vs cIAP1/2-BIR3 selectivity.

MM-PBSA Binding Free-Energy Prediction
To reach our goal, we started by setting up a method for rapid and reliable predictions of synthetic ligand-binding affinities to the XIAP-BIR3 domain. To assess the predictive performance, we selected four XIAP-BIR3 synthetic ligands for which the X-ray structures of their complexes with XIAP-BIR3 (residues 248-352) were known, and the binding constants were measured experimentally. At the time of our work, 36 3D structures with XIAP-BIR3 were available in the PDB database [19]; 4 NMR solution structures and 32 structures were solved using X-ray diffraction. Among them, only 16 were co-crystallized with a non-peptide or non-peptidomimetic ligand. For our study, we selected four complexes from sixteen available ones with the following PDB IDs: 5C7C [20], 5M6M [21] 5OQW [22] and 5M6L [21]. The choice of these complexes was based on the ligand structural diversity/similarity, as well as on the difference/similarity in their binding activities to the XIAP-BIR3 domain determined by fluorescence polarization assay (see Figure 1).   To predict the experimental ligand binding affinities, we chose to use the MM-PBSA method. Firstly, short molecular dynamics (MDs) simulations of 50 ns (NAMD 2.13) using all-atom CHARMM36m force field for protein and CGENFF for ligands were carried out for each tested XIAP-BIR3/ligand complex. During these simulations, the possible impact of different additional parameterizations, as the repartitioning of the hydrogen mass (HMR) [23] and the effect of parametrization of cation π-interactions for tryptophan, tyrosine and phenylalanine residues (WYF) [24], was tested. In HMR parametrization, the heavy atom mass is redistributed on the bound hydrogen atoms, making it possible to slow down the high-frequency movements of the macromolecule and increase the simulation time step. The cation-π interaction is a noncovalent binding force occurring between the aromatic ring, providing both a negative electrostatic potential surface and cations through electrostatic interaction. It was highlighted that cation-π interactions can play an important role in protein-ligand recognition. That is why we wanted to evaluate the influence of these two additional parametrizations in our study. In this way, each complex was simulated four times: (i) without any additional parametrization, (ii) with HMR, (iii) with WYF and (iv) with both WYF and HMR additional parametrizations. Finally, MM-PBSA methodology was applied on each simulation to predict ligand binding free energy ∆G MM−PBSA .
First, protein and ligand stability in the XIAP-BIR3 binding site in all simulations was checked visually and then through the Root Mean Square Deviation (RMSD) calculations. Visually, we observed that, during some simulations, the XIAP-BIR3 N and C termini fluctuated a lot (see Supplementary Materials Figure S1). Thus, to assess the stability of the binding site in each simulation, we aligned the protein conformations excluding both termini (without residues 248-258 for N-terminus and without residues 336-352 for C-terminus) and counted the backbone atom RMSD of the protein core (residues 259-335-see Supplementary Materials Figure S2). The RMSD did not exceed 1Å average in all simulations, thus confirming the XIAP-BIR3 binding site's stability. To check the ligand stability during each simulation, a protein alignment without both termini was applied at first, and then the ligand RMSDs were calculated. On average, the ligand RMSDs were smaller than 2Å in all simulations (see Supplementary Materials Figure S3), thus confirming the ligand stability in the simulated complexes (with only the exception of ligand in complex 5C7C simulated with WYF parametrization-see Supplementary Materials Figure S5). During this simulation, the 1-methyl-3-dimethyl-6-chloroindoline moiety of the ligand moves from its initial position around 15 ns, and it is positioned rather outside the XIAP-BIR3 binding site at the end of the simulation (RMSD = 3.97 ± 1.81Å).
Then, MM-PBSA analyses were carried out on the generated MD trajectories. Firstly, we estimated ∆G MM−PBSA values without taking into account the entropic term, −T∆S, as usually applied in the literature (Table 1). Among the four parametrizations, the simulation without any additional parametrization usually applied in theoretical studies has given rather poor correlations (r = 0.56884 see Figure 2). On the other hand, the simulations with HMR and HMR + WYF additional parametrizations made it possible to predict ligand affinities in the correct order and with very good correlation compared to the experimental data. The best correlation between experimental and predicted ∆G was observed for HMR parametrization (r ≈ 0.98 see Figure 2A). Even if the correlation with HMR additional parametrization was very satisfactory, predicted ∆G MM−PBSA values showed a significant difference (of about 30 kcal/mol) compared to the measured experimental data. We therefore decided to evaluate the contribution of the entropic part −T∆S, and we tested two methods described in the literature: the normal mode (NM) [25] method and the interaction entropy method (IE) [26,27]. The calculated ∆G MM−PBSA values, including the entropic contribution for each tested simulation protocol, are summarized in Table 1. We observed that the adding of the entropic member calculated using the NM decreased the prediction accuracy in all four protocols. The non-tested protocol with NM entropic term predicted the correct order of ligand affinities. For all of the tested protocols, we obtained a satisfactory correlation (r ≤ 0.81, Figure 2B). In contrast, the ∆G MM−PBSA values, which were calculated using the IE method, gave overall better results than those obtained using NM. A new MD simulation with the HMR additional parametrization produced a very good correlation (r = 0.97436; Figure 2C). In conclusion, the HMR additional parametrization produces the best correlation between predicted and experimental free energies when the entropic term is not taken into account ( Figure 2). Moreover, this parametrization also maintains a very good performance when the entropic term calculated using the IE method is included (r = 0.98264 without entropic term compared to r = 0.97436 with entropy from IE). Moreover, the predicted ∆G MM−PBSA values with the entropic term calculated using the IE method are closer to the experimental ∆G exp (Table 1); the differences between the predicted and experimental data were decreased to 10 kcal/mol.
We also probed the influence of the simulation temps on the correctness of the MM-PBSA prediction by computing ∆G MM−PBSA values at 10 ns and 25 ns of the simulations (see Supplementary Materials Table S1). We were able to identify certain trends. In the prediction without the entropic term, the lengthening of the simulations decreased the prediction accuracies for the simulations with any additional parametrization and with WYF additional parametrization, but it increased the exactness in the simulations with HMR additional parametrization. The same trend was also detected for the predictions with the entropic temp calculated using the IE method, except for the simulations without additional parametrization, for which the lengthening of the simulations increased the prediction accuracies this time. Nevertheless, in the case of predictions with the entropic term computed using NM, we did not succeed in deriving rules.
Next, we tested the prediction performance protocols on a ligand not belonging to the synthetic ligand category. We chose to use the AVPI tetrapeptide of Smac. To be able to compare the predicted free energies for the synthetic ligands to that of AVPI, the CGENFF force field was also applied to parameterize the AVPI tetrapeptide. As the XIAP-BIR3 3D structure co-crystallized with only AVPI tetrapeptide is not available in the PDB databank, we generated this complex by using the docking (see Supplementary Materials Figure S4) and submitted it to four MD simulation protocols, like for the previous complexes with synthetic ligands. The prediction accuracy of the binding free energy for the AVPI complex was also successful in the simulations with HMR additional parametrization, without the entropic term, as well as with the IE entropic term (Table 1). Very good results were also observed for the HMR + WYF additional parametrization protocol without the entropic term. This confirms that the use of the HMR additional parametrization significantly increases the exactness of the binding free-energy prediction in the MM-PBSA strategy. The inclusion of the entropic term calculated using the IE method maintains the correct accuracy for the HMR additional parametrization protocol, generating predicted values that are closer to the experimental ones.
To annotate the key interactions established between the studied ligands and the XIAP-BIR3 domain, contact maps were calculated for all simulated complexes, as well as their dynamical evolution along MD trajectories (see Supplementary Materials Figures S5 and S7). Afterward, we chose to focus our analysis on the 5M6L complex. Indeed, the 5M6L ligand is the most affine one, and we analyzed the contribution of each residue of XIAP-BIR3 to the ligand binding. To do so, an average ligand-binding free energy per residue was calculated from the MD simulation with HMR additional parametrization ( Figure 3). The polar and non-polar contribution to binding free energy par residue was then analyzed separately ( Figure 3C,D).
Molecules 2023, 28, 5155 6 of 19 Figure S4) and submitted it to four MD simulation protocols, like for the previous complexes with synthetic ligands. The prediction accuracy of the binding free energy for the AVPI complex was also successful in the simulations with HMR additional parametrization, without the entropic term, as well as with the IE entropic term (Table 1). Very good results were also observed for the HMR + WYF additional parametrization protocol without the entropic term. This confirms that the use of the HMR additional parametrization significantly increases the exactness of the binding free-energy prediction in the MM-PBSA strategy. The inclusion of the entropic term calculated using the IE method maintains the correct accuracy for the HMR additional parametrization protocol, generating predicted values that are closer to the experimental ones. To annotate the key interactions established between the studied ligands and the XIAP-BIR3 domain, contact maps were calculated for all simulated complexes, as well as their dynamical evolution along MD trajectories (see Supplementary Materials Figures S5 Figure 2. Experimental free-energy correlations with predicted ones using MM-PBSA approaches (A) without entropic term, (B) with entropic term calculated using the normal mode strategy and (C) with entropic term calculated using interaction entropy method (gray curves-without additional parametrization; red-HMR; blue-WYF; green-HMR and WYF additional parametrizations). and S7). Afterward, we chose to focus our analysis on the 5M6L complex. Indeed, the 5M6L ligand is the most affine one, and we analyzed the contribution of each residue of XIAP-BIR3 to the ligand binding. To do so, an average ligand-binding free energy per residue was calculated from the MD simulation with HMR additional parametrization ( Figure 3). The polar and non-polar contribution to binding free energy par residue was then analyzed separately ( Figure 3C,D). The analysis results showed that seven residues contribute to the ligand binding in the XIAP-BIR3 active site with a |ΔG| > 3.0 kcal/mol: Thr308, Glu314, Trp323, Leu307, Asp309, Trp310 and Gly306. Among these residues, the major polar contribution, calculated as the sum of electrostatic interaction energy and the polar solvation term, was assigned to the Glu314 (−8.0 ± 8.9 kcal/mol) and Thr308 (−5.6 ± 2.1 kcal/mol). These two residues, indeed, engage hydrogen bonds with the ligand: Glu314 side chain interacts through a salt bridge with one hydrogen atom of the protonated nitrogen of the piperazine ring of the ligand, and Thr308 interacts through its backbone carbonyl group with the second hydrogen atom. Thr308 established a second hydrogen bond by its backbone nitrogen atom with the ligand-carbonyl group. The occupation of these three hydrogen bonds at the course of the trajectory is globally in the same order, 87.4% for Glu314, 89.0% for the H-bond with Thr308 backbone nitrogen and 84.1% for the H-bond with Thr308 The analysis results showed that seven residues contribute to the ligand binding in the XIAP-BIR3 active site with a |∆G| > 3.0 kcal/mol: Thr308, Glu314, Trp323, Leu307, Asp309, Trp310 and Gly306. Among these residues, the major polar contribution, calculated as the sum of electrostatic interaction energy and the polar solvation term, was assigned to the Glu314 (−8.0 ± 8.9 kcal/mol) and Thr308 (−5.6 ± 2.1 kcal/mol). These two residues, indeed, engage hydrogen bonds with the ligand: Glu314 side chain interacts through a salt bridge with one hydrogen atom of the protonated nitrogen of the piperazine ring of the ligand, and Thr308 interacts through its backbone carbonyl group with the second hydrogen atom. Thr308 established a second hydrogen bond by its backbone nitrogen atom with the ligand-carbonyl group. The occupation of these three hydrogen bonds at the course of the trajectory is globally in the same order, 87.4% for Glu314, 89.0% for the H-bond with Thr308 backbone nitrogen and 84.1% for the H-bond with Thr308 backbone oxygen atom (see Supplementary Materials Figure S5). In the course of the dynamics, the fluctuations of the ligand position in the XIAP-BIR3 binding groove occasionally caused the ligand to approach Asp309 (−1.7 ± 6.4 kcal/mol). Indeed, in 28.9% of the trajectory, a hydrogen bond is also formed between the protonated nitrogen atom of the piperazine ring of the ligand and Asp309 main chain carbonyl group. The analysis of the non-polar contribution (sum of van der Waals interaction energy with the non-polar solvation term) revealed several residues as stabilizing the ligand binding through non-polar contacts, and the most important contributions were detected for Trp323 (−7.9 ± 0.8 kcal/mol), Leu307 (−6.6 ± 0.8 kcal/mol), Thr308 (−5.3 ± 1.0 kcal/mol), Gly306 (−3.3 ± 0.5 kcal/mol), Trp310 (−3.0 ± 0.5 kcal/mol) and Asp309 (−2.9 ± 1.1 kcal/mol). Interestingly, Thr308 and its neighboring residue Asp309 contribute to the ligand binding through polar and non-polar interactions.
The published structural and mutation data on the XIAP-BIR3 complex determined the following residues as being critical in binding the Smac N-terminus: Glu314, Gln319 and Trp323 [12]. The residues annotated in our study showed that the synthetic ligands interact through the same residues. However, the Gln319 contribution to the synthetic ligand binding is less important (−2.6 ± 1.3 kcal/mol) according to our study.

Construction and Validation of a 3D Pharmacophore of Synthetic XIAP-BIR3 Inhibitors
In the final step, we built a representative 3D pharmacophore for the synthetic ligands binding to the XIAP-BIR3 binding groove. Indeed, the pharmacophore approach has become one of the important ligand-based tools in in silico methodologies for drug research and development process. The interest of this approach is to highlight essential motifs for the ligand-target interaction, which can then be applied to guide the design of new ligands. The LigandScout software [28] used in this study allowed us to generate the 3D pharmacophore either using only the ligand structures (ligand-based approach) or using the target structure in addition to the ligand structure (structure-based pharmacophore) [29]. In our study, in order to take into account the essential structural motifs for the interaction with the target, we generated XIAP-BIR3 pharmacophores by using the structure-based approach. Our 3D pharmacophore model was built from 5M6L (X-ray structure co-crystallized with the most affine synthetic ligand to date). To generate the 3D pharmacophore, we exploited not only the crystallographic complex structure but also the previously generated MD trajectory simulated with HMR additional parametrization. We built one pharmacophore from the X-ray complex structure and another from the 2500 structures representing the entire MD trajectory (see Supplementary Materials Figure S6). The pharmacophore built from the X-ray structure had one hydrophobic and one aromatic-hydrophobic feature more and one hydrogen bond acceptor at least compared to the pharmacophore generated from MD trajectory (see Supplementary Materials Figure S6). The two pharmacophores were then aligned and merged, and the resulting pharmacophore comprised 8 features and 22 exclusion volumes (see Figure 4). Among the eight pharmacophoric features, Lig-andScout annotated four hydrophobic (H) and one aromatic-hydrophobic features (HAr), one hydrogen bond acceptor (HBA) and two hydrogen bond donors (HBD). The identified features of the pharmacophore reflect the ligand interaction with the XIAP-BIR3 binding site through Lys297, Thr308, Asp309, Glu314, Trp323 and Tyr324 (Figure 4). Thr308, Asp309, Glu314 and Trp323 residues already emerged as crucial from our free-energy decomposition per residue (Figure 3). In addition, LigandScout software also highlighted Lys297 and Tyr324 as important anchoring residues by positioning the hydrophobic features in their proximity. For Lys297 and Tyr324, we detected via our MM-PBSA analyses their dominant contribution non-polar one, which is in good agreement with the proposed hydrophobic characters of these residues. However, these two residues' contributions to the global ligand-binding free energy were rather moderate: Lys297 (−2.4 ± 4.1 kcal/mol) and Tyr324 (−1.7 ± 0.8 kcal/mol).
To assess the generated 3D pharmacophore accuracy to discriminate active and inactive compounds, we used it to screen our test chemical library (see Materials and Methods Table 2) composed of 173 active compounds and 5417 inactive compounds. A good pharmacophore model would be able to identify a significant portion of active compounds and as few inactive ones as possible. Thus, an efficient screening should have high sensitivity, specificity, enrichment factor and area under the ROC curve values. The screening on our testing library showed that the pharmacophore n • 1 performed rather weakly in terms of its sensitivity value (see Table 3). The sensitivity was only 2.9% when we overlaid all pharmacophoric features on the library ligands and 30.6% when we allowed the software to randomly discard one feature during alignment on the ligands. In view of these results, we decided to optimize it by omitting one or more features in order to test the impact on the overall screening performance. In view of the MM-PBSA results, we first tested the influence of features linked to Lys297 (first optimization round). We generated three pharmacophores, starting by deleting one (H4) of the two hydrophobic features linked to Lys297 (pharmacophore n • 2, see Figure 4). Next, we removed the second remaining hydrophobic feature (H3), and we conserved the hydrophobic-aromatic feature (HAr) (pharmacophore n • 3 on Figure 4). Finally, we removed the hydrophobic-aromatic feature (HAr) instead of the hydrophobic feature H3 (pharmacophore n • 4; see Figure 4).  To assess the generated 3D pharmacophore accuracy to discriminate active and inactive compounds, we used it to screen our test chemical library (see Materials and Methods Table 2) composed of 173 active compounds and 5417 inactive compounds. A good pharmacophore model would be able to identify a significant portion of active compounds and as few inactive ones as possible. Thus, an efficient screening should have high sensitivity, specificity, enrichment factor and area under the ROC curve values. The screening on our testing library showed that the pharmacophore n°1 performed rather weakly in terms of its sensitivity value (see Table 3). The sensitivity was only 2.9% when we overlaid all pharmacophoric features on the library ligands and 30.6% when we allowed the software to randomly discard one feature during alignment on the ligands. In view of these results, we decided to optimize it by omitting one or more features in order to test the impact on the overall screening performance. In view of the MM-PBSA results, we first tested the influence of features linked to Lys297 (first optimization round). We generated three pharmacophores, starting by deleting one (H4) of the two hydrophobic features linked to Lys297 (pharmacophore n°2, see Figure 4). Next, we removed the second remaining hydrophobic feature (H3), and we conserved the hydrophobic-aromatic feature (HAr) (pharmacophore n°3 on Figure 4). Finally, we removed the hydrophobic-aromatic feature (HAr) instead of the hydrophobic feature H3 (pharmacophore n°4; see Figure 4).  We observed that the removal of any of the three Lys297-related features did not notably improve the performance of the pharmacophore when the alignment of all pharmacophore features on the chemolibrary ligands was imposed: the sensitivity increased only by 1.1-6.9% (pharmacophore n • 2-n • 4; Table 3). In contrast, we noted that the sensitivity performances improved significantly for each tested pharmacophore when the screening allowed for the random exclusion of one feature during the alignment on ligands (pharmacophore n • 1-n • 4, Table 3). For example, the sensitivity increased from 9.8% to 77.5% for the pharmacophore n • 4. Among the four tested pharmacophores in the first round of optimization, the n • 4 one produced the best overall performance. We therefore selected it for the second optimization round. To perform this second optimization round, we analyzed, in detail, the discarded features in the screening (allowing for the random exclusion of one feature during the chemical-library alignment on pharmacophore n • 4). We thus identified that the most frequently omitted feature among the six initial ones was the H1 hydrophobic feature. In view of these results, we generated a pharmacophore n • 5 without the hydrophobic feature H1 (therefore composed of only five features). The library screening against this pharmacophore without the exclusion of any feature during alignment on the library revealed a very good performance (see Table 3). Altogether, the suppression of the H1 hydrophobic feature was unexpected because it was one of the features located near Trp323. Indeed, Trp323 is a residue highlighted as a one interacting strongly with the ligand in our per-residue free-energy analysis (non-polar contribution). To clarify the role of every part of the ligand in the interaction with the XIAP-BIR3 binding site, we carried out a decomposition of the interaction energy per ligand group (see Figure 5). The results confirmed the outcome of pharmacophore optimization. Indeed, we observed that Trp323 interacts strongly with 2,3-dihydro-pyrrolo [3,[2][3][4][5][6]pyridine scaffold (G1 on Figure 5) and with carbonyl-linker (G3, specially with its hydrophobic CH 2 group) and rather weakly with the morpholine part (G5 in Figure 5). Thus, the suppression of a feature placed on this morpholine will not affect the synthetic ligand interaction modelling with the XIAP-BIR3 binding site and, in particular, with Trp323. This analysis also highlighted that Tyr324 only plays a minor role in ligand binding (see Figure 5). The H2 presence of the hydrophobic feature of our pharmacophore is rather related to the interaction with Trp323 than with Tyr324.  with the XIAP-BIR3 binding site and, in particular, with Trp323. This analysis also highlighted that Tyr324 only plays a minor role in ligand binding (see Figure 5). The H2 presence of the hydrophobic feature of our pharmacophore is rather related to the interaction with Trp323 than with Tyr324.  In conclusion, this cross-analysis pointed out as important groups for ligand binding to XIAP-BIR3 the protonated piperazine moiety (G4) and carbonyl-linker (G3), with the two groups forming hydrogen bonds between ligand and XIAP-BIR3. We also highlighted the important contribution of the 2,3-dihydro-pyrrolo [3,[2][3][4][5][6]pyridine scaffold (G1) and 4-fluorobenzyl moiety (G2), which interact mainly in a hydrophobic manner. We observed that the G1 group interacts mostly with Trp323, while the G2 group interacts simultaneously with several residues on a moderate level: Lys297, Gly306, Leu307 and Thr308 (see Figure 5B). These results lead us to believe that this pocket of the BIR3 binding site hosting the group G2 is not completely filled by published synthetized ligands.
Interestingly, comparing XIAP-BIR3 residues to cIAP1-BIR3 and cIAP2-BIR3 residues revealed crucial residues that bind to synthetic ligands (see Figure 5C). Only three residues are not conserved between XIAP and cIAP1-2: Lys297/Asp, Thr308/Arg and Asp309/Cys. It can therefore be concluded that targeting these three residue side chains could increase the ligand selectivity for XIAP versus cIAP1-2. Among them, the Lys297 side-chain interaction through pharmacomodulations of the G2 group seems to be the most promising way to improve the selectivity.
A pharmacophore for XIAP-BIR3 inhibitors has already been published in the literature [30]. Opo and co-workers also generated their pharmacophore by using LigandScout software, but from the 5OQW X-ray structure, one of the structures used in our MM-PBSA study. Unlike us, their pharmacophore was generated by also considering the co-crystallized water molecules in the XIAP-BIR3 binding site, and it is composed of twelve features: four hydrophobic features, one positive ionizable group, three hydrogen bond acceptors, five hydrogen bond donors and fifteen exclusion volumes. The performance of our optimized pharmacophore n • 5 is of the same level compared to the AUC value, and it is even better from the point of view of the enrichment factor, 20.4, with respect to 10.0. Indeed, our simpler 3D pharmacophore, with only five features, represents the chemical groups necessary for synthetic ligands to bind in the XIAP-BIR3 binding groove.

Conclusions
The aim of this study was to establish a computational tool that would allow us to design and/or optimize the XIAP-BIR3 selective synthetic ligands. Firstly, we demonstrated that the MM-PBSA approach for free-energy calculation from molecular dynamics predicts ligand-binding affinities in the right order (CHARMM36 force field and HMR additional parametrization). The inclusion of the entropic term calculated using the interaction energy methodology maintains good predictability of the binding affinities, and the free-energy predicted values tend to be closer to the experimental ones. Our study showed that the key residues for synthetic ligand binding are Thr308 and Glu314. They both form hydrogen bonds with the ligand protonated nitrogen of the piperazine ring and with the ligand carbonyl-linker part. Therefore, the presence of double H-bond donors, as well as the Hbond acceptor on ligands, is required to bind to XIAP-BIR3. Indeed, the 3D pharmacophore for XIAP-BIR3 inhibitors, built and optimized by us, reflects these observations. Among the retained five features that are necessary for the binding of synthetic ligands to XIAP-BIR3, we revealed two hydrogen-bond donors and one hydrogen-bond acceptor. Nevertheless, these H-bonds are formed either through Thr308 main-chain atoms or through Glu314 side-chain carboxyl groups. Glu134 is a conserved residue (Glu/Asp) in the IAPs family and does not allow us to reach XIAP-BIR3 ligand selectivity.

Structure Preparation
During this study, X-ray structures of the following complexes with XIAP-BIR3 retrieved from the PDB database were used (Table 4): 5C7C [20], 5M6M [21], 5OQW [22] and 5M6L [21]. Then, to calculate the experimental value of ligand binding free energy ∆G exp , the following equation was applied: where the gas constant value was R = 1.985 8775 × 10 −3 kcal·K −1 ·mol −1 and T = 303, 15 K. As the XIAP-BIR3 complex with Smac AVPI tetra peptide was not solved to date, to prepare this complex, we chose to apply the docking strategy. The AVPI tetra peptide, which was retrieved from the PDB structure of XIAP-BIR3 co-crystallized with full SMAC (PDB ID: 1G73 [12]), was docked into the X-ray crystal structure 5C7C of XIAP-BIR3, without waters and bounded synthetic ligand. The docking was carried out using Glide software (Schrödinger ® ) with the default parameters [31]. The protein structure was prepared before using Schrödinger Protein Preparation Wizard (PrepWizard) [31]. The Glide score was used to evaluate the generated AVPI poses, and the pose with the best score was selected as a starting structure for MD simulations.

Molecular Dynamics Simulations
All dynamics simulations were performed using NAMD 2.13 [32]. The all-atom CHARMM36m forcefield [33,34] was used for the XIAP-BIR3 protein, and CGENFF [35] was used for ligands. Two additional force fields were applied in this order: (i) hydrogen mass repartition (HMR) [23] and (ii) additional parametrizations of π-cation interaction for tryptophan, tyrosine and phenylalanine (WYF) residues [24]. The starting systems were generated by the CHARMM-GUI server [36]. The three Cys residues, Cys300, Cys303 and Cys327, coordinating the zinc cation were parametrized as anionic cysteine (CYM), and the four histidine residues were modeled along the PropKa suggestion [37,38]: 320 and 346 as HSD, and 302 and 343 as HSE. Each system was solvated using the TIP3P explicit water model [39] within a rectangular box; the box size ensured that the simulated complex was at a minimum distance of 10 Å from the edge. To neutralize the total charge system, 0.15 M of KCl was added. The vacuum dielectric constant was used during all calculations. Cubic periodic boundary conditions were applied to the systems by using the IMAGE algorithm. The applied cutoff distance was 16 Å, and van der Waals interactions were truncated using a force-switching function between 10 and 12 Å. The Particle Mesh Ewald (PME) was used to calculate long-range electrostatic interactions [40]. The SHAKE algorithm was applied to restrain all bonds involving hydrogen atoms.
Each complex was subjected to 4 different dynamics simulations to assess the impact of additional parameters of the force field (HMR and WYF). Firstly, each underwent energy minimization in 10,000 steps, with harmonic restraints applied on heavy atoms (1 kcal/mol/Å 2 force constant for backbone atoms and 0.5 kcal/mol/Å 2 force constant for sidechain atoms), followed by 10,000 steps of minimization without any restraints. Next, the minimized systems were heated to 303.15 K, and the dynamics were temperatureequilibrated during 250 ps via heating reassignment under NVT conditions, with harmonic restraints applied to heavy atoms. Finally, the systems ran freely for 50 ns under NPT conditions, with a time step of 2 fs for the simulations without HMR additional parametrization and 4 fs with HMR parametrization. Langevin dynamics with a damping coefficient of 1 ps −1 was used to maintain the system temperature, and the Nosé-Hover-Langevin piston method was used to control the pressure at 1 atm. Production trajectories were saved every 20 ps, and subsequent analyses were performed using the CHARMM program version c40b2 [41].

MM-PBSA Prediction
To estimate the binding free energy of a ligand to a receptor, the MM-PBSA method was applied. MM-PBSA calculates the binding free energy as the sum of the classical enthalpic contributions (bound, van der Waals and electrostatic energies), the solvation-free energies and the entropic contribution [42]. In our protocol, we simulated only the complex and the ensemble of data for the free receptor and ligand for each snapshot created by removing the appropriate atoms. Then, the binding free energy was calculated using the following equation [43,44]: where the first term, E inte , is the standard Molecular Mechanics ligand-receptor interaction energy calculated from electrostatic and van der Waals interaction energies; and G polar and G non−polar are the polar and non-polar contributions to the solvation free energies, respectively. G polar was determined by solving the Poisson-Boltzmann equation, whereas the G non−polar term was estimated from a linear relation to the solvent accessible surface area (SASA), using the equation 0.00542 × SASA + 0.92. The value of G polar and G non−polar in each trajectory snapshot was calculated as follows: To calculate all contributions to the MM-PBSA CHARMM, homemade scripts were applied.
For the last term, entropic contribution, three approaches were tested: firstly, (i) the entropic term was omitted, and then the entropic term was approximated using either (ii) normal mode analysis of harmonic frequencies calculated at the Molecular Mechanics level [25] or using the recently published approach (iii) interaction entropy that investigates the entropy change upon binding [26,27]. The interaction entropy (IE) contribution to binding free energy is defined as follows: where ∆E inte = E inte − E inte is the fluctuation of the receptor-ligand interaction energy. The ensemble average of E inte was extracted from MD simulations.

3D Pharmacophore Building
The 3D pharmacophore was generated using the 3D structure of human XIAP-BIR3 co-crystallized with a synthetic inhibitor (PDB ID code 5M6L [21]) and using the structural data resulting from its MD simulations with additional HMR parametrization. LigandScout software [28] was applied for the detection and interpretation of crucial interaction patterns between XIAP-BIR3 and the ligand. LigandScout extracts and interprets ligands and their macromolecular environment from PDB files and automatically creates and visualizes advanced 3D pharmacophores. Besides the features representing the interaction of the ligand with the target (as hydrogen bond donor, hydrogen bond acceptor, hydrophobic group, etc.), we also generated our pharmacophore exclusion volumes.
In order to compare the pharmacophore performance during the optimization rounds, we referred to the evaluation metrics, such as sensitivity, specificity, enrichment factor and area under the ROCs (Receiver Operating Characteristics) curve. The purpose of these statistical tests is to clearly differentiate the pharmacophore performances from each other. Sensitivity determines the model's ability to retain active molecules. In contrast, specificity assesses the model's ability to reject inactive molecules. The enrichment factor (EF) is the ratio between the number of active ligands aligned with the pharmacophore and the total number of ligands. The enrichment factor therefore measures the enrichment of active compounds in the hit list compared to a pure random selection. A method that is superior to random selection results in an EF value greater than 1. The area under the ROC curve (AUC), on the other hand, evaluates the probability that, for two screened ligands, one active and the other inactive, the score value is higher for the active than for the inactive. The AUC value between 0.5 and 1 is mandatory for a model to be validated (AUC = 1 perfect prediction; AUC = 0.5 random prediction).
These tests are defined by the following equations [45]: where TP is the number of active ligands aligned with the pharmacophore, TN is the number of discarded inactive ligands and FP is the number of inactive ligands aligned with the pharmacophore. Moreover, Htot = (TP + FP) is the total number of aligned ligands, A is the number of active ligands in the initial dataset, I is the number of inactive ligands in the initial dataset and D is the total number of ligands in the initial dataset. In general, an efficient screening should have high values of sensitivity, specificity, enrichment factor and area under the ROC curve. During the screening, we applied an alignment of all features on the ligands and the extraction mode based on the most stable ligand energy conformation. For each pharmacophore, we carried out two screenings: (i) in the first one, we imposed that all features of the 3D pharmacophore model had to be aligned to ligands; and (ii) in the second screening, the software could randomly discard one feature during the alignment on ligands.

Testing Chemical Library
To test our 3D pharmacophore performance, we created a validation dataset. To build it, we first detected all compounds annotated with XIAP-BIR3 biological activities in CHEMBL database [46,47]. From this dataset, the redundant ligands and the peptide ligands or peptide-mimetic ones were deleted. This resulted in a dataset of 487 synthetic XIAP-BIR3 ligands. Then, it was necessary to determine a cutoff (threshold) to classify the ligands in terms of biological activity (active/inactive). As the number of ligands available for XIAP-BIR3 was limited, we applied the approach that is usually used in the literature for targets that are not perfectly annotated [48]. In this strategy, a pIC 50 (pIC 50 = −log 10 IC 50 ) value of 6 is usually used as a cutoff. All molecules with a pIC 50 value greater than 6 are considered active, and those with a value under 6 are considered inactive. Moreover, molecules with a pIC 50 value between 4.5 and 6 are excluded from the dataset (see Table 2). This led us to a dataset with a ratio of active to inactive ligands of 1/1.3 (i.e., 173 active ligands with respect to 218 inactive ligands). This ratio was not sufficient to obtain statistically significant values during the screening and correctly assess its performance.
Thus, an increase in the number of inactive compounds in our dataset was necessary to improve the screening performance, and to do so, we included XIAP-BIR3 decoys from the DUD-E database [48] (compounds were assumed to be inactive and built from fragments of active molecules). For our target, there were 5199 decoys available that we added to our starting dataset, resulting in a library of 5686 compounds. This library was then used to evaluate the performance of our 3D pharmacophores (see Table 2).
The 3D structure of library compounds was generated and protonated at pH = 7.4, using LigandPrep tool of Schrödinger ® software [31].
Supplementary Materials: The following supporting information can be downloaded at https:// www.mdpi.com/article/10.3390/molecules28135155/s1, Figure S1: Alignment of the final structure from molecular dynamics simulation on the initial one; Figure S2: RMSD of the protein backbone atoms along the MD simulations; Figure S3: The RMSD of the ligands along the MD simulations; Figure S4: (A) Selected pose from docking for AVPI into 5C7C structure. (B) The alignment of the end structure from molecular dynamics simulation on the initial one for four studied complexes with synthetic ligands; Figure S5: Key interacting residues of XIAP-BIR3 with ligands and their dynamical evolution along MD trajectories for all simulated complexes; Figure S6: Two pharmacophores used to create the final pharmacophore; Figure S7: ∆G MM-PBSA per residue for the five studied complexes; Table S1: Experimental binding free energy ∆G exp and calculated ones using the MM-PBSA method ∆G MM-PBSA for different parametrizations of the force field; Table S2: Pearson correlation coefficient between ∆G exp free energy and calculated one, using the MM-PBSA method ∆G MM-PBSA for different parametrizations of the force field and at different times of simulations.