In silico screening of small molecule modulators and their binding studies against human sirtuin-6 protein

Abstract Sirtuin-6 (SIRT6), class III family of deacetylase regulates several biological functions, including transcriptional repression, telomere maintenance, and DNA repair. It is unique among sirtuin family members with diverse enzymatic functions: mono-ADP-ribosylase, deacetylase and defatty-acylase. The studies so far implicated SIRT6 role in lifespan extension, tumor suppression, and is considered as an attractive drug target for aging-related disease. In this study, we have carried out in silico screening for human SIRT6 modulators using NCI Diversity Set III library, molecular dynamic (MD) simulations to analyze the protein-ligand interaction, and validated their binding-affinity ( ) using MicroScale Thermophoresis. This study yielded two novel compounds, ((3Z)-3-((4-(dimethylamino)phenyl)methylidene)-5-(5,6,7,8-tetrahydronaphthalen-2-yl)furan-2-one and 5-phenyl-2-(5-phenyl-2,3-dihydro-1,3-benzoxazol-2-yl)-2,3-dihydro-1,3-benzoxazole showing high-affinity interaction for SIRT6. The structural analysis from MD simulation suggests both compounds might act as substrate-analogs or mimic the nicotinamide binding. On considering the uniqueness of SIRT6 substrate binding acyl channel among sirtuin family member, binding of both compounds to the above site suggesting their specificity for SIRT6 isoform. Therefore, it may form the basis for the development of potential modulators for human SIRT6. Communicated by Ramaswamy H. Sarma


Introduction
Sirtuin's are evolutionarily conserved, NAD þ -dependent, class III family of deacetylase enzymes. They are extensively studied for their role in various human diseases (Kugel & Mostoslavsky, 2014;Sanders et al., 2010;Sauve et al., 2006). Mammals have seven sirtuin family members (SIRT1-7) having distinct functions and subcellular localization (Michan & Sinclair, 2007;Sacconnay et al., 2016). They have conserved catalytic core region consisting of $275 amino acids at the center and is flanked with N-and C-terminal regions varying in length and sequence (Sanders et al., 2010). The crystal structures of sirtuin catalytic core regions are known, and it consists of two domains; a large and structurally homologous Rossmann-fold domain, and a more structurally diverse, smaller, zinc-binding domain (Finnin et al., 2001;Min et al., 2001;Pan et al., 2011). The loops connecting these domains form a cleft, where the co-factor NAD þ and acetyl-lysine containing peptide substrate enter from opposite sides and binds to the enzyme. Moreover, the amino acids within the cleft between the two domains share the highest sequence conservation among the sirtuin family (Sacconnay et al., 2016;Sanders et al., 2010).
The SIRT6 has diverse enzymatic activities among sirtuin family members: deacetylation, defatty-acylation and mono-ADP-ribosylation (Feldman et al., 2013;Jiang et al., 2013;Mao et al., 2011;Zhang et al., 2016). The full-length SITR6 is about 355 amino acids in length, consists of the conserved catalytic domain at N-terminus (1-296 amino acids), and the flanking tail at C-terminus (297-355 amino acids) is prolinerich and hydrophobic in nature with unknown function (Michan & Sinclair, 2007). The crystal structure of the SIRT6 with H3K9-myristoylated peptide revealed that occupation of the substrate in the hydrophobic pocket could lead to ordered N-terminal region (1-12 residues) otherwise above residues density are not visible in the SIRT6-ADPr complex structure (Jiang et al., 2013;Pan et al., 2011). SIRT6 is localized in the nucleus and is known to bind tightly to chromatin (Kugel & Mostoslavsky, 2014;Mostoslavsky et al., 2006). Unlike SIRT1, the SIRT6 deacetylation activity is substrate-specific on histone H3K9 and H3K56 (Michishita et al., 2008(Michishita et al., , 2009Yang et al., 2009). Moreover, SIRT6 hydrolysis longchain fatty acyl-lysine more efficient compared to acetylatedpeptide (Feldman et al., 2013;Zhang et al., 2016). SIRT6 is highly expressed in brain tissues in cortical and hippocampal regions and enriched in the synaptosomal membrane fraction (Cardinale et al., 2015). SIRT6 knockout mice exhibits severe metabolic abnormality, accelerated aging phenotype and die prematurely. Brain-specific SIRT6deficient mice survive, but lack of SIRT6 provokes neurodegeneration by increasing DNA damage, apoptosis and toxic Tau phosphorylation (Kaluski et al., 2017;Mao et al., 2011). SIRT6 is also shown to associate specifically with telomeres, and its depletion leads to telomere dysfunction with end-toend chromosomal fusions and premature cellular senescence (Michishita et al., 2008). Moreover, SIRT6-depleted cells exhibit abnormal telomere structure that resembles Werner syndrome defects, a premature ageing disorder.
So far, sirtuin modulators are identified mainly by screening chemical library, catalytic mechanism-based design approaches, often combined with structure-activity relationship investigations. During last five years, the crystal structures of SIRT6 (12-296) catalytic domain consisting nonhydrolyzable co-factor analogs ADP-ribose (ADPr) in complex with modulators structure are reported (Klein & Denu, 2020;You et al., 2017You et al., , 2019. Nevertheless, only a few reported SIRT6 modulators have reached to the stage of animal studies; hence, identifying novel human SIRT6 modulators are essential (Mautone et al., 2020). In this study, we carried out in silico screening of human SIRT6 modulator using NCI library, followed by molecular dynamics simulation studies, and validated their binding affinity (K d ) using MicroScale Thermophoresis (MST). This study yields two novel compounds with good binding affinity for SIRT6 both in the presence or absence of ADPr.

Ligand preparation
The Diversity Set III of DTP program from NCI compound library was used to screen the potential competitive inhibitors of SIRT6. The compounds' 2 D structures were drawn to get the smiles format using ChemDraw Ultra version 8.0 (Mendelsohn, 2004). The smiles of respective compounds were converted to their pdb format, and their 3 D coordinates were generated using Open Babel version 2.4.1 (O'Boyle et al., 2011). The ligands were then preprocessed for docking using prepare_ligand4.py python script in AutoDock MGL Tools version 1.5.6 (Morris et al., 2009). This script formats the ligand file in pdbqt format necessary for AutoDock, the format includes partial charges, correct AutoDock atom types and torsional flexibility of atoms.

Receptor preparation
The target protein structure PDB ID 5MFP (1.98 Å) present in the complex with ADPr, a moiety of the co-substrate NAD þ was retrieved from Protein Data Bank (www.rcsb.org). The selected SIRT6 structure was co-crystallized with an activator; therefore, prior to preprocessing the ligand, solvent molecules and ions were removed. The crystallographic structure was minimized, gasteiger charges were computed, and nonpolar hydrogens were merged using the AutoDock Tools.

Docking
Analyzing the protein ligand interactions is imperative for drug designing and discovery and to aid this, computational methods like docking and simulation are used. Docking helps in the prediction of different binding modes and binding affinity of the ligand toward protein. For this study, the compounds were subjected to blind docking using AutoDock Vina version 1.1.2 (Trott & Olson, 2010). The geometric center of grid box was placed at center of the macromolecule with co-ordinates x ¼ 113.416, y ¼ 9.064, and z ¼ À17.919. The grid box of 60 Å Â 70 Å Â 50 Å points in x, y, z-direction respectively with 1 Å spacing, was generated to cover the entire protein. The Lamarckian genetic algorithm was used to predict 9 conformations per ligand. The exhaustiveness of the global search was kept at 9 and rest of the parameters were kept at default.
The output ligand pdbqt files were processed using process_VinaResults.py script, thereby yielding each model into separate files with details of their close contacts, ligand efficiency and so on. Poses are ranked based on their AutoDock scoring function with the best binding mode having the least (most negative) binding affinity. Finally, the compounds were evaluated and filtered based on their binding site, affinity and interactions with the receptor molecule.

Molecular dynamics (MD) simulations
MD simulations were performed to investigate the structural dynamics of protein-ligand complexes on Desmond 4.4 software (Bowers et al., 2006) incorporating OPLS_2005 (Shivakumar et al., 2010) force field. The SIRT6 þ ADPr (apo protein) and SIRT6 þ ADPr þ optimal conformation of the ligand compound from docking studies (protein-ligand complex) were subjected to MD simulations. Before building the model system, the protein was prepared by addition of hydrogen atoms, assigning bond orders, and optimizing the hydrogen bonding network. The complex systems were solvated in an orthorhombic box with a buffer size of 10 Å between the solute structure and the box edge. The solutes were reoriented by minimizing the volume of the box. The predefined TIP3P water model (Jorgensen et al., 1983) was used to simulate the explicit water molecules. Along with adding an Na þ ion to make electrically neural system before simulations, 150 mM salt concentration was set up under periodic boundary conditions. The solvated and neutralized system was relaxed into a local energy minima by minimization using steepest decent and limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithms for 100ps. The temperature was controlled at 300 K by Nose-Hoover Thermostat. Similarly, the pressure was maintained at 1.01325 bar by Martina-Tobias-Klein method. The SHAKE algorithm was applied to constraint the bonds involving hydrogen atoms. Short-range coulombic interactions were cut off after 9 Å, and long-range interactions were handled by Particle Mesh Ewald summation method. This was followed by a production run of 100 ns with trajectory recording every 50ps. The generated 2000 frame trajectories were analyzed using Simulation Interaction Diagram (SID), Simulation Event Analysis (SEA) modules of Desmond and Maestro version 12.3.013 (Release 2020-1) modules from the Schr€ odinger Suite.

Binding free energy estimation
For the binding free energy calculations, the MMPBSA.py script from the freely available AmberTools20 (Case et al., 2020) suite was applied on the equally spaced snapshots generated from 200 ns MD trajectories of each complex. The snapshots were striped to remove solvent molecules aiming to replace the explicit solvent with the implicit Generalized Born model. The ff14SB (Maier et al., 2015) and GAFF force fields were used for protein and compound parameterization respectively. The force field parameter modification (.frcmod) files with custom parameters for nonstandard residues (i.e. the heteroatoms -ADPr, Zn 2þ and compounds) were created using antechamber which were then fed to LEaP program for further processes. The necessary input AMBER topology (prmtop) and coordinate files (inpcrd) were generated for the receptor, ligands and respective complexes using tleap with mbondi2 radii. The calculations were performed with Generalized Born implicit solvent model igb ¼ 2 and salt concentration of 150 mM using the Molecular Mechanism-Generalized Born Surface Area (MM-GBSA) (Genheden & Ryde, 2015; approach. Rest of the parameters were kept default in the MM-PBSA input file. The free energy difference for each complex is calculated using the following Equation (1) : DG gas is binding free energy between receptor and ligand in gas phase; DG solv,complex , DG solv,protein DG solv,ligand are solvation free energies of complex, protein and ligand respectively. The DG gas considers the average molecular mechanical energy between protein and ligand while also considering the absolute entropy change, Equation (2): Since the aim of the study was to obtain an estimated binding free energy, the TDS term is ignored in Equation (2). DE electrostatic is the electrostatic energy, while DE vdw is the van der Waals interaction energy. The DG solv solvation energy term for each system is calculated here using Generalized Born equation for electrostatic contributions and an empirical term for non-polar contributions, Equation (3):

Cloning, expression and purification of human SIRT6
The human full-length SIRT6 construct was obtained from Addgene, and was used as a template to PCR amplify the truncated SIRT6 (13-296). It was sub-cloned into pET28a vector at NdeI and BamHI restriction sites and sequences verified before the bacterial expression. The truncated SIRT6 constructs with (His) 6 -tag at N-terminus was expressed in E. coli expression strain BL21(DE3) cells. The cells were grown at 37 C until density reaches OD 600 0.6. The cells were then induced with 0.5 M IPTG and grown at 16 C post induction for 18 hrs. The cells were then harvested and resuspended in lysis buffer A (50 mM Tris-HCl, pH 8.0, 500 mM NaCl, 2 mM b-mercaptoethanol, PMSF and Benzamidine), and lysed by sonication. The SIRT6 in the supernatant was affinity purified using HiTrap IMAC FF column (GE Healthcare, USA), which is followed by overnight thrombin (Sigma) digestion to cleave the (His) 6 -tag at 4 C. The SIRT6 was further purified via HiTrap Heparin HP (GE Healthcare, USA) column using a linear gradient from 0 to 1000 mM NaCl in 50 mM Tris-HCl, pH 8.0 and 2 mM b-mercaptoethanol. The fractions containing purified SIRT6 were pooled, concentrated, and loaded into size-exclusion chromatography (HiLoad Superdex S-75 16/60; GE Healthcare, USA). The concentration of final purified SIRT6 protein was determined by the NanoDrop TM OneC Spectrophotometer (Thermo Scientific TM ) using extinction co-efficient at OD 280 , concentrated the protein, and stored in À80 C until further use.

Microscale thermophoresis (MST)
The MST experiments were performed according to the NanoTemper technologies protocol and affinities calculated using Monolith NT.115 (Red/blue) instrument (NanoTemper Technologies GmbH, Munich, Germany). For MST experiments, the SIRT6 (13-296) construct was labeled using lysin reactive Monolith Protein Labeling Kit RED-NHS 2nd Generation (Amine Reactive) (NanoTemper Technologies) according to the supplier labeling protocol. The dye carries a reactive NHS-ester group that reacts with primary amines (lysine residues) to form a covalent bond. The final concentrations of NT-647 NHS labeled SIRT6 was 125 nM, chosen on the bases of fluorescence intensities from the pretest assay setup. The pretest essay helps to choose the optimal concentration of the labeled molecule and the buffer in which the sample is homogeneous. The unlabeled individual NCI Compounds #18, #12, #4, #32 and #35 were chosen for binding studies with the above-mentioned SIRT6 construct based on their energy scores from docking analysis. Quercetin and Catechin gallate, known SIRT6 binding compounds, were used as positive controls for the study (You et al., 2019). For this study, the SIRT6 was freshly prepared in the assay buffer (20 mM Tris-HCl pH 8.0, 150 mM NaCl, 1.0% DMSO) supplemented with 0.05% Tween-20 to improve the homogeneity of the sample. The experiments were performed by incubating varying concentrations of target compounds with a constant concentration of respective labeled SIRT6 construct at room temperature for 10 min before recording the measurement using NT.115 standard treated capillaries (NanoTemper Technologies). The data were acquired at 25 C using LED power in the range of 90% with MO.Control 1.5.3 (NanoTemper Technologies GmbH, Munich, Germany). MST setup allows detection of protein aggregation and adhesion effects of the assay samples on bases of the MST traces and capillary scan. The assay setup for each binding affinity check was analyzed, and no sign of aggregation was observed. Further analyzed the recorded data with MO.Affinity Analysis 2.2.7 NanoTemper Technologies GmbH) and determined the K d values. The measured values reported as MST dependent responses of the fluorophore to the temperature jump. The manuscript figures were prepared using Graphpad Prism 8.4.3 (GraphPad, San Diego, California).

Results and discussion
The SIRT6 is an attractive therapeutic target for various human diseases, including cancer and neurodegenerative diseases (Kaluski et al., 2017;Kugel & Mostoslavsky, 2014;Kugel et al., 2015). It has a unique, elongated hydrophobic substrate-binding pocket that accommodates long-chain acyl substrates (Jiang et al., 2013;Klein & Denu, 2020). In this study, we carried out in silico screening, MD simulation, in vitro binding studies, which yielded two novel compounds with good binding affinity for SIRT6. Intriguingly, both compound's binding sites are different from the SIRT6 inhibitor and activator molecules binding site known so far. We have also modeled in the six missing residues in the loop region (R2 region) of the SIRT6 using ModLoop (Fiser et al., 2000;Fiser & Sali, 2003) which uses the loop modeling routine in MODELLER (Webb & Sali, 2016) and predicts the different conformations of the loop. For the prediction of coordinates of the six missing residues PDB ID 5 Â 16 was used as the template. The best fit output from ModLoop was used to redo the simulation studies of 100 ns (see below). Finally, MD simulation was also carried out for loop-modelled structure for cross-validation.

Analysis of molecular docking results
In this study, we have carried out in silico screening of SIRT6 (PDB ID: 5MFP; Figure 1) with the NCI Diversity Set III comprising 1600 compounds using AutoDock Vina (Trott & Olson, 2010). The blind docking covers the entire receptor for ligand docking without specifying a single binding site. This helps identify the potential binding sites of ligands on the whole protein and select accordingly without bias toward already known sites. The results were examined by visualizing the interaction fields with the amino acids surrounding the active site and substrate-binding groove. Based on the energy scores and their interaction with SIRT6, the identified top five compounds preferentially binding to the substrate-binding site were selected for further biophysical studies. The details of the selected compounds are shown in Table 1. Among those chosen compounds, four compounds were bound in the substrate-binding groove, and one compound binds adjacent to the adenine binding region of NAD þ co-factor (Site 'A'; Figure 1). Out of nine poses generated for each ligand most bound to the substrate binding and the co-factor binding sites; except very few poses binding to different regions with low energy. We observed that most of the interactions were principally van der Waals and hydrophobic interactions in nature ( Figure S1). For ligand efficiency the python script process_VinaResults.py from MGLTools-1.5.6 of AutoDock was utilized. Since the ligand efficiency calculations depend on the heavy atom count in a ligand; the size independent ligand efficiency (SILF) calculations (Nissink, 2009) were also performed for more accurate ligand efficiency determination ( Table 1).
The docked structural pose of the compound #18 (NSC127133) forms two hydrogen bonds with the backbone atoms of Arg220 via the hydroxyl group of the carboxy phenyl ring, and a weak hydrogen bond occurs between the oxygen atom of phenanthridinone and Trp188 indole side chain. Trp188 engages in p-p T-shaped interactions with the phenanthridinone group from one side while from opposite site His133, Leu186 and Val115 involve in different aromatic interactions (Figure 2A). The long compound 12 (NSC84100) consisting of 2 dihydro-5-phenylbenzoxazole moieties is majorly held in the hydrophobic groove by multiple aromatic interactions with His133, Leu186 and Trp188 ( Figure 2B). However, half of the compound protrudes toward the peptide substrate entry site thus exposing to the solvent. Here, it orients to form hydrophobic interactions with residues Leu192 and Pro221 lining the substrate entry channel. The compound #4 (NSC33575) was relatively deep-seated and spanning the entire channel ( Figure 2C). The 2-furanone ring surrounded laterally by His133, and Trp188 sits above ribose moiety of ADPr, with the tetrahydronaphthalene moiety occupying the probable NAD þ nicotinamide binding site i.e. the 'B' site of NAD þ . Toward the substrate exit site, edge-toface and p-alkyl aromatic interactions of tetrahydronaphthalene with Phe64, Trp71 and Val115 were observed. The dimethylamino group attached to the phenyl ring shows potential polar interaction with the backbone oxygen of Arg220 ( Figure 2C). Analysis of the predicted orientation of compound #32 (NSC326182) reveals interactions with the nonpolar residues, Val115, Met157, Ile185, Leu186, Trp188 and Ile219 lining the aperture ( Figure 2D). It is anchored at the position by a network of hydrophobic interactions. For the phenanthroimidazole derivative, compound #35 (NSC332670), most of the poses occurred near the substrate binding channel though the best pose was observed at the region adjoining NAD þ adenine binding site ( Figure 2E). We observed that the other predicted poses for compound #35 shows binding at the substrate channel like other mentioned compounds. Therefore, we performed simulations using this conformation. The hydroxy group of phenyl ring forms a strong triad of hydrogen bonds with backbone atoms of Asp187, Trp188 and Asp190 ( Figure 2F). The phenanthroimidazole ring secured in the peptide entry site by aromatic interactions with Ile219 and Arg220, the residues flanking the acyl channel. We speculate that the colossal ring system might be the reason for the hindrance of the compound's entry inside the channel. Comparisons between these modulators reveal firm anchorage in the hydrophobic pocket due to multiple aromatic interactions with the respective ring systems.

Analysis of molecular dynamics simulation results
The docking studies were combined with more extensive MD simulation studies to investigate the SIRT6-ligand complex's structural dynamics. Along with the examination of fundamental parameters, detailed interaction analysis was conducted to identify the best modulator.
For the 100 ns simulation run generating 2000 frames of the protein þ ADPr (apo protein) and the protein þ ADPr þ compound (protein-ligand complex) structures, backbone Root Mean Square Deviation (RMSD) was calculated ( Figure 3A). Upon comparison with the apo protein, the complexes' RMSD values reveal the protein structure to be unaffected by binding of the modulators. SIRT6 complexes with compound #4 and compound #12 attain conformational stability at $3.6 Å than the apo form which attains stability at $2.5 Å ( Table 2). All the complexes show a plateau in RMSD values and convergence of structure within the initial 5 ns of the 100 ns simulation run.
To investigate the change in ligand binding site over the entire simulation run, we calculated ligand RMSD with respect to the protein structure ( Figure 3B). The compound #18, compound #12, compound #4, and compound #32 shows high stability in their initial binding site throughout the simulation. Whereas compound #35 deviates highly after 26 ns of simulation and tends to deviate more after 80 ns where the RMSD value goes up to 9 Å ($5 Å deviation from initial frame) indicating that it may not be stable compared to other discussed compounds.
Root Mean Square Fluctuations (RMSF) indicate the local structural flexibility of a protein. The RMSF profiles of the backbone atoms of protein-ligand complexes ( Figure 3C) indicate high flexibility in regions corresponding to residues 65-90 and 160-182 in apo-and the complex structures. The region between 65 and 90 residues contains the a3 helix, which is involve in NAD þ binding ( Figure 1). As observed from the RMSF profiles, the fluctuation in this region is higher than the apo protein for compound #4 and compound #12, indicating that these compounds might be influencing the co-factor binding region. The second region comprising 160-182 residues form a part of the flexible loop region of the Zn 2þ binding motif leading to acyl channel residues of the large Rossmann fold domain. Strikingly the fluctuations are not caused due to binding of the modulators. Interestingly, it was noticed that residues 277-287, though corresponding to the loop region of C terminal, demonstrated higher stability in protein-ligand complexes compared to the apo protein.  (Min et al., 2001;Sacconnay et al., 2016).
The protein-ligand free energies of binding calculated on the MD trajectories are presented in Table 3 for comparative analysis. The observation indicates compound #4 has the strongest binding (-39.73 kcal/mol) toward SIRT6 followed by compound #32 (-39.34 kcal/mol). The free energy values differ from the in vitro studies which show no binding for compound #32. It was interesting to observe that compound #32 has lowest DG solvation value (6.5 kcal/mol) which leads to its overall low binding energy. Compound #12 which shows best experimental K d value (approx. 18 times lower than compound #4) falls behind compound #4 in free energy calculations but this should be noted that the energy difference between these compounds is very less. According to the calculated values compound #35 shows least favorable binding free energy (-31.92 kcal/mol) for SIRT6 compared to other compounds; this agrees with our in vitro and in silico observations. The major contributor in MMGBSA binding free energy values is van der Waal energy term. With exception of compound #32 rest of the compounds give calculated free energy values in accordance with the trend seen in experimentally determined K d value. The deviations are expected while comparing predicted vs experimental values.
Interaction analysis of the trajectories revealed that His133 involves in hydrogen bonding with compound #18's carboxyl group of the benzene ring for initial 50 ns and afterwards for last 40 ns of simulation demonstrates p-p stacking interactions with carboxy phenyl and phenanthridinone groups. The phenanthridinone ring is sandwiched by aromatic interactions on both sides with His133 and Trp188. The aromatic interactions with Trp188 can be seen throughout the 2-(2-((6-oxo-5H-phenanthridin-3yl)carbamoyl)phenyl)benzoic acid -9.3 -0.3100 À3.35 12 5-phenyl-2-(5-phenyl-2,3-dihydro-1,3-benzoxazol-2-yl)-2,3-dihydro-1,3-benzoxazole  simulation time. Additionally, polar interactions were observed with backbone of Asp187 and linker amine group of the compound. Brief water mediated interactions can be observed with Asp190 and Ser191 which are maintained for less than 20% duration. Surrounded by other non-polar amino acids, the compound shows non continuous aromatic interactions with Ile185, Leu192 and Ile219 as well. Initially, the distal end of compound #12 exposed toward the substrate entry side shows high fluctuations in an attempt to make interactions with the N terminal residues Pro10 and Phe11. The brief interaction was of p-cation and p-p stacking nature which was distorted as the compound slips deeper in the cleft and further stabilizes. We observed that the hydrophobic interactions with His133 and Trp188 of the acyl channel and Trp71 of a3 helix were responsible for pulling the compound ( Figure 4A,B). Compound #12 also shows hydrogen bonding with Asp187 via oxazole nitrogen atom ( Figure S3). Hydrogen bonding with Asp190 was also observed in the beginning of simulation for a short duration of $15ns. Analysis of the simulation trajectory exhibited continuous interactions of compound #4 with His133 and Phe64 residues causing its high stability at the binding site. However, the compound also displayed aromatic interactions with Trp188 ( Figure 4D). Hydrophobic interactions were also observed with Ile61, Trp71, Phe82, Met157, Ile185, Ile186, Leu192, Pro193, Ile219 and Pro221. These residues surround the compound and form inconsistent aromatic interactions which are scattered throughout simulation. The tetrahydronaphthalene system present at the acyl exit site is surrounded by non-polar residues Pro80, Phe82, Met136 and Met157 ( Figure 4C). The His133 engages in two types of interactions: a hydrogen bond with the furanone ring and simultaneous aromatic interactions with the dihydroxybenzenamine and furanone rings. The strong interactions with Phe64 and Trp71 only formed with compound #4 and #12 binding explain the fluctuations observed in the a3 helix of respective complexes compared with the apo structure. In contrast to rest of the compounds, simulation trajectory of compound #32 showed an additional polar interaction between the oxygen atom of benzodioxole and Ser222 residue. Inspection of the run shows hydrogen bonding with Leu186 as well but after 50 ns. Leu186 can be seen interacting with the compound throughout the simulation. It is involved not only with hydrogen bonding but also with aromatic and water mediated interactions as well. While the polar contact sums to only 20% time of trajectory, the aromatic interaction is retained for $70% duration. Compound #32 also shows aromatic interactions with His133.
As explained in the docking analysis part the fourth binding mode of compound #35 was used for dynamics studies. It forms hydrogen bonding with the backbone oxygen atom of Asp190 and Trp188 via the hydroxy group of phenol ring. The phenol ring is involved in the inconsistent aromatic hydrogen bonding with Asp190 and Arg220. The aromatic hydrogen bonding interaction with Arg220 is very short lived (C) Protein RMSF curves in apo and complex forms. R1 corresponds to region from Arg65 to Arg90; R2 corresponds to region from Lys160 to Arg182; and R3 corresponds to region from Asp277 to Pro287. (D) Positions of R1, R2 and R3 highlighted in green, blue and purple, respectively in structure. seen at the beginning of simulation only and is replaced by water mediated interaction for the last 40 ns ($20% of run). The same type of interaction is also observed with Glu189 and Leu192 in brief intervals after 40 ns of run. Compound #35 does show aromatic interactions with His133, Leu186 and polar contact with Asp187 which continues initial 30 ns run afterwards are lost. Over the simulation, the compound moves out of the channel toward the peptide entry site and is held at the entry site by its phenanthroimidazole ring system which forms hydrophobic and hydrogen bond interaction with Leu192 ( Figure S4). Thus, indicating that the compound #35 exhibits non-stable binding with the protein.
MD simulation was also carried out for loop-modelled structure for cross-validation. As speculated, the compounds bind to the same pocket and show similar interactions. The RMSD values showing the convergence of complexes also remain in the same range (not shown).
Based on structural analysis from molecular dynamics simulations, we speculate the high affinity of compound #4 and #12 might be because of their interactions with the residues of a3 helix and residues of the hydrophobic channel especially His133 and Trp188. This single helix is replaced by a loop region in other isoforms of sirtuin making it highly flexible and prone to a high degree of fluctuations (not shown). The stable helix replacement might lead to the specificity of these compounds to SIRT6. Hence, we speculate that the predicted compound #4 and #12 may be substantial modulators of SIRT6 compared to the other predicted compounds.

Binding studies of the compounds with SIRT6
For this study, we have expressed and purified SIRT6 truncated (13-296 residues) proteins using a bacterial expression system and stored in À80 C ( Figure S5). For MST experiments, the SIRT6 protein was freshly thawed and labeled using lysin reactive Monolith Protein Labeling Kit RED-NHS 2nd Generation (Amine Reactive) (NanoTemper Technologies) according to the supplier's labeling protocol. The MST studies were performed with varying concentration of the compounds (highest concentration 2 mM) incubated with SIRT6 (125 nM) in assay buffer (20 mM Tris pH 8.0, 150 mM NaCl, 1.0%DMSO) supplemented with 0.05% Tween-20. The capillary scan showed no sign of aggregation and adhesion of sample with laser on/off times for 30 s using 90% MST power. All five compounds were measured for their binding affinity using MST both in the presence and absence of ADPr. This study showed that the ADPr moiety, the hydrolyzed NAD þ product, induced co-factor binding site closure to support better binding affinity for compounds #12 and #4. Among the five compounds studied, compound #18 has a K d value of 800 mM, compounds #4 and #12 with a K d value 27 mM and 1.6 mM. Intriguingly, compound #12 and #4 showed high-affinity binding compared to compound #18, which had higher energy scores during docking studies. The compound #12 has a longer aromatic ring structure than compound #4, which likely contributed significantly to the higher affinity against SIRT6. The compound #32 and #35 had no binding with SIRT6 ( Figure 5). Therefore, the predicted compounds from the simulation data show stable binding and correlates with our biophysical assay. Furthermore, the comparative analysis with known SIRT6 modulators; Quercetin and Catechin gallate as a K d value of 49.7 mM and 100 mM, respectively ( Figure 5) (You et al., 2019). Intriguingly, the identified compounds #4 and #12 show a higher binding affinity toward the SIRT6.

Conclusion
In conclusion, we have used in silico and in vitro methods to identify low micromolar modulators of a therapeutic SIRT6 protein. We performed in silico screening of compounds from the NCI database, structural dynamics, and followed by binding affinity studies using purified bacterially expressed human SIRT6. Our simulation results revealed strong binding for compound #12 and #4 to the substrate channel and the active site, whereas non-stable binding for compound #35, which are in accordance with the MST data. Comparative interaction analysis reveals interactions with His133, Trp188 and the a3 helix might contribute to observed strong binding affinity. The compounds #4 and #12 described here can act as substrate-analogs or mimic the nicotinamide binding of this protein since engaging in interactions with residues involved in co-factor binding. Since both compounds bind to SIRT6's acyl channel, we extrapolate that both compounds might have specificity to the SIRT6 isoform. Therefore, it may form the basis for the development of potential inhibitors or activators of human SIRT6.