Discovery of Novel Tankyrase Inhibitors through Molecular Docking-Based Virtual Screening and Molecular Dynamics Simulation Studies

Tankyrase enzymes (TNKS), a core part of the canonical Wnt pathway, are a promising target in the search for potential anti-cancer agents. Although several hundreds of the TNKS inhibitors are currently known, identification of their novel chemotypes attracts considerable interest. In this study, the molecular docking and machine learning-based virtual screening techniques combined with the physico-chemical and ADMET (absorption, distribution, metabolism, excretion, toxicity) profile prediction and molecular dynamics simulations were applied to a subset of the ZINC database containing about 1.7 M commercially available compounds. Out of seven candidate compounds biologically evaluated in vitro for their inhibition of the TNKS2 enzyme using immunochemical assay, two compounds have shown a decent level of inhibitory activity with the IC50 values of less than 10 nM and 10 μM. Relatively simple scores based on molecular docking or MM-PBSA (molecular mechanics, Poisson-Boltzmann, surface area) methods proved unsuitable for predicting the effect of structural modification or for accurate ranking of the compounds based on their binding energies. On the other hand, the molecular dynamics simulations and Free Energy Perturbation (FEP) calculations allowed us to further decipher the structure-activity relationships and retrospectively analyze the docking-based virtual screening performance. This approach can be applied at the subsequent lead optimization stages.


Introduction
The tankyrase enzymes (TNKS1 and TNKS2, also known as PARP5a and PARP5b) belong to the poly(ADP-ribose) polymerase (PARP) superfamily. They play vital roles in mitosis control, telomere maintenance and regulation of the canonical Wnt pathway [1]. Since aberrant Wnt signaling

•
Semi-rigid molecular docking into the tankyrase binding site was performed using the Smina software and the compounds were selected based on the values of the Vinardo scoring function. The previously developed machine learning-based scoring function was also employed as an additional screening filter.

•
Compounds that have acceptable molecular weight, lipophilicity (LogP), aqueous solubility and human intestinal absorption as well as low risk of hERG-mediated cardiac toxicity were selected (the properties were predicted using previously developed QSPR/QSAR models). • Expert analysis of the resulting compounds was performed to eliminate potentially unstable, reactive or excessively complex structures.

•
For the seven selected compounds, molecular dynamics simulations and MM-PBSA calculations were carried out in order to provide additional independent assessment of their potential activity. • Biological evaluation of inhibitory activity of the selected compounds was carried out.
Even with steady improvement in the accuracy of computational methods over the years, it is not uncommon when only a fraction of the compounds predicted to be active shows some real activity. To minimize these risks, we used consensus scoring including molecular docking, ML scoring, QSAR models for the physico-chemical profile prediction and MM-PBSA method for binding energy estimation.
Although the MM-PBSA binding energy estimates show a broad range of correlations to the experimental values [18], they are widely used in practice and could, in our opinion, provide useful complement to the docking scores. In order to estimate the binding energies of tankyrase inhibitors, a preliminary molecular dynamics simulation of 30 ns was performed. The resulting system state was used as a starting point for ten independent runs of 5 ns each as suggested in the work [19]. The mean and confidence interval RMSD (root mean square deviation) values were estimated using the bootstrap procedure for each run and aggregated using mean and L2-norm, respectively.
The molecular docking and the closely related ML-based scoring served as primary screening filters reducing the initial library to the relatively small focused library of 174 compounds. It is worth noting that the distribution of docking scores for the screening library was close to normal with the mean value of −8.5 kcal/mol and the standard deviation of 1.7 kcal/mol. Then the QSAR/QSPR models were used to select 17 compounds for further expert assessment. Seven compounds selected by this virtual screening workflow are shown in Figure 1. These compounds were further evaluated in vitro against the tankyrase enzyme.
Molecules 2020, 25, x 3 of 15  Semi-rigid molecular docking into the tankyrase binding site was performed using the Smina software and the compounds were selected based on the values of the Vinardo scoring function. The previously developed machine learning-based scoring function was also employed as an additional screening filter.  Compounds that have acceptable molecular weight, lipophilicity (LogP), aqueous solubility and human intestinal absorption as well as low risk of hERG-mediated cardiac toxicity were selected (the properties were predicted using previously developed QSPR/QSAR models).  Expert analysis of the resulting compounds was performed to eliminate potentially unstable, reactive or excessively complex structures.  For the seven selected compounds, molecular dynamics simulations and MM-PBSA calculations were carried out in order to provide additional independent assessment of their potential activity.  Biological evaluation of inhibitory activity of the selected compounds was carried out.
Even with steady improvement in the accuracy of computational methods over the years, it is not uncommon when only a fraction of the compounds predicted to be active shows some real activity. To minimize these risks, we used consensus scoring including molecular docking, ML scoring, QSAR models for the physico-chemical profile prediction and MM-PBSA method for binding energy estimation.
Although the MM-PBSA binding energy estimates show a broad range of correlations to the experimental values [18], they are widely used in practice and could, in our opinion, provide useful complement to the docking scores. In order to estimate the binding energies of tankyrase inhibitors, a preliminary molecular dynamics simulation of 30 ns was performed. The resulting system state was used as a starting point for ten independent runs of 5 ns each as suggested in the work [19]. The mean and confidence interval RMSD (root mean square deviation) values were estimated using the bootstrap procedure for each run and aggregated using mean and L2-norm, respectively.
The molecular docking and the closely related ML-based scoring served as primary screening filters reducing the initial library to the relatively small focused library of 174 compounds. It is worth noting that the distribution of docking scores for the screening library was close to normal with the mean value of −8.5 kcal/mol and the standard deviation of 1.7 kcal/mol. Then the QSAR/QSPR models were used to select 17 compounds for further expert assessment. Seven compounds selected by this virtual screening workflow are shown in Figure 1. These compounds were further evaluated in vitro against the tankyrase enzyme.

