Synthesis of Novel Halogenated Heterocycles Based on o-Phenylenediamine and Their Interactions with the Catalytic Subunit of Protein Kinase CK2

Protein kinase CK2 is a highly pleiotropic protein kinase capable of phosphorylating hundreds of protein substrates. It is involved in numerous cellular functions, including cell viability, apoptosis, cell proliferation and survival, angiogenesis, or ER-stress response. As CK2 activity is found perturbed in many pathological states, including cancers, it becomes an attractive target for the pharma. A large number of low-mass ATP-competitive inhibitors have already been developed, the majority of them halogenated. We tested the binding of six series of halogenated heterocyclic ligands derived from the commercially available 4,5-dihalo-benzene-1,2-diamines. These ligand series were selected to enable the separation of the scaffold effect from the hydrophobic interactions attributed directly to the presence of halogen atoms. In silico molecular docking was initially applied to test the capability of each ligand for binding at the ATP-binding site of CK2. HPLC-derived ligand hydrophobicity data are compared with the binding affinity assessed by low-volume differential scanning fluorimetry (nanoDSF). We identified three promising ligand scaffolds, two of which have not yet been described as CK2 inhibitors but may lead to potent CK2 kinase inhibitors. The inhibitory activity against CK2α and toxicity against four reference cell lines have been determined for eight compounds identified as the most promising in nanoDSF assay.


