Thioredoxin sequence forms the basis of our study as this is a good target for cancer chemotherapy. We selected the gene sequences from UniProt KB, which is commonly used as knowledge base for molecular sequences. Most of the sequences in UniProt KB are derived from the conceptual translation of nucleotide sequences. It plays an important role by providing a stable, comprehensive, freely accessible central resource on protein sequences and functional annotation. We used computational analysis for the functional annotation for the gene sequences.

Position Specific Iterative BLAST (PSI-BLAST) profile was generated from local alignments of the most highly scoring hits in the initial BLAST results by calculating position-specific scores for every position in the alignment in the sequences. Five template sequences were generated based on this alignment program. This iterative procedure increased the sensitivity of the BLAST search and helped us to identify new relationships between the query and database entries. Clustal W was used for Multiple Sequence Alignment Program (MSA) to determine the conserved sequences among the templates.

As the PSI- BLAST similarity obtained was less than 30%, we preferred the fold recognition method as the option to build 3-D model. Usually fold recognition methods are so efficient especially in the following cases: first, when the sequence has little or no primary sequence similarity to any sequence with a known structure. Second, when some model from the structure library represents the true fold of the sequence.

Our study falls into the first category in which we tried to recognize the structural fold of the target protein from a structure template library, given its sequence information and then generate an alignment between the query and the recognized template protein, from which the structure of query protein was predicted. We used I-TASSER web server (8) which has generated five predicted 3D models for the requests. A scoring function (C-score) based on the relative clustering structural density and the consensus significance score of multiple threading templates were obtained to estimate the accuracy of the I-TASSER predictions.

The output of the I-TASSER server for our query protein included the prediction of secondary structure, top five full-length models with confidence scores, the estimated TM-score, RMSD, standard deviation of the estimations and top ten templates. The binding site predicted by I-TASSER server suggests 26 amino acids residues as the possible binding site residues.

Modeller 8 was used to construct the model based on the generated templates for the following PDB ids 1G4MA, 2WTRB, 2R51A, 1CF1C, 2FAUA. All templates were taken based on their folds to construct the model.

Functional characterization of a protein sequence is a common goal in biology, and is usually facilitated by having an accurate three-dimensional (3-D) structure of the studied protein. In the absence of an experimentally determined structure, comparative or homology modeling can sometimes provide a useful 3-D model for a protein that is related to at least one known protein structure. Comparative modeling predicts the 3-D structure of a given protein sequence (target) based primarily on its alignment to one or more proteins of known structure (templates). The prediction process consists of fold assignment, target-template alignment, model building, and model evaluation (9).

The following threading programs were used to collect the templates:

MUSTER

This model was optimized using GROMACS (10) which is a molecular dynamics package primarily designed for biomolecular systems such as proteins and lipids. The minimized model obtained was used for virtual screening in order to filter the compounds using GOLD score.

The use of virtual screening to discover new inhibitors is becoming a common practice in modern drug discovery (11).

Receptor-based virtual screens seek to “dock” members of a chemical library against a given protein structure, predicting the conformation and binding affinity of the small molecules (12). All possible compounds which obey rule of 5 were collected from PUBCHEM, which is the main resource for obtaining freely-available bioassay data provided by the National Center for Biotechnology Information [NCBI] (13).

ZINC is a free database for virtual screening that contains over 4.6 million compounds in ready-to-dock, 3D formats, available at the URL http://zinc.docking.org. Molecules in ZINC are annotated by molecular property that include molecular weight, number of rotatable bonds, calculated LogP, number of hydrogen-bond donors, hydrogen-bond acceptors, chiral centers, chiral double bonds (E/Z isomerism), polar and a polar desolvation energy (in kcal/mol), net charge and rigid fragments. The database contains 494,915 Lipinski compliant molecules and 202,134 ‘lead-like’ molecules, having molecular weight in the range 150 to 350 with calculated LogP <4, number of hydrogen-bond donors ≤3 three, and number of hydrogen-bond acceptors ≤6. ZINC provides several search criteria such as molecular property constraint, ZINC codes, vendor based, and molecular substructure search (14).

The sorted ligands again were docked using GOLD and their interactions were studied. The protein was loaded into GOLD wizard, tailed by addition of hydrogens, then selection of binding site residues, choosing score (GOLD), and then the program. It is the first algorithm to be evaluated on a large dataset of complexes (15). The result gives an empirical free energy scoring function that estimates the free energy of binding of molecules.

The Chemscore function was implemented as a scoring function for the protein–ligand docking program GOLD, and its performance compared to the original Goldscore function in terms of docking accuracy, prediction of binding affinities, and speed. In terms of producing binding energy estimates, the Gold-score function appears to perform better than the Chemscore function and the two consensus protocols, particularly for faster search settings. Even at docking speeds of around 1–2 min/compound, the Goldscore function predicts binding energies with a standard deviation of ∼10.5 kJ/mol (16). Based on the interaction studies and the Lipinski's rule of five, four molecules were selected as the hit molecules for further analysis in the wetlab. ChemDraw is the software used to sketch and build the molecular models.

Utilizing PSI-BLAST in the National Center for Biotechnology Information (NCBI), we searched the sequences based on the query (UniprotKB ID: Q9H3M7) with which we generated the five templates for Thioredoxin which are shown in (Figure 1). PSI-BLAST program increased the sensitivity of the BLAST search and helped us to generate the possible alignment between the thioredoxin query sequences and database sequences. It also constructed a position-specific substitution matrix for the query sequence during the search. The obvious advantage of this approach was that we could obtain more accurate statistical estimates, along with the application of BLAST/PSI-BLAST in terms of speed of aligning the sequences. This gave us five corresponding PDB hits and alignments. Furthermore, by using the Multiple Sequence Alignment Program we determined the conserved positions which are shown in figures 1A and 1B.