Biological Evaluation
The inhibitory activity of the compounds was determined in vitro by measuring the tankyrase enzyme activity using immunochemical assay to detect the accumulation of poly(ADP-ribose) (PAR) in the course of the PARP enzymatic reaction. The initial screening results of the compounds A1-A7 at the concentration of 20 µM and NAD + at 1 µM are shown in Figure 2. It can be seen that PAR is absent only in two positions corresponding to the compound A1. In positions containing the compound A3, the product of the enzymatic reaction is present in a significantly smaller amount than in the absence of inhibition. These data suggest that compounds A1 and A3 likely act as inhibitors of the tankyrase enzyme. These two compounds based on similar scaffolds were selected for further evaluation.
The inhibitory activity of the compounds was determined in vitro by measuring the tankyrase enzyme activity using immunochemical assay to detect the accumulation of poly(ADP-ribose) (PAR) in the course of the PARP enzymatic reaction. The initial screening results of the compounds A1-A7 at the concentration of 20 μM and NAD + at 1 μM are shown in Figure 2. It can be seen that PAR is absent only in two positions corresponding to the compound A1. In positions containing the compound A3, the product of the enzymatic reaction is present in a significantly smaller amount than in the absence of inhibition. These data suggest that compounds A1 and A3 likely act as inhibitors of the tankyrase enzyme. These two compounds based on similar scaffolds were selected for further evaluation.
In order to measure their inhibitory activity, the concentration-response curves for the compounds A1 and A3 were obtained at 100 μM NAD + concentration (see Appendix A, Figure A1 and Table A1). Taking into account small number of repeated experiments and high data variance, the IC50 values can be cautiously estimated as less than 10 nM for compound A1 and less than 10 μM for compound A3. Moreover, subsequent preliminary experiments suggest that the inhibition for the compound A1 is not competitive and that inhibition of the enzyme by compounds A1 and A3 is independent of the duration of preincubation and is not associated with irreversible inhibition.

Retrospective Analysis of the Virtual Screening Results
Additional studies were performed in order to evaluate to what extent the methods used for the hit-oriented virtual screening could be employed during further lead optimization. In particular, the docking and ML model predictions were recalculated to collect better statistics and compared to the MM-PBSA results.
The molecular docking procedure as implemented in the Autodock Vina and Smina software is stochastic in nature and its different runs result in different ligand positioning in the binding site and, consequently, different scoring function values. To analyze the effect of this uncertainty, one hundred docking runs were performed for each of the A1-A7 compounds starting from various initial conformations and random seeds. The averaged resulting scores and standard deviations are shown in Table 1.
The docking poses obtained on the previous step were also scored using the previously developed target-oriented machine learning scoring approach [13]. The model employed in the present study was refined taking into account additional data available in the most recent version of In order to measure their inhibitory activity, the concentration-response curves for the compounds A1 and A3 were obtained at 100 µM NAD + concentration (see Appendix A, Figure A1 and Table A1). Taking into account small number of repeated experiments and high data variance, the IC 50 values can be cautiously estimated as less than 10 nM for compound A1 and less than 10 µM for compound A3. Moreover, subsequent preliminary experiments suggest that the inhibition for the compound A1 is not competitive and that inhibition of the enzyme by compounds A1 and A3 is independent of the duration of preincubation and is not associated with irreversible inhibition.