Introduction
Protein kinase CK2 is considered as one of the essential human proteins, which is involved in numerous cellular processes. This is a highly pleiotropic kinase, which can phosphorylate hundreds of protein substrates including transcriptional factors, proteins affecting the structure of DNA/ RNA, and/or involved in RNA synthesis and translation, signaling proteins and also pathogenic virus replication machinery proteins [1]. The numerous natural substrates of kinase CK2 make it indispensable for many cellular functions including cell viability, apoptosis, cell proliferation and survival, angiogenesis, DNA damage and repair, the ER-stress response, the regulation of carbohydrate metabolism, and in the nervous system [2]. Kinase CK2 is the holoenzyme that acts in vivo as a heterotetramer, consisting of two catalytic (α-and/or α') and two regulatory (β) subunits, but also may act as a monomeric kinase alone, the regulatory subunit only improving the activity [3].
Since kinase CK2 was found to be overexpressed in many forms of human cancers [4], its role in cancer pathogenesis is thoroughly investigated. CK2 could promote cell survival by phosphorylation of sites adjacent to caspase cleavage sites, which blocks its activity [4][5][6][7]. Kinase CK2 is also involved in DNA damage response and DNA repair pathways, which may additionally contribute to the regulation of cancer cell survival [8]. Despite the fact that the exact function of hCK2 in oncogenesis remains unknown and also its contribution pathways, which may additionally contribute to the regulation of cancer cell survival [8]. Despite the fact that the exact function of hCK2 in oncogenesis remains unknown and also its contribution may be different among various carcinomas [9][10][11], widespread interest has been expressed in targeting this kinase therapeutically [12].
Our long-term systematic studies on protein kinase CK2 ligands have demonstrated that halogenated benzotriazole derivatives showed moderate inhibitory activity against CK2, which is closely related to the physicochemical properties of the ligand [31][32][33][34][35][36]. We have also analyzed the effect of substitution pattern, demonstrating that benzotriazole derivatives carrying halogen atoms at positions 5 and 6 interact with CK2 much stronger than those substituted at positions proximal to the triazole ring (positions 4 and 7) [34,35]. What is more, the presence of the iodine atom in the structure significantly enhances the affinity [37], which is in line with results obtained by Janeczko et al. for massively iodinated benzimidazole derivatives [16]. Both benzotriazoles and benzimidazoles possess o-phenylenediamine fragment in their structure (see Figure 1). Looking for the novel groups of compounds as putative kinase CK2 inhibitors, we decided to synthesize groups of similar compounds, based on o-phenylenediamine configuration.
Here, we synthesized and described six series of heterocyclic halogenated ligands derived from 4,5-dihalogeno-benzene-1,2-diamines. These ligands were selected to enable the separation of the scaffold effect from the hydrophobic interactions attributed directly to the halogen atoms. We measured the hydrophobicity using reverse-phase HPLC. The binding affinities to the catalytic subunit of human protein kinase CK2 (hCK2α) were initially assessed using low-volume differential scanning fluorimetry (nanoDSF), and the dominating modes of binding were analyzed using a molecular modeling approach. Cytotoxicity against four reference cell lines was tested for the most promising compounds along with their inhibitory activity against CK2.

Hydrophobicity Data
RP-HPLC method, a fast and efficient alternative for standard logD determination in octan-1-ol/water system [42], was used to assess the hydrophobicity of all synthesized compounds. This method has been successfully used in our laboratory for halogenated benzotriazole derivatives, including the evaluation of the effect of hydrophobicity on their binding to the catalytic subunit of protein kinase CK2 [24,43].
As expected, the substitution of two hydrogen atoms with the larger halogen atoms: fluorine, chlorine, bromine, or iodine, respectively, always results in an increase of apparent hydrophobicity (Figure 2A). This trend is generally independent on the scaffold of the molecule ( Figure 2B).

Hydrophobicity Data
RP-HPLC method, a fast and efficient alternative for standard logD determination in octan-1-ol/water system [42], was used to assess the hydrophobicity of all synthesized compounds. This method has been successfully used in our laboratory for halogenated benzotriazole derivatives, including the evaluation of the effect of hydrophobicity on their binding to the catalytic subunit of protein kinase CK2 [24,43].
As expected, the substitution of two hydrogen atoms with the larger halogen atoms: fluorine, chlorine, bromine, or iodine, respectively, always results in an increase of apparent hydrophobicity ( Figure 2A). This trend is generally independent on the scaffold of the molecule ( Figure 2B).

Hydrophobicity Data
RP-HPLC method, a fast and efficient alternative for standard logD determination in octan-1-ol/water system [42], was used to assess the hydrophobicity of all synthesized compounds. This method has been successfully used in our laboratory for halogenated benzotriazole derivatives, including the evaluation of the effect of hydrophobicity on their binding to the catalytic subunit of protein kinase CK2 [24,43].
As expected, the substitution of two hydrogen atoms with the larger halogen atoms: fluorine, chlorine, bromine, or iodine, respectively, always results in an increase of apparent hydrophobicity (Figure 2A). This trend is generally independent on the scaffold of the molecule ( Figure 2B).  Indeed, regardless of the scaffold, the increase of hydrophobicity upon halogen substitution in distant positions of the benzene ring of 2-trifluoromethyl-1H-benzimidazole (3a), 2-hydroxymethyl-1H-benzimidazole (4a), and 1H-benzo[d]imidazol-2(3H)-one (5a) is identical. A similar trend is also observed in the case of quinoxaline-2,3-diol (7a) and benzo-1,2-diamine (1a), however, for the latter, the effect of halogenation is shifted downwards by 0.1 log(τ) units. The largest changes in hydrophobicity experience quinoxaline (6a) derivatives, which should be directly attributed to the lack of the effect of electronegative halogen atom on the proximal proton-donating groups. Such effect, partially compensating halogen hydrophobicity, was evidenced in protein-ligand structures, in which hydrogen bonds proximal to a halogen atom were statistically shortened (so must be regarded as stronger) when a nitrogen atom close to the halogen atom (less than six chemical bonds) was a hydrogen bond donor [44]. No such configuration is present in 6a, so the effect of halogen substitution is larger because it cannot be compensated in 6c-e by the halogen-induced strengthening of solute-solvent interaction. Finally, the effect of halogen substitution in the most hydrophobic 2,1,3-benzothiadiazole (2a) remains irregular, which putatively results from an uncontrolled aggregation of these solutes.

Binding Affinity-Thermal Shift Assay
Binding affinities of the newly synthesized compounds towards hCK2α were assessed with the low-volume differential scanning fluorimetry (nanoDSF). The shift of the temperature of thermal denaturation in the presence of 10-fold excess of a particular ligand relative to the apo form of hCK2α (∆T m ) can be used as a semi-qualitative measure of the binding affinity [35] or inhibitory activity [45].
Groups of ligands with different scaffolds vary visibly in their affinity to the target protein (see Table 1 and Figure 3). Interestingly, two of the most hydrophobic groups of compounds, 2,1,3-benzothiadiazoles (2a-e) and quinoxalines (6a-e), virtually do not interact with hCK2α-for all these compounds, ∆T m value never exceeds 1 • C. Benzene-1,2-diamine derivatives (1a-e) also interact weakly with hCK2α (∆T m value 1.7 • C for 1e), but this effect should be attributed directly to the size of the ligand-the absence of the second aromatic ring in the scaffold that appears to be necessary for the efficient ligand-binding to the kinase [20]. There is no general correlation between hydrophobicity and binding affinity, however, when one analyzes individual series representing the same scaffold. A comparison of log(τ) and ∆T m for the protein-ligand complexes indicates that any increase in ligand hydrophobicity is clearly reflected by an increase in binding affinity, as long as the ligand binds to the target (see Figure 3). Moreover, it seems that the effect of Cl, Br, and I substitution on the stabilization of the protein-ligand complexes not only does not depend on the structure of the heterocyclic ring but is directly proportional to the hydrophobicity of the halogen. The latter is evidenced by almost parallel lines representing the contribution of halogen atoms in individual ligands. It could be thus anticipated that for the studied series of ligands, the thermodynamic contribution of non-specific hydrophobic interactions of halogen atoms to the free energy of ligand binding remains separated from that of the direct electrostatic interactions involving a more polar heterocyclic ring. This hypothesis was tested with in silico modeling of protein-ligand structures. . Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ΔTm determined using nanoDSF.

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct measurement of IC50-solely compounds with the thermal shift above 2 °C may be considered as at least moderate inhibitors of hCK2α.
Comparison of thermal shift assay data (∆Tm) with the directly measured IC50 con- Figure 3. Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ∆T m determined using nanoDSF. and I substitution on the stabilization of the protein-ligand complexes not only does not depend on the structure of the heterocyclic ring but is directly proportional to the hydrophobicity of the halogen. The latter is evidenced by almost parallel lines representing the contribution of halogen atoms in individual ligands. It could be thus anticipated that for the studied series of ligands, the thermodynamic contribution of non-specific hydrophobic interactions of halogen atoms to the free energy of ligand binding remains separated from that of the direct electrostatic interactions involving a more polar heterocyclic ring. This hypothesis was tested with in silico modeling of protein-ligand structures. and I substitution on the stabilization of the protein-ligand complexes not only does not depend on the structure of the heterocyclic ring but is directly proportional to the hydrophobicity of the halogen. The latter is evidenced by almost parallel lines representing the contribution of halogen atoms in individual ligands. It could be thus anticipated that for the studied series of ligands, the thermodynamic contribution of non-specific hydrophobic interactions of halogen atoms to the free energy of ligand binding remains separated from that of the direct electrostatic interactions involving a more polar heterocyclic ring. This hypothesis was tested with in silico modeling of protein-ligand structures.  . Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ΔTm determined using nanoDSF.

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct  3. Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ΔTm determined using nanoDSF.

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct  3. Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ΔTm determined using nanoDSF.

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct  3. Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ΔTm determined using nanoDSF.

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct  3. Correlation between HPLC-derived hydrophobicity index (log(τ)) and the ΔTm determined using nanoDSF.

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct

Inhibitory Activity
Since in the thermal shift assay, most of the studied compounds seem to have a negligible affinity towards hCK2α, only eight compounds were further subjected to the direct measurement of IC 50 -solely compounds with the thermal shift above 2 • C may be considered as at least moderate inhibitors of hCK2α.
Comparison of thermal shift assay data (∆Tm) with the directly measured IC 50 confirms that nanoDSF is an attractive method for inhibitors screening, however, for ligands with different scaffolds, hence different modes of binding, this might not be directly reflected in the biochemical assay.

Cell Toxicity
The same four bromo and four iodo-derivatives that, according to nanoDSF assay, interact with the catalytic subunit of CK2α (scaffolds 3, 4, 5, 7) were tested against the reference cancer cell lines: epidermoid carcinoma A-431, colorectal carcinoma HCT116, HCT116p53 −/− , and normal fibroblasts BJ. All of them displayed moderate activity to reduce cell viability (<100 µM). Interestingly, compound 3e was identified as even more active than 3d, the promising activity of which has been already reported [46]. In this study compound, 3d displayed IC 80 values between 2 and 29.1 µM tested for five cancer cell lines (IM-9, MOLT3, U-937, MCF-7, and PC-3), and 8.7 µM on phytohemagglutinin-stimulated normal human lymphocytes. It was especially active against two cell lines: MCF-7 (breast carcinoma) and PC-3 (prostate carcinoma), with IC 80 values below 3 µM. Our data show the selectivity index of 3e~10 (relative to BJ fibroblasts), which makes this compound worth further investigations. We were unable to determine the cytotoxicity for 5e, for which the observed nonmonotonic dose-response effect was most probably caused by uncontrolled aggregation. All cell toxicity data are summarized in Table 2, and the analyses of viability are shown in Supplementary  Figures S30-S37. The toxicity is generally uncorrelated with the inhibitory activity; however, this effect can be directly assigned to membrane permeability. This hypothesis is strongly supported by a clear correlation between compound hydrophobicity (log (τ)) and IC 50 values estimated from viability data (Figure 4).

Molecular Modeling
Molecular modeling was initially performed to assess whether a given scaffold could be harbored at the ATP binding site. Free energy of ligand binding (ΔGbind) was estimated with the aid of the VINA-AutoDock algorithm implemented in Yasara Structure [47]. The obtained ensembles of structures agreed qualitatively with our knowledge concerning the Interestingly, contrary to the enzymatic assay, the activity of 7e remains marginal (IC 50~1 00 µM), but such a low value should be assigned to a moderate hydrophobicity of this compound. In this context, one can expect that esterification of the two hydroxyl groups of 7e would improve its apparent activity in viability tests.

Molecular Modeling
Molecular modeling was initially performed to assess whether a given scaffold could be harbored at the ATP binding site. Free energy of ligand binding (∆G bind ) was estimated with the aid of the VINA-AutoDock algorithm implemented in Yasara Structure [47]. The obtained ensembles of structures agreed qualitatively with our knowledge concerning the binding of halogenated ligands at the ATP-binding site of CK2α [21,48,49]. Thus, all chlorinated, brominated, and iodinated ligands preferably adopt the orientation, in which both halogen atoms are solvent-protected, being close to the hinge region, while the heterocyclic ring points towards the polar region of the ATP-binding site formed by sidechains of Lys68, Glu81, and Asp175. The competition between hydrophobic and electrostatic interactions is visible for the majority of tested ligands, again clearly identifying the balance of hydrophobic and electrostatic interactions as the main driving force.
Interaction between a fluorinated ligand (1b, 2b, 3b, 7b) and the protein is found generally overestimated ( Figure 5, red points). Contrarily, the binding affinity of the majority of iodinated and some brominated ligands is underestimated ( Figure 5, blue and green points). However, the latter effect may reflect the absence of the parameterization of halogen bonding interactions in the classical force-field used in the scoring function [50]. On the other hand, these discrepancies can be treated as indicators for the existence of halogen bonding in these complexes, which mainly applies to 3e, 4e, and 7e. Interestingly, despite the discrepancy in the estimated binding affinity, most of the low-energy poses of 5,6-diiodo-, 5,6-dibromo, and 5,6-dichloro-derivatives preserve halogen atoms pointing toward the hinge region, thus, in the orientation, enabling halogen-bonding, with the heterocyclic ring close to Lys68 (see Supplementary Table S1 and Figures S40 and S41). with the heterocyclic ring close to Lys68 (see Supplementary Table S1 and Figures S40 and  S41).

Discussion
We have synthesized six groups of halogenated heterocycles, based on benzene-1,2,diamine derivatives. To analyze their applicability as the precursors of potential hCK2 inhibitors, we used molecular modeling together with the nanoDSF method, which ena-

Discussion
We have synthesized six groups of halogenated heterocycles, based on benzene-1,2,diamine derivatives. To analyze their applicability as the precursors of potential hCK2 inhibitors, we used molecular modeling together with the nanoDSF method, which enabled qualitative comparison of their binding affinity. Many studies have been carried out on CK2 Type I inhibitors and some preferences of ligands in CK2 ATP-binding sites are already well known [20,33,35,48,51]. Besides the obvious steric hindrance of the scaffold, the two main driving forces of binding ligand at the ATP-binding site of CK2α are clearly identified: the electrostatic contribution together with hydrophobic effect [20,33,35]. Negatively charged ligands commonly bind deeper in the CK2α pocket, preferably interacting with Lys68 [20]. Additionally, the substitution with hydrophobic halogen atoms significantly increases the affinity towards the CK2α by a favorable change in the physico-chemical properties of the ligand. The latter also includes the possible formation of halogen bonds with the hinge region [48].
Using these different groups of compounds, we analyzed the interplay between the scaffold and the effect of substitution with halogen atoms. For smaller scaffolds, like benzotriazole/benzimidazole, unsubstituted on triazole/imidazole ring, a salt bridge or hydrogen bond formation with Lys68 and halogen bonding with the hinge region are largely mutually exclusive because these two regions are in the hCK2α, too far apart for these interactions to occur simultaneously [20,48]. This is why we tested larger scaffolds with either a larger heterocyclic ring (6a-e, 7a-e) or substitution in the imidazole, both of them to preserve possible interactions with a hinge. We analyzed solely dihalogenated compounds of the same pattern of halogen atoms at the benzene ring to assure that predominating interactions with the hinge region (occurring with more halogen atoms) and configuration effect (depending on the halogenation pattern) will not predominate or interfere with the interdependence we were analyzing [24,34,43]. Interestingly, for the tested halogen configuration, we can notice that only ligands with uneven charge distribution on the non-benzene ring (groups 3, 4, 5, 7) exhibit measurable affinity towards CK2, while ligands even much hydrophobic but without proton-donating groups (series 2 and 6) do not stabilize the complex (i.e., ∆T m < 1 • C). This observation indicates that two halogen atoms are insufficient to stabilize the ligand in the complex with a position close to the linker region, and the formation of hydrogen bonds with the lysine seems to be necessary. Such dependence of ligand position on the number of halogen atoms was observed for the halogenated benzotriazoles, where 5,6 dibromobenzotriazole interacts with Lys68 (pdb 6TLP), while tetrabromobenzotriazole interacts either with the hinge region or with Lys68 (pdb 6TLL). What is more, an inspection of in silico determined structures for complexes with ligands that bind to the protein (∆T m < 1 • C) points out to this unique interplay, as ligands with no halogen atom or with fluorine (which is not capable for the formation of halogen bond [52]) display in almost all cases different orientation relative to the hinge region (see Figures S40 and S41), possibly forming hydrogen bonds. The substitution with other halogen atoms switches the preferable orientation of the ligand, which has an effect that is somehow explicable by hydrophobic interactions. However, since the force-field used in the docking procedure does not include the explicit potential for a halogen bond, these interactions are underestimated.
These different orientations of chlorinated, brominated, and iodinated compounds might explain why the relation of measured log(τ) and ∆T m values for series representing the same compound scaffold is linear solely for these three substituents, but not for hydrogen and fluorine (see Figure 3).
It is worth noting that for these two 5,6-halogenated scaffolds, the observed tendency is opposite to that observed for their tetrahalogenated analogs. Thus, 2-trifluoromethyl-1H-benzimidazole interacts stronger than the 1H-benzo[d]imidazol-2(3H)-one counterpart (7 vs. 13 µM). This again points out the subtle interplay between the scaffold and the substitution with halogen atoms, indicating that the interaction of 2-trifluoromethyl-1Hbenzimidazole with Lys68 is slightly stronger than for its isostructural analog-carrying carbonyl group.
Finally, the viability tests performed for reference cell lines showed moderate activity for some compounds. Interestingly, IC 50 values determined in this test are correlated with a compound's hydrophobicity rather than the inhibitory activity in the CK2α assay. This indicates a dominant role of membrane permeability, but it may also indicate that CK2α is not a major target for studied test compounds.

Chemistry
All starting materials and solvents for reactions were purchased from Sigma Aldrich (now Merck KGaA, Darmstadt, Germany), Fluorochem (Hadfield, UK), ABCR (Karlsruhe, Germany), or Chempur (PiekaryŚląskie, Poland). 1 H-NMR spectra were recorded with Varian 500 MHz spectrometer. TMS or the residual solvent signal were used as the internal standard. High-resolution mass spectrometry spectra were recorded on a TQ OrbitrapVelos instrument (Thermo Scientific, Waltham, MA, USA). The reaction progress was monitored by the thin-layer chromatography analysis using silica gel plates (Kieselgel 60F 254 . Merck, Darmstadt, Germany). Column chromatography was performed on Silica Gel 60 (0.040-0.063 mm. Merck, Darmstadt, Germany).
Detailed experimental data, as well as spectral data for all synthesized compounds, are available in the Supplementary Materials.

Synthesis of 2,1,3-benzothiadiazole derivatives (2a-e)
Benzene-1,2-diamine derivative (5 mmol) was dissolved in dry DCM (15 mL). Triethylamine (3 mL) was added, then thionyl chloride (15 mmol, 1.1 mL) dissolved in 3 mL of dry DCM was added dropwise at 0 • C. The reaction mixture was allowed to reach room temperature, then it was heated at 40 • C overnight. Afterward, the reaction mixture was filtered. The filtrate was evaporated, the crude product was purified by column chromatography on silica gel using DCM as eluent. Synthesis of 2-trifluoromethylo-1H-benzimidazole derivatives (3a-e) Benzene-1,2-diamine derivative (2 mmol) was dissolved in trifluoroacetic acid (2 mL) and a catalytic amount of concentrated HCl was added. The reaction mixture was heated at reflux overnight. The reaction was quenched by the addition of 50 mL of concentrated NaHCO 3 solution; afterwards, it was extracted by ethyl acetate (3 × 30 mL). Combined organic layers were dried by MgSO 4 , then the solvent was evaporated. The crude product was purified by column chromatography on silica gel using 85:15 (hexane:AcOEt) as eluent (compounds 3d and 3e) or by recrystallization from hexane (compounds 3a, 3b, and 3c).

Synthesis of 2-hydroxymethylo-1H-benzimidazole derivatives (4a,c-e)
Benzene-1,2-diamine derivative (5 mmol) and hydroxyacetic acid (15 mmol, 1.14 g) were dissolved in water (1.5 mL) and concentrated HCl (0.5 mL) was added. The reaction mixture was heated at reflux overnight. The reaction was cooled to RT and then 20% NaOH solution was added until pH = 13. The formed participate was filtered and washed several times with water. Compound 4e was additionally purified by recrystallization from MeOH.

Molecular Modeling
Docking of all tested ligands at the binding site of the catalytic subunit of human protein kinase CK2 was performed with VINA algorithm implemented in Yasara Structure package (www.yasara.com, version 20.12.24) using amber14 force-field with 8 Å cutoff for long-range electrostatic interactions. We used 8 just-resolved structures of complexes of hCK2α with variously brominated benzotriazoles (Protein Data Bank records 6TLW, 6TLV, 6TLU, 6TLS, 6TLR, 6TLP, 6TLO, and 6TLL) as the reference protein structures. This, including alternative protein conformations, stated 18 structures to be tested as the 'receptor'. Such an extended set of template protein structures was chosen to allow the sampling of a larger conformational space that allows multiple ligand orientations at the binding site (see Supplementary Figure S39). In silico screening procedure was restricted to the cuboid region that contained all residues from the ATP binding site, further additionally extended by 2 Å in each direction. During simulations, the sidechain atoms of protein residues proximal to the location of the original ligand (4 Å threshold) were flexible, while coordinates of all other protein atoms were kept fixed. Finally, each of 35 ligands was subjected to 25 independent cycles of the docking procedure using 18 structural templates. The resulting structures were then clustered, and for each ligand, the highest scored one was used as the representative. This score was used as a rough estimate of the free energy of ligand binding.