Identification of novel thiadiazin derivatives as potentially selective inhibitors towards trypanothione reductase from Trypanosoma cruzi by molecular docking using the numerical index poses ratio Pr and the binding mode analysis

Chagas disease is a serious health problem in Central and South America for which effective treatment is not currently available. This illness is caused by the protozoa Trypanosoma cruzi, a species that relies on a thiol-based metabolism to regulate oxidative stress. Trypanothione reductase enzyme plays a central role in the metabolic pathway of the parasite. In this work, a virtual screening of a library of novel thiadiazine derivatives against trypanothione reductase using molecular docking was performed. Four different series of hybrid ligands having in the structure one or two peptoid moieties (series I and II) or the tetrazole ring (series III and IV) were considered. An ad hoc numerical index called poses ratio was introduced to interpret the results of the docking analysis and to establish relevant structure-interaction relationships. In addition, six binding modes were found for the ligands with the highest populated conformational clusters after applying contact-based analysis. The most regular and relevant were binding modes I and II, found mainly for ligands from series I. A subsequent molecular docking on human glutathione reductase enzyme allowed to assess the possible cytotoxicity of the ligands towards human cells. A selective binding profile was found for ligands with interactions in the Hydrophobic cleft, the spermidine and the Z subsites inside the active site of trypanothione reductase. At the end of the study, new thiadiazine-based compounds were identified as plausible candidates to selectively inhibit the parasitic enzyme.


Introduction
Chagas disease or American trypanosomiasis is still considered a serious health problem in Central and South America, according the Worldwide Health Organization (WHO), despite it was discovered more than a century ago. Furthermore, imported Chagas disease is increasingly viewed as a global emergent health issue, due to the migration of infected individuals to other continents, mainly North America and Europe. Around six to seven millions of people worldwide are affected with this illness, with 12 000 fatalities reported each year and around 65 million of human beings are in risk of contracting it. Chagas disease is caused by the protozoa agent Trypanosoma cruzi (T. cruzi) and has vector-borne transmission and nonvectorial transmission mechanisms [1][2][3]. In the absence of effective anti-trypanosomal treatment the sickness has acute and chronic phases. The acute phase, that only lasts 8-12 weeks, is characterized by circulating trypomastigotes in blood; while in the chronic phase parasitemia levels become undetected by microscopy. Around 30% of the chronic patients develop, over a period of years to decades, a clinically evident cardiac and/or gastrointestinal disease, probably related to the presence of T. cruzi in smooth muscle and autonomic ganglia in heart, esophagus or colon [1][2][3].  [9] from T. cruzi and glutathione reductase (right; PDB:1GRA) [10] from H. sapiens. Both proteins are depicted in the same scale and the surface of the corresponding active sites are highlighted. b Structures of trypan-othione disulfide (TS 2 ) and glutathione dimer (GSSG), the natural substrates for TR and GR, espectively. The positively charged spermidine (Spm) residue, distinctive of TS 2 is colored red. c Biochemical reactions catalyzed by TR and GR. NADP + : nicotinamide adenine dinucleotide phosphate and its reduced form NADPH C-terminus by a basic spermidine (Spm) chain. Because the spermidine component is not symmetrical, the molecule is not really a dimer and the two tripeptide portions are distinguishably named I or II. Hence, trypanothione structure is often abbreviated γ-GluI-CysI-GlyI-Spm-GlyII-CysII-γ-GluII. Mammals have an equivalent enzyme designated as glutathione reductase (GR), which exerts a similar catalytic function than TR (Fig. 1a). GR reduces glutathione dimer (GSSG) to its thiol monomeric form L-γ-glutamyl-L-cysteinyl-glycine (GSH). GSSG and GSH from mammals are the molecular analogs of TS 2 and T(SH) 2 from T. cruzi, respectively (Fig. 1b). GSH scavenges potentially dangerous oxidizing reagents which are the result of aerobic metabolism within the cell (Fig. 1c) [8]. Because GSSG lacks the Spm residue, its overall charge is negative, in contrast with the larger and positively charged TS 2 .
Enzymes GR and TR belong to the family of flavincontaining disulfide oxidoreductases. They are active as homo-dimers of about 52 kDa per subunit and have an overall 30% of sequence identity, with 14 conserved residues in the substrate binding site [8,10]. In both molecules, the catalytic site is rather complex, comprising two separate regions, i.e., the NADP site (N-site) and the active site (G-site), connected to each other by the flavin ring of flavin adenine dinucleotide (FAD) and a redox active disulfide bridge [8,11]. Despite all these similarities, the specific amino acid alterations among both enzymes affect the size and the chemistry of the disulfide-substrate-binding sites. In the literature, the difference in the substrate specificity of TR and GR has been explained in terms of steric and electrostatic factors [11]. The active site of TR is rich in glutamic acid and non-polar residues; while in GR basic residues (such as arginine) prevail. The orientations of the domains in the enzymes are also different, with the consequence that the disulfide-substrate-binding sites don't have the same size. Therefore, the parasitic TR binds the larger positively charged TS 2 , while human GR binds smaller double negatively charged GSSG (Fig. 1b). The fact that, in general, the enzymes bind different substrates, has been the cornerstone of many studies pursuing selective compounds against Chagas disease [11].
Thiadiazines, more specifically 3,5-disubstituted thiadiazine-2-thione (THTT), have been studied by the Laboratory of Organic Synthesis (LSO) for more than 20 years [12][13][14][15][16]. Several biological studies in vitro have demonstrated the antiparasitic activity of THTT, which is related to a non-specific high cytotoxicity [12,13]. In the last years, our group has focused on obtaining molecules containing the thiadiazine core with higher activity and/ or lower cytotoxicity. Although the exact mechanism of action of thiadiazines remains obscure, the possibility of a selective inhibition of TR by some of thiadiazine-based compounds might reduce their inherent cytotoxicity. With that purpose, multicomponent procedures such as the Ugi-four-components reaction have been used for diversity-oriented synthesis of hybrid molecules with improved properties [15]. This versatile reaction enables the modification of molecules by introducing a dipeptoid residue [16] or a tetrazole ring and additional substituents-the diversity points-in the molecular scaffolds, which ease the design of libraries of new compounds. Similar synthetic approaches has been extensively used for the design of thousands of biologically active molecules [17]. In particular, tetrazoles have found widespread use in the context of medicinal chemistry because of their simple preparation, bioisosteric properties and high chemo and enzymatic stability [18].

Fig. 2
Chemical structure of the four thiadiazine-based series of ligands designed for molecular docking. The distinguishing molecular scaffolds appear colored: the thiadizine ring is bur-gundy, the N-substituted dipeptide is green and the tetrazole is ultramarine.*In this position, both possible stereoisomers were considered As a continuation of our previous studies [15,16], we designed a virtual library of 2094 hybrid THTT (Fig. 2) to assess their potential role as TR inhibitors through an in silico approach. The chemical structures of all the substituents (see Supplementary Material SA1) within the ligands were chosen so that any compound in this analysis may be easily synthesized from readily available commercial chemicals, including natural or synthetic mimetics of amino acids [15,16,[18][19][20]. In order to select the best candidates for a future experimental study, the compounds were screened using a filtering protocol based on molecular docking. It is widely known that this computational tool effectively predicts the affinity and binding mode of a given ligand towards a particular macromolecular target [21,22]. Molecular docking has been previously reported as a useful tool in the context of searching potential drugs against T. cruzi using TR, GAPDH, trans-sialidase, polyamine transport and cruzain enzymes as receptors [6,7,9,11,[23][24][25][26][27][28][29][30][31][32][33][34][35]. In the present theoretical study, TR was selected as target taking into account the decisive role of this enzyme in T. cruzi metabolism. Therefore, the ligands were designed to be structurally similar to known inhibitors of TR-having a heterocyclic and peptidic nature-. Ligands with the highest potential to act as TR inhibitors were included in an additional docking simulation towards GR to analyze their potential cytotoxicity towards human cells and therefore their profile as selective inhibitors. As the examination of molecular docking results is often cumbersome because of the large amount of ligands considered and different enzyme-ligand interaction patterns, we have introduced a novel parameter called Pose ratio (Pr). This is a numerical index that eases the interpretation of the clustering by assessing the propensity of an individual ligand towards a specific binding mode and the effect of a particular substituent in the overall binding with the enzyme.

Preparation of the ligands and TR protein
The designed library of potential ligands comprises four series or families of new hybrid thiadiazine-based structures. Series I and series II contain one or two peptoid moieties, respectively; while series III and IV have a tetrazole ring and differ mainly in the position of this heterocycle. Polar, charged, aliphatic and aromatic substituents were included in the structures. Neutral media (pH = 7.0) was considered to set the protonation states of the ligands, i.e. protonated-positively charged-amines, guanidine and imidazoles; and deprotonated-negatively chargedcarboxylic acids and 1H-tetrazoles. The three-dimensional structure coordinates were generated using Avogadro 1.1.1 [36] software and were pre-optimized by molecular mechanics with the Steepest Descent Algorithm and force-field MMFF94 [37]. Further optimization with PM6-DH2 semi-empirical Hamiltonian [38,39] implemented in MOPAC package [40] followed. All ligand PDB files were converted into PDBQT format, the input of AutoDock Vina program [41], in the guest user interphase AutoDockTools [42]. Charges on the ligand atoms were generated using the Gasteiger model, non-polar hydrogens were merged and default active torsions (rotatable bonds) were retained using the TORSDOF utility (Supplementary Material SA0).
The crystal structure of the protein TR (PDB: 1BZL) [9] and GR (PDB: 1GRA) [10] resolved at 2.40 Å and 2.00 Å, respectively, were retrieved from the Protein Data Bank (PDB, http://www.rcsb.org). BIOMT transformations were applied to generate the coordinates for the multimer representing known biologically significant oligomerization state of the enzyme 1GRA. These structures were optimized using pdb2pqr.py (Version 2.1.0) online server [43] with AMBER force field [44] and the protonation states of ionizable groups at pH = 7.0 were assigned by using PROPKA [45]. The ligands and waters were eliminated and the PDBQT files were generated with AutoDockTools [42]. The size of the simulation grid was selected at 36 × 44 × 44 Å 3 and 26 × 30 × 28 Å 3 for TR and GR, respectively. The center of the simulation box was set at the geometrical center of the active site.

Docking procedure and analysis of enzyme-ligand complexes
For the validation of the docking run, previously randomized structures of the natural substrates of the enzymes (TS 2 and GSSG) were submitted to re-docking experiments. The obtained poses were compared with the experimental geometries in the enzyme-ligand complexes through the total Root-Mean-Square-Deviation (RMSD). A contact-based analysis allowed to match the enzymeligand interactions present in the obtained theoretical poses with those present in the crystal pose. Docking simulations were performed using AutoDock Vina 1.1.2 program [41], which will be referred as Vina in the upcoming explanations. The active sites of TR and GR were kept rigid during the simulation. The docking parameters were set to default with the exception of the following: exhaustiveness = 32 and num_modes = 2. The Vina code predicts the adopted conformation of an enzymeligand complex and organizes the solutions based on a scoring function, which correlates the obtained conformations with the binding affinity (kcal/mol). The best 2 docked conformations or poses (num_modes = 2) from 10 independent runs were analyzed, to produce 20 final docked poses per ligand. These poses were then clustered based upon RMSD using a cutoff of 2.5 Å with an ad hoc Python (Version 3.7) [46] script that implements Daura's clustering algorithm [47].
The binding energy of each cluster was set as the average binding energy of all the poses. For a particular ligand, the three clusters with the highest number of poses were selected as its binding modes; the remaining poses were considered as non-clustered. A binding mode described by a cluster with 10 or more poses was considered as a representative binding mode. This threshold, that includes at least half of the poses generated by the program for each ligand, was selected to be in agreement with the criterion of Rosenfeld, Olson and co-workers [48]. For several ligands, two representative binding modes were found, i.e. two clusters with 10 poses each.
Having the clusters for each binding mode, the parameter P r (as from Poses ratio) was defined as the ratio between the number of representative poses (those forming a representative binding mode) and the number of non-representative poses (any other pose, whether clustered or not) for ligands having a common substituent in their structure. P r values were calculated for all the substituents, while the mean P r and the corresponding standard deviation were calculated for each position (see Supplementary Material SA1) within each series. Supplementary Material SD contains an Excel spreadsheet exemplifying the calculation of P r for the substituents in compounds from Series IV.
Further analysis of interatomic contacts between the ligands with representative binding modes and the enzymes TR and GR was performed with the freely available Python-implemented computer algorithm BINANA [49]. Key binding characteristics such as hydrogen bonds, hydrophobic contacts, salt bridges (electrostatic attractions among opposed charged fragments) and π interactions were identified. All adjustable parameters by BINANA were set as default (see Supplementary Material SA0). When a ligand showed contacts with at least 60% of the amino acid residues forming a subsite, that interaction was considered as important. Ligands were regarded as potential inhibitors of TR and GR enzymes when they showed interactions with at least six (out of ten) or four (out of six) subsites, respectively. Ligands with two representative binding modes were considered as potential inhibitors of an enzyme if any of the two analyzed poses fulfilled the previous requirement. Figures depicting the binding modes of the ligands were prepared using PyMOL™ Version 2.3.1 [50].

Results and discussion
In order to evaluate the accuracy of the selected docking parameters, the natural substrates TS 2 and GSSG were docked in the active sites of TR and GR, respectively. Two clusters with favorable energy were obtained for both natural substrates ( Table 1). The poses showed significant differences in position (RMSD > 2 Å) when compared to the experimental conformations of the co-crystallized ligands. However, contact-based analysis using BINANA showed that more than half of the interactions found in the crystallized protein complexes obtained from the PDB were reproduced. These results are in agreement with similar studies in the literature [25,[31][32][33][34][35]. The fairly good coincidence between the binding mode of TS 2 and GSSG in the two clusters encouraged us to retain the same docking parameters for the subsequent simulation runs.

Unravelling the role of the substituents on the binding
According to Rosenfeld-Olson, binders may be distinguished from non-binders, to some extent, using a combination of two trends, namely, predicted binding energy and level of clustering [48]. Conversely, it is known that the energy results derived from docking programs must be treated with caution because values rely on the size (number and type of atoms) of each particular molecule. Thus, absolute energy values are somewhat meaningless, unless they are supported by experimental data, which is not the case at this time [33]. The designed hybrid thiadiazine-based ligands were docked into the active site of TR using Vina program [41]. The binding energy of each cluster was calculated and negative values were obtained in all cases. Insofar the Hence, the distinction among potential inhibitors and non-inhibitors was determined taking into account the poses population of the clusters and the enzyme-ligand interactions as described below. Table 2 summarizes the results of the clustering following the docking towards TR enzyme. Increasing the number of active torsions in the ligands scaffolds rises the amount of non-clustered poses and less populated (non-representative) binding modes are achieved. Thus, the percentage of clustered poses and amount of representative binding modes followed the order Series II < < Series I < Series III < Series IV. Evidently, ligands from series I, III and IV have higher potentialities as TR binders. In the case of ligands with many active torsions (Series II), non-representative binding modes were mostly obtained, except for two cases. In total, 858 ligands from the whole library showed non-representative binding modes, thus having less potential for any inhibitory activity. Consequently, these compounds were ignored in the following examination.
Furthermore, aiming to get insights into important structure-interaction relationships, the parameter Poses ratio (P r ) was introduced (see the section Material and Methods for the relevant descriptions). Figure 3 illustrates the P r values calculated for all the substituents of the ligands from series I, III and IV, as well as the mean and  . 3 The influence of the substituents on cluster formation. The structures of all the substituents (from R 1 to R 11 in X axis) can be found in Supplementary Material SA1. Substituents with a positive ( +) or negative ( −) statistically significant difference with mean values are highlighted. *For the substituent in this position, the values of P r encompass the two possible stereoisomers standard deviation found in each position. Similar P r mean values were obtained for different positions in the same series. Nevertheless, high dispersion (standard deviation) of the values was found in all cases, which it is associated to the different nature of the substituents and their particular structural requirements to bind in a well-defined (representative) binding mode to the enzyme target TR. As a general trend, positive statistical significance on P r values were found for voluminous substituents with no active torsions, where the steric factor (volume) is determinant. These substituents increase the enzyme-ligand interactions while only adds few active torsions to the structure. Negative statistical significance was obtained for long chains, groups with several routable bonds. In Series I, the largest P r was observed for compounds with a tert-butyl (4a) in position R 4 . Large P r were also encountered for compounds with the substituents 2e, 2m, 2o and 3g which correspond to rather diverse amino acid side chains of Glu, Trp, His and D-Phe, respectively. Low P r were obtained for 1b (isobutyl), 3b (isobutyl) and 3f (CH 2 CO 2 CH 3 ); whereas the lowest values corresponded to the highly flexible pentylene and propylene aliphatic chains of 2b and 2c, respectively. In the case of ligands from Series III, only those with very bulky substituents were favored. The largest P r were obtained for substituents 7a (cyclohexyl), 8b (benzyl) and 9e (benzyl). In contrast with the previous case, the presence of the tert-butyl (7b) at the R 7 position did not exert a statistically significant difference. Noteworthy, the absence of a substituent at positions R 8 and R 9 (8a, 9a = H) substantially diminished the formation of representative binding modes, even when there is no variation in the overall amount of active torsions of the structures. This last particular example illustrates the diversity of factors influencing P r , such as the strength of the enzyme-ligand interactions. Regarding compounds from series IV, ligands with a cyclohexyl substituent (10a) in position R 10 were favored over other possibilities. Reasonably, this ring has the ability to pack and interact strongly with the residues within the catalytic site. In position R 11 the most suitable substituent was the isopropyl (11c) group.

Identifying the potential inhibitors via contact-based analysis
Crystallographic and molecular modelling data has allowed to get insights into the interaction between TR enzyme and its natural substrate TS 2 [9,27,28,30]. Table 3 contains the amino acid residues from TR implicated in the binding with TS 2 -that is, the contacts-via hydrogen bonding, van der Waals and other non-covalent interactions that have been identified [4]. These amino acid residues have been grouped in seven subsites, based on their interaction distance restriction of 4.0 Å with each building block from the natural substrate trypanothione, namely, γ-GluI, CysI, GlyI, Spm, γ-GluII, CysII and GlyII (Fig. 4) [9]. Other subsites (Z, Y and Hydrophobic cleft) relevant for the interaction have also been identified within TR active site (Fig. 4b) [9,11,26]. These subsites have been previously used as a template for contact-based analysis in the rational design of potential TR compounds by molecular docking. For example, analogs of clomipramine bearing a tricyclic dibenzosuberyl moiety that showed antitrypanocidal activities in vitro were recognized as competitive inhibitors of TR [24]. Docking studies identified five interaction subsites (γ-GluI, Z, CysI, GlyI and Spm sites) as the binding locations of these tricyclic inhibitors. Likewise, a comprehensive analysis of the binding modes of novel triterpenic chemotypes from the Natural Products Database (NatProDB) from Bahia State allowed the identification of Table 3 Relevant subsites of interactions within the active sites of enzymes TR and GR and their composition.
A and B stand for the names of independent polypeptide chains from the enzyme homodimers Investigation that relies exclusively on the predictions of the scoring function provided by the docking software is not able to accurately find all active compounds. Therefore, to discriminate potentially active from non-active ligands, one model pose from the 1355 representative binding modes of the 1236 ligands (see Table 2) obtained after the clustering were considered for a contact-based examination. The purposes of this step were: (1) to make a further filtering of the ligands, (2) to understand the kind of interactions present in the enzyme-ligand complexes and (3) to identify potential inhibitors of TR. First, the contacts of the ligands with the 28 amino acid residues constituting the 10 subsites within TR active site were obtained using BINANA (see Supplementary Material SA2 and SA4). Interactions consisted mainly in van der Waals contacts, though several hydrogen bonds, π-interactions and salt bridges were also obtained (Supplementary Material SA4). Then, the potentially active ligands were identified (see Computational Methods section for the criteria used). It was found that 277 ligands fairly mimic the interaction pattern of the natural substrate TS 2 , and also include contacts in the distinctive binding subsites of known inhibitors (Hydrophobic cleft, Z and Y). Taking a general picture, compounds from series I and III proved Fig. 4 Interaction subsites within the active site of TR enzyme. a Active site of TR in complex with trypanothione disulfide (PDB:1BZL). The structure of TS 2 is depicted with colored building blocks. b Different subsites of interaction within the disulfide-binding site of TR. Each subsite-building block couple are colored the same the more suitable binders. Though ligands from Series IV rendered highly populated clusters, their interaction profile was mostly poor, which could be due to their small size. It was also reinforced the idea that neutral bulky substituents with few active torsions in the ligands scaffold contributed more effectively to achieve favorable enzymeligand interactions.
For a better understanding of the binding of the novel THTT scaffolds to TR enzyme, the 1355 poses-from the 1236 ligands obtained from the clustering-were sorted into six average binding modes (BM TR ), based on similarities in the interaction pattern found in enzyme-ligands complexes. Table 4 summarizes the important interactions of the BM TR with the subsites and the most recurrent amino acid contacts. Figure 5 depicts the orientation of the ligands forming each of the BM TR inside the active site. Binding mode I (BM TR I) comprises ligands that showed interactions in the Hydrophobic cleft and Spm subsites of TR, the vast majority from Series I. Hence, BM I is located in the left corner of the TR active site and involves also -GluI, GlyI and CysI subsites, where most ligands also showed contacts. Ligands with binding mode II (BM TR II) showed interactions in the Z subsite of TR. This subsite has proven of great importance in the design of selective inhibitors over GR human enzyme [26,51]. Contacts in the -GluII, -GluI, CysI and GlyI subsites were also achieved. Therefore, BM TR II encompass a large region in the middle and bottom sections of TR active site. Binding mode III (BM TR III) is composed by ligands from series I and series III. It is located in the middle top of TR active site. These ligands displayed contacts in the region I of TR natural substrate site (encompassing -GluI, CysI and GlyI), as well as in the Spm and Y subsites. Binding mode IV (BM TR IV) encloses ligands from series I and series III and is located in the Table 4 Results of binding mode analysis with enzyme TR as receptor. For each binding mode (column 1) the relevant subsites of the protein (column 2) and the recurrent contacts (columns 3-5) are given. Contacts were classified in three groups according to the percentage of the total of poses in which they appear. 17

Inspecting selectivity against human glutathione reductase enzyme
The design of selective inhibitors of TR over GR has been main topic of several molecular docking-based studies aiming to identify potential non-cytotoxic lead compounds for the treatment of Chagas disease. For example, nitrofurans and phenothiazine derivatives have been previously identified as non-selective inhibitors of TR based mostly on the binding energy obtained from molecular docking simulations for the ligands complexes with both enzymes [32][33][34]. It has also been found in enzymatic  studies that the selective inhibitors bind at Z and Hydrophobic cleft subsites from TR active site [11,26,29]. As in the analogous TS 2 -TR system, crystallographic and molecular modelling data has allowed to identify the amino acid residues involved in the binding of GSSG with the active site of GR (Fig. 6) [32][33][34]. This information was applied in the present study for the identification of potential selective inhibitors. A further docking study towards GR was carried out in order to recognize selectivity in the 277 potential inhibitors of TR (Supplementary Material SB2). By this approach, less cytotoxic compounds with a superior profile as chemotherapeutics for Chagas disease might be found. In a similar fashion to the preceding study, Rosenfeld-Olson criterion was applied to discriminate binders from nonbinders [48]. As shown in Table 5, 74 ligands presented non-representative binding modes and were considered automatically non-binders of GR. On the contrary, 203 ligands had representative binding modes, which makes them binders. As binding does not necessary implies inhibition, a second contact-based analysis followed, applying the same approach used previously. The contacts of the ligands with the amino acid residues constituting GR active site were obtained using BINANA (see Supplementary Material SA3 and SA5).Then, the potentially active ligands were identified (see Computational Methods section for the criteria used). On the other hand, the poses from 29 ligands showed interactions with less than three (out of six) subsites. Consequently, these compounds were also regarded as potentially non-inhibitors of GR. Summarizing, at the end of the study, a total of 103 ligands were selected as potential selective inhibitors of TR enzyme (see Supplementary Material SC). These compounds bind strongly to TR active site through representative binding modes, whereas anchor at GR active site with non-representative binding modes (74 examples) or by interacting in a minor amount of the natural substrate subsites (29 examples).    I  21  0  II  0  III  3  IV  0 We were interested in the interaction pattern of the THTT ligands with the six natural subsites in GR enzyme (see Supplementary Material SA3 and SA5). Two average binding modes (BM GR ) were identified (Table 6 and Fig. 7). BM GR I included ligands with contacts in CysII and -GluI subsites, while ligands from BM GR II displayed contacts in the CysII and -GluII subsites. Regardless, both BM GR had additional recurrent interactions with other subsites. Many hydrogen bonds were found, commonly with R37(A), Y114(A) and R347(A) residues; and a few salt bridges were achieved, almost exclusively with E472(B) residues (see Supplementary Material SA5). Additionally, 24 ligands could not be clustered in any defined binding mode and were labeled as BM GR III.
Among the 103 selective TR inhibitors found, 95 ligands belonged to series I, 1 ligand belonged to series II and 7 ligands belonged series III. Their more important binding modes were BM TR I (28 ligands) and BM TR II (59 ligands), which have interactions in -GluI, CysI, GlyI and -GluII subsites corresponding to the natural substrate. This binding modes also display important contacts in the Hydrophobic cleft, the spermidine and Z subsites (see Table 4). These later subsites do not have analogous in the active site of GR and therefore are essential for the selectivity in the binding of the ligands [11,26,29]. Table 7 contains some example of ligands from Series I that do not form any clusters with GR. The substituents present in their structure are neutral or positive charged, which accounts for the selectivity towards the negatively charged binding site of TR, rich in non-polar and acidic residues, as mentioned above. Positions R 1 and R 4 are occupied by a variety of substituents, which are all neutral and bulky. Position R 2 admits a neutral voluminous or a basic positively charged group (guanidino or amino), with the notable exception of I-1521, that carries a hydrogen atoms in this place.
The specific contacts of three ligands (I-0103, I-0332 and I-1001) with the active site of trypnothione reductase are represented in Fig. 8. Though the compounds are bound to TR mostly by van der Waals forces, hydrogen bonds, -stackings, cation-interactions and salt bridges are also present. Compounds I-0103 and I-0332 accommodate in a very similar manner, as expected from their high structural similarity. The benzyl ring in position R 1 of these ligands protrudes towards the Z subsite, affording contacts with P462(B) and L399(B) in the case of I-0103, or with P398(B) and F396(B) in the case of I-0332. This conformation is further stabilized by non-classic interactions involving the benzene ring: I-0103 has a cation-contact with H461(B), whereas I-0332 has a --stacking with F396(B). The thiadiazine ring lays out in the middle, just above the -Glu II subsite and close to V54(A) and H461(B). The basic substituent in position R 2 (3-guanidinopropyl in I-0103 and 4-aminobutyl in I-0332) penetrates to the Gly I subsite, the more flexible 4-aminobutyl chain reaching deeper, down to L18(A) and Y111(A). The basic guanidine and amine groups form hydrogen bonds with G50(A) and S15(A), Fig. 7 Representation of the average binding modes within the active site of GR enzyme. All the ligands were included to illustrate the spatial extension of the binding mode. For each BM GR a selected pose is depicted and the code of the ligand is indicated

Conclusions
Summarizing, in this study novel potential selective inhibitors of TR enzyme from T. cruzi, the pathogen of Chagas disease, were identified in silico. With that purpose, a library containing 2094 hybrid molecules based on the THTT heterocyclic core was designed. Such thiadiazines are molecules with a broad spectrum of biological activities that has been neglected in the last years due to their general cytotoxicity; however, we believe that is plausible to modulate their biological action by chemical modification. Having this hypothesis in mind, we considered ligands modified with peptoid or tetrazole moieties, which can be easily introduced with multicomponent reactions. The molecules were docked towards TR using the freely available software AutoDock Vina. A further molecular docking on human enzyme GR was carried out to explore the possible cytoxicity of the most promising candidates. The parameters chosen to analyze the results  Hydrogen bonds (green), cation-(magenta), -stacked (violet) and salt bridge (yellow) contacts are also highlighted. Residues with van der Waals interactions are displayed as cyan sticks. All distances are given in Å and select the most promising ligands were the population of conformational clusters and contact-based analysis of the enzyme-ligand complexes. Based on the Rosenfeld-Olson criterion, the novel numerical index called Poses ratio (P r ) was proposed to analyze the clustering results in a straightforward manner. When P r was calculated for the different substituents within the structures of interest, fragments with a tendency to form highly populated binding modes were identified, and hence, valuable substituent-interaction relationships could be determined. We envision that this new parameter may be particularly useful as a support of other established tools in high throughput virtual screenings where thousands of compounds are considered for lead identification. In the present particular case, it was found that medium size ligands with a few rotatable bonds and large substituents are the most promising inhibitors of TR. A contact-based analysis carried out with BINANA software, permitted to sort the ligands into six binding modes, according to their interaction pattern with ten relevant subsites within the active site of TR. These subsites are described in the literature to play a crucial role in the binding of the natural substrate and known inhibitors of TR. The most recurrent amino acid contacts were found, encompassing hydrogen bonds, polar and other non-covalent interactions. To the best of our knowledge, no other docking study filtering the ligands according to their interaction profile has been previously performed in the quest of active compounds against T. cruzi. At the end of the study, 103 ligands that have a plausible selective profile in their inhibitory action were recognized. The filtering docking-based protocol presented here may be extended to other enzymes of interest.

Author contributions
The project was conceived and headed by JC-B, whereas the computational work was supervised by EM-C. JEA-R performed the calculations, optimizations and the preparation of the supplementary material. All the authors collected, organized and analyzed the data. ERL-R and GMO-C prepared the first draft of the manuscript, while all the authors revised and contributed to following versions. DA made the artwork. All the authors read and approved the final manuscript prior to submission.

Availability of data and material
The online version of this article contains supplementary material, which is available to authorized users.

Conflicts of interest The authors declare that they have no conflicts of interest
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.