Retrospective Analysis of the Virtual Screening Results
Additional studies were performed in order to evaluate to what extent the methods used for the hit-oriented virtual screening could be employed during further lead optimization. In particular, the docking and ML model predictions were recalculated to collect better statistics and compared to the MM-PBSA results.
The molecular docking procedure as implemented in the Autodock Vina and Smina software is stochastic in nature and its different runs result in different ligand positioning in the binding site and, consequently, different scoring function values. To analyze the effect of this uncertainty, one hundred docking runs were performed for each of the A1-A7 compounds starting from various initial conformations and random seeds. The averaged resulting scores and standard deviations are shown in Table 1.
The docking poses obtained on the previous step were also scored using the previously developed target-oriented machine learning scoring approach [13]. The model employed in the present study was refined taking into account additional data available in the most recent version of the ChEMBL 26 database [20]. Similar to the built-in scoring functions, we collected the predicted activity probabilities and their standard deviation values (see Table 1).
As can be seen from Table 1, for all the compounds previously selected for the in vitro testing, the docking scores and MM-PBSA energies are quite similar and promising. The binding probabilities predicted by the ML model show a greater variance-the model predicts that the compounds A1, A2, A3, A6 and A7 could be binders with decent probability while the probabilities for the A4 and A5 compounds are relatively lower. However, taking into consideration the results from all three components of the consensus scoring approach, the data would suggest significant likelihood of binding for all the compounds.

Molecular Dynamics Studies
Looking retrospectively at the values presented in Table 1, one can see that almost all the predictions are quite similar although only two compounds have shown some level of inhibitory activity in the experiment. Obviously, such scoring techniques cannot be used to predict the effect of structural modification on the binding affinity or to rank the compounds based on their binding energies.
A predictive model capable of suggesting and evaluating such modifications would be highly useful and desirable in the subsequent lead optimization stages. Thus, the binding of the discovered active compounds was analyzed in more detail to get quantitative insight into the relations between the estimated binding energy and the inhibitory activity level. Then the binding energies were estimated via the free energy perturbation (FEP) method that is believed to be more accurate than molecular docking or MM-PBSA method.

Binding Modes
The binding modes of compounds A1 and A3 predicted by molecular docking and refined by molecular dynamics are shown in Figure 3. As can be seen, despite the scaffold similarity, A1 and A3 have been predicted to bind differently. Although localized roughly in the same pocket, the compounds take different conformations with apparent intramolecular π-π stacking in A3. Thus, compound A1 is expected to form hydrogen bonds with the backbone atoms of Gly1032, Asp1045, Tyr1060 and the side chain of Lys1067 and interact with His1048 through π-π stacking. The compound A3 could interact with Tyr1060 and Tyr1071 via π-π stacking and form hydrogen bonds with the backbone atoms of Gly1032, Tyr1050 and Tyr1060.

FEP Calculations
The alchemical free energy perturbation (FEP) method is based on a non-physical thermodynamic cycle where the binding free energy is computed as the sum of multiple steps involving ligand coupling/decoupling in the bound and unbound states. This method could be seen as an enhanced sampling approach as it allows one to preserve the binding mode of the decoupled ligand. It is broadly acknowledged as more robust compared to molecular docking and overall is in the best agreement with the experiment [21].
The energy differences calculated for each stage of the alchemical thermodynamic cycle (see Section 3.5.3 for detailed explanation) are listed in Table 2. The absolute binding energies for