We constructed 3-D models of these templates utilizing the I-TASSER web server. The target sequences were first threaded through a representative PDB structure library to select the best possible models from PDB, based on the templates. Using I-TASSER we obtained the scoring values as C-score, TM-score (0.08) as well as RMSD (2 Å) values. Based on these models we could find the binding site residues. By aligning the templates and refining the iterative threading assembly (I-TASSER), Web server was used to integrate the platform for automated protein structure and function prediction. And we could generate five complete models with Combining confidence score (C-score) and protein length.

The advantage of using I-TASSER is that it first generates three-dimensional (3D) atomic models from the amino acid sequences, and then it does multiple threading alignments and iterative structural assembly simulations. Hence, the function of the protein could be inferred by structurally matching the 3D models with other known proteins. The output from the webserver run was important because it contained the information of full-length secondary and tertiary structure predictions, and functional annotations on ligand-binding sites, enzyme commission numbers and Gene Onto-logy terms. The accuracy of the predictions was based on the C-score of the modeling which depicted the best equivalent residues of two proteins based on the structural similarity and the output of TM-score.

Our search for the binding sites in the modelled structure was done by using I-TASSER server with which we could locate the exact positions of various amino acid residues at their respective binding sites. We made sure that these residues were in the binding pocket within the vicinity of the active site of the modelled protein, as shown in the results below:

PHE:14 ASN:15 ASP:16 PRO:17 GLU:18 VAL:20 VAL:146 ASP:147 VAL:149 PRO: 331 PRO:332 CYS:333 TYR:334 HIS:342 ARG:343 LEU:344 GLU:345 SER:346 TYR:366 GLU:369 PHE:370 MET:373 PRO:374 PRO:376 TYR:378 THR:379

The identification of putative ligand-binding sites on proteins is important for the prediction of protein function. Knowledge-based approaches using structure databases have become interesting, because of the recent increase in structural information (Figure 2).

Verification of the built model was done to ensure whether the model was programmed correctly and the algorithms were implemented. Validation results determined that the distribution of amino acid residues were at the most favourable region in the Ramachandran plot (more than 90%). This is an indication of the stereochemical quality of the model taken for the structural analysis, it and also validated the target-ligand binding efficacy of the structure. Ramachandran plot displays the main chain torsion angles phi, psi (φ, ψ); (Ramachandran angles) in a protein of known structure (Figure 3).

Dihedral angle checks Ramachandran plot shows phi-psi distribution. Each residue is classified according to its region: ‘core’, ‘allowed’, ‘generous’, or ‘disallowed’. Residues in the generous and disallowed regions are high-lighted on the plot. A log-odds score shows how normal or unusual the residue's location is on the Ramachandran plot for the given residue type. Procheck results gave us the value of 96.3% residues in most favoured regions in R-Plot which suggests that they predict thioredoxin model of good quality (Figure 4A).

ERRAT is another program which we used for verifying crystal structures. Error values in this program are plotted as a function of the position of a sliding residue in the window. The error function is based on the statistics of non bonded atom-atom interactions in the modeled structure. ERRAT prompts the models to have the overall quality factor to be above 95%. And our results have shown that the value of overall quality factor was 95.506%. This confirmed that our developed model had reliable high resolution and quality compared to a database structures considered for the study (Figure 4B).

Further we analyzed the compatibility of our predicted (3D) model by utilizing Verify-3D program with its own amino acid sequence (1D). For each residue of the amino acid the scores of a sliding 21-residue window (from -10 to +10) were added and plotted. The returned 3D-1D profile showed 3D-1D score of above 0.2 and consisted of 92.33% of the residues in the predicted model (Figure 4C).

Statistical analysis in Ramachandran plot compared well with the observed and expected distributions of experimental observables and provided powerful tools for the quality control of our protein structure. The distribution of backbone dihedral angles (‘Ramachandran plot’) have often been used for such quality control, but without a firm statistical foundation. Hence the output for a protein structure is a Ramachandran Z-score, expressing the quality of the Ramachandran plot relative to current state-of-the-art structures.

Loop optimization for generating our model was done by using the software Modeller 9V8 (17). The model had initial potential energy=55550610.692 and initial RMS gradient=1115700.680. For energy minimization of our model we used GROMACS program. This was done by using steepest descent algorithm for 1000 steps, and 5000 steps for conjugate gradient algorithm we obtained its potential energy as-931094.12296 and RMS gradient as 0.63272. The models are presented in (Figures 5A and 5B).

We selected the ligands from PUBCHEM and ZINC data bases (almost 400), and by virtual screening using MERCURY and MARVIN VIEW we could shortlist 4 ligands as the best fit ligands. These are listed below and are presented in table 1.

3-hydroxy-2,3-diphenylbutanoic acid,

Since the results depend on the choice of scoring functions obtained in GOLD (Table 2), an analysis was performed based on the ligand binding score greater than 50. Also, there was no overlap between the top-scoring compounds from protein-ligand versus ligand-based scoring (Table 2). The small molecules which we determined involved compounds with similar chemical structures, similar modes of action, or drug interactions.