FEP Calculations
The alchemical free energy perturbation (FEP) method is based on a non-physical thermodynamic cycle where the binding free energy is computed as the sum of multiple steps involving ligand coupling/decoupling in the bound and unbound states. This method could be seen as an enhanced sampling approach as it allows one to preserve the binding mode of the decoupled ligand. It is broadly acknowledged as more robust compared to molecular docking and overall is in the best agreement with the experiment [21].
The energy differences calculated for each stage of the alchemical thermodynamic cycle (see Section 3.5.3 for detailed explanation) are listed in Table 2. The absolute binding energies for

FEP Calculations
The alchemical free energy perturbation (FEP) method is based on a non-physical thermodynamic cycle where the binding free energy is computed as the sum of multiple steps involving ligand coupling/decoupling in the bound and unbound states. This method could be seen as an enhanced sampling approach as it allows one to preserve the binding mode of the decoupled ligand. It is broadly acknowledged as more robust compared to molecular docking and overall is in the best agreement with the experiment [21].
The energy differences calculated for each stage of the alchemical thermodynamic cycle (see Section 3.5.3 for detailed explanation) are listed in Table 2. The absolute binding energies for compounds A1 and A3 having very similar scaffolds were estimated as −10.8 ± 0.2 kcal/mol and −4.0 ± 0.2 kcal/mol respectively. These correspond to the pK d values of 7.9 and 2.9, in reasonable agreement with the approximate experimental pIC 50 values of 8.0 and 5.0 for compounds A1 and A3, respectively. These results suggest a better explanation of the observed difference in the inhibitory activity of the compounds A1 and A3 compared to the MM-PBSA method. According to the calculated FEP values, the free energy of decoupling and restraining the ligand in a complex is the major source of the overall differences in binding energies. In contrast to the docking results, these estimates suggest that compound A1 likely has more favorable interactions in the binding site compared to compound A3 which can be attributed to the differences in the binding modes and the presence of additional polar groups. The energy differences for the decoupling from solvent are dominated by the Coulomb term and have an opposite effect on the total binding energy (probably also attributable to additional polar interactions for A1) but their absolute size is significantly smaller. Finally, the difference in the restraint energies is negligible, apparently reflecting similar conformational flexibility of the compounds.
As an additional check of the FEP applicability to compounds ranking, we also performed the calculations for compound A2 and A7. Despite their good scores estimated by molecular docking, ML model and MM-PBSA approaches, the compound A2 shows only very weak inhibitory activity while for the compound A7 no inhibition was detected. The resulting free energy values are listed in Table 2. The total free energy of binding can be estimated as −8.2 ± 0.2 kcal/mol for compound A2 and as 7.6 ± 0.2 kcal/mol for compound A7. Surprisingly, the calculated binding energy for the virtually inactive compound A2 appears to be better than for the discovered inhibitor A3. This issue is most probably explained by the insufficient sampling during the FEP calculations, which could in principle be addressed by extension of the simulation. The observed inhibitory activity may also be affected by other processes besides the ligand-target binding. Nevertheless, even with such discrepancy, the binding energies predicted by FEP are in better agreement with the experiment compared to the ones provided by the previously applied methods, since they correctly rank compound A1 as much more potent than either A2 or A3. On the other hand, although the compound A7 shares substantial scaffold and structure similarity with compounds A1 and A3, its total free energy of binding was predicted to be positive indicating unfavorable thermodynamics of binding, primarily attributable to the ligand decoupling in solution. This is also in good agreement with the experimental inhibition data.
Overall, the predictions based on FEP calculations are not quite perfect but more reliable compared to the docking and MM-PBSA methods. That could be related to the explicit accounting for the solvation effects, in contrast to the docking-based methods or continuum solvation models used by the MM-PBSA approach. These preliminary results suggest that the FEP-calculated binding energies could potentially be used to guide lead optimization of novel TNKS inhibitors. However, a more detailed and thorough study is required in order to choose optimal parameters of the simulation. Also, it is highly desirable to conduct a separate study to evaluate the correlation of FEP energies with available experimental activity values of the known inhibitors. Alternatively, the novel machine learning approach to the analysis of the molecular dynamics trajectories [22] could be applied to predict the inhibitor activities.

Virtual Screening Library
The library of compounds for virtual screening was formed based on the ZINC 15 database [17]. Non-reactive compounds (anodyne reactivity as defined in the database) were selected from the predefined drug-like subset. These compounds were further filtered with respect to potential suppliers in order to ensure their availability for ordering in our region. The ChemAxon Instant JChem 19.23 software package (https://chemaxon.com, ChemAxon, Budapest, Hungary) was used for structure database management. The 3D coordinates of the downloaded chemical structures were generated from the SMILES format using RDKit version 2019.03.4 Python library (http://www.rdkit.org). The resulting set of the compounds included about 1.7M distinct chemical structures.

Molecular Docking
Based on our previous study [13], the structure of the TNKS2-ligand complex (PDB: 4N4V) was chosen for optimal performance of molecular docking-based virtual screening. Semi-rigid docking procedure was performed by means of Smina version 2019.10 software [23]. The search box was centered on the bound ligand and extended by default offset of 10 Å from the ligand in every direction. The Vinardo scoring function was used to find a binding pose and the exhaustiveness parameter was increased to 20 from the default value of 8. The compounds having the scoring function value below −12 kcal/mol were selected.

Prediction of Physicochemical and ADMET Properties
The assessment of a number of key physicochemical properties, pharmacokinetic parameters and toxicity endpoints was carried out for the compounds that passed the initial docking-based filter. The following predicted properties were considered as acceptable at this early stage of the search for lead compounds-molecular weight MW < 600; lipophilicity LogP < 6 (estimated using our previously developed model and the ALogPS 3.0 model available on the OCHEM (online chemical modeling environment) platform [24]); aqueous solubility > 10 −5 M (estimated using the ALogPS 3.0 model available on the OCHEM platform [24]); human intestinal absorption HIA > 75% (estimated using our previously developed model [25]); low risk of hERG-mediated cardiac toxicity (channel affinity pK i < 6 and inhibitory activity pIC 50 < 6 as estimated using our previously developed models [26]). The blood-brain barrier permeability (LogBB) values were also estimated using our previously developed model [27] to provide additional data for expert evaluation and selection of the compounds. The predictive HIA, hERG and LogBB models based on accurate and representative training sets, fragmental descriptors and artificial neural networks are available in the integrated online Molecules 2020, 25, 3171 9 of 15 service for ADMET (absorption, distribution, metabolism, excretion, toxicity) properties prediction (ADMET Prediction Service, http://qsar.chem.msu.ru/admet/). The predicted properties for the selected compounds are listed in the Appendix B (Table A2).

Biological Evaluation
The tankyrase enzyme activity was determined in vitro using the immunochemical assay based on measuring the accumulation of poly(ADP-ribose) (PAR) in the course of the PARP enzymatic reaction. PAR sorbed on nitrocellulose binds murine anti-PAR antibodies that, in turn, bind to peroxidase-conjugated antibodies against murine immunoglobulins. Peroxidase activity is determined by the chemiluminescent method via luminol oxidation, with the resulting luminescent product illuminating the X-ray film.
The tankyrase 2 catalytic fragment (human, recombinant, residues 849-1166) and its known inhibitor XAV939 were purchased from Sigma-Aldrich (St. Louis, MO, USA). Poly(ADP-ribose) polymerase 1 (PARP1) from the HT F Homogeneous PARP Inhibition Assay Kit (Amsbio) was used as a positive control. The control inhibitor or the test compounds (InterBioScreen, Moscow, Russia) were dissolved in dimethyl sulfoxide (DMSO). An enzyme was incubated in the presence or absence of the compounds in a reaction mixture comprising 100 mM Tris-HCl (pH 8), 2 mM MgCl 2 , 2% DMSO and NAD + for 30 min at 25 • C. The total volume of the mixture was 5 µL. Upon completion of the enzymatic reactions, the reaction mixture was applied in duplicates as dots on nitrocellulose matrix (Protran Nitrocellulose Transfer Membrane BA85, Whatman) that was then successively treated with the bovine serum albumin (BSA) in 0.9% NaCl and the Anti-PAR (Ab-1) Mouse mAb (10H) (Sigma-Aldrich) diluted 1:600 in the same solution. The antigen-antibody reaction was carried out for at least 1 h at room temperature or overnight at 4 • C. Next, the primary antibodies were washed off with the TBST buffer (10 mM Tris-HCl, 0.9% NaCl, 0.1% Tween 20, pH 7.4) and nitrocellulose was treated with the peroxidase-conjugated goat anti-mouse IgG antibodies (Fab specific) (Sigma-Aldrich) diluted 1:60,000. After one-hour incubation, the secondary antibodies were also washed off with TBST, nitrocellulose matrix was transferred to the solution of SuperSignal West Pico Chemiluminescent Substrate (Thermo Scientific) and exposed to an X-ray film for signal detection. The signal was quantified using ImageJ 1.52a software [28].
The IC 50 values were estimated from the concentration-response curves using the in-house Python script based on the SciPy library 1.4.1 (scipy.optimize.curve_fit function) by fitting the four parameter logistic (4PL) regression model [29]: where Y is the measured response, C is the inhibitor concentration, MIN and MAX are the minimum and maximum observed response values, IC 50 is the half-maximal inhibitory concentration and H is the Hill coefficient. The model parameters are listed in the Appendix A (Table A1).

Basic Molecular Dynamics Simulation Protocol
Initial structures of the tankyrase complexes with inhibitors were obtained by molecular docking of the compounds into the tankyrase-2 enzyme (PDB: 4TK5). The missing residues in the experimental protein structure were reconstructed by homology modeling using the MODELLER v.9.19 program [30] and the unresolved amino acid side chains were modeled using the Dunbrack rotamer library [31]. The atomic charges for the ligand molecules were calculated by the AM1-BCC method in the antechamber program from the AmberTools 16 package [32]. The ligand parameters were modeled using the General Amber Force Field (GAFF) [33]. Amber topology and coordinate files using explicit solvent (TIP3P water model, 0.15M NaCl) and dodecahedron periodic box with 10 Å margin were created by the parmchk and tleap programs from the AmberTools package using the AMBER99SB force field [34]. The data files were converted to the GROMACS format using the acpype script [35]. All subsequent molecular dynamics simulations were performed using GROMACS 2019.4 [36] software. After the initial relaxation by the steepest descent minimization algorithm, the system was subjected to temperature stabilization (300 K) in the NVT ensemble for 100 ps (50,000 steps of 2 fs each) using the modified Berendsen thermostat. Then the molecular dynamics simulation was performed in the NPT ensemble for 100 ps (50,000 steps of 2 fs each) using the Berendsen barostat in order to stabilize the pressure at 1 atm. After the temperature and pressure were equalized, a production molecular dynamics simulation was performed using v-rescale thermostat and Parrinello-Rahman barostat on the NVIDIA GTX 1080Ti high-performance GPU.

MM-PBSA Calculations
In order to estimate the binding energies of the tankyrase-inhibitor complexes, a preliminary molecular dynamics simulation of 30 ns was performed according to the basic protocol described above. The resulting system state was used as a starting point for ten independent runs of 5 ns each as suggested in the work [19]. The bootstrap procedure [37] was used to estimate the mean and confidence interval RMSD values of the calculated binding energies for each run that were then aggregated using mean and L2-norm, respectively.
The binding energy was calculated by the single-trajectory MM-PBSA method using the g_mmpbsa standalone software [37]. To reduce the computational cost, 251 equidistant frames from a time interval of 1-5 ns in the 5 ns simulation trajectory were used. The solute dielectric constant was set to 4. The binding free energy was estimated as: where ∆G bind,aq is the binding free energy, ∆E MM is the molecular mechanics energy change consisting of the electrostatic energy change and the van der Waals energy change, ∆G solv is the solvation free energy change consisting of polar and non-polar terms (in this study, the Poisson-Boltzmann solver was used for the polar term and the SASA model for non-polar contribution) and T∆S is the entropy term. Due to high computational demand, the entropy term was omitted during this study (moreover, a number of works suggest that its inclusion does not lead to any improvement in the ranking of relative binding affinities [38,39]).

Absolute FEP Calculations
The alchemical free energy perturbation (FEP) method is based on a non-physical thermodynamic cycle comprising the following states-(1) physical unbound state, (2) alchemical state where the ligand is decoupled from the solution, (3) alchemical state where the ligand is decoupled and restrained in the binding site, (4) physical bound state. This method could be seen as an enhanced sampling approach as it allows to preserve the binding mode of the decoupled ligand. The binding free energy is calculated as: where ∆G 0 bind is the free energy of binding, ∆G prot elec+vdw+restr is the free energy of decoupling and restraining the ligand in a complex, ∆G solv elec+vdw is the free energy of decoupling the ligand in solution and ∆G solv restr is the free energy for restraining the decoupled ligand in solution (due to the closed-form analytical expression, it demands very little computational overhead).
The initial states of the tankyrase-inhibitor complexes for the FEP calculations were obtained from the same preliminary 30 ns molecular dynamics simulations that were used for the MM-PBSA calculations. Absolute binding free energy calculations were performed using standard GROMACS tools according to the methodology [40]. The free energy differences across the intermediate states were calculated using the Multistate Bennett Acceptance Ratio (MBAR) approach [41] implemented in the alchemical_analysis.py Python script [42].
Ligand decoupling from a complex involves gradually turning off the Coulomb and Lennard-Jones interactions as well as restraints over a total of 30 windows controlled by linearly spaced lambda values. For each lambda value, the energy minimization, NVT and NPT equilibration and 1 ns production run were performed, followed by the MBAR-based estimation of the free energy. Similarly, ligand decoupling from solvent comprises a total of 20 windows with linearly spaced lambda values (in this case, only Coulomb and Lennard-Jones interactions were turned off). The length of the production runs for this step was set to 5 ns. For consistency reasons, the MBAR free energy estimation was also performed. The set of the binding site restraints was defined by one distance, two angle and three dihedral harmonic potentials according to work by Boresch et al. [43].

Conclusions
In the course of this study, the virtual screening of more than 1.7 million commercially available compounds was performed based on the molecular docking results, predicted physico-chemical and ADMET properties and molecular dynamics simulations. In the in vitro biological evaluation, two out of seven selected compounds showed a decent level of inhibitory activity against the TNKS2 enzyme with the IC 50 values of less than 10 nM and 10 µM.
The in-depth retrospective analysis revealed that the molecular docking-based scoring functions or MM-PBSA calculations are not suitable for accurate ranking the compounds according to their affinity values and can only be used for binary classification to select potential hits. In contrast, an alternative approach based on the free energy perturbation (FEP) calculations, although not perfect, allowed us to establish a more reliable correlation between the binding energy and inhibitory activity of four compounds and to explain some differences in their activity in spite of very similar scaffolds. That could be related to the explicit accounting for the solvation effects, in contrast to the docking-based methods or continuum solvation models used by the MM-PBSA approach.
Thus, instead of the fast and cheap docking-based scoring, we propose the free energy calculations as a possible method for the subsequent lead optimization stages. However, it should be applied only after additional detailed study aiming to select optimal parameters and evaluate its ranking performance on a large set of known inhibitors.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

Appendix A. Concentration-Response Curves for Compounds A1 and A3
The concentration-response curves for compounds A1 and A3 are shown in Figure A1. Measurements were carried out at 100 µM NAD + concentration according to the procedure described in Materials and Methods (Section 3.4). The response values are the means of duplicate points. The IC50 values were estimated using in-house Python script based on the SciPy library 1.4.1 (scipy.optimize.curve_fit function) by fitting the four-parameter logistic (4PL) regression model: where is the measured response, is the inhibitor concentration, and are the minimum and maximum observed response values, is the half-maximal inhibitory concentration and is the Hill coefficient. The model parameters are listed in Table A1. Taking into account small number of repeated experiments and high data variance, the IC50 values can be cautiously estimated as less than 10 nM for compound A1 and less than 10 μM for compound A3.

Appendix B. Prediction of Physicochemical and ADMET Properties of Compounds A1-A7
The prediction of key physicochemical and ADMET properties during virtual screening was performed using the methods and procedure described in Materials and Methods (Section 3.3). The predicted values for compounds A1-A7 selected by virtual screening are listed in Table B1.  The IC 50 values were estimated using in-house Python script based on the SciPy library 1.4.1 (scipy.optimize.curve_fit function) by fitting the four-parameter logistic (4PL) regression model: where Y is the measured response, C is the inhibitor concentration, MIN and MAX are the minimum and maximum observed response values, IC 50 is the half-maximal inhibitory concentration and H is the Hill coefficient. The model parameters are listed in Table A1. Taking into account small number of repeated experiments and high data variance, the IC 50 values can be cautiously estimated as less than 10 nM for compound A1 and less than 10 µM for compound A3.

Appendix B. Prediction of Physicochemical and ADMET Properties of Compounds A1-A7
The prediction of key physicochemical and ADMET properties during virtual screening was performed using the methods and procedure described in Materials and Methods (Section 3.3). The predicted values for compounds A1-A7 selected by virtual screening are listed in Table A2.