Graves' disease (GD) is an autoimmune Thyroid disorder (AITD), typically characterized by the presence of thyrotoxicosis, goiter, and ophthalmopathy (Grave’s ophthalmopathy) [1]. It is one of the major causes of hyperthyroidism in the geographical areas with Iodine abundance [2]. GD has reported a frequency of 20-30 cases per 100,000 individuals each year, with the possibility 3% of women and 0.5% of men acquiring GD during their lifetime [3]. The main ground cause of GD is the binding of circulating IgG antibodies to G-protein– coupled thyrotropin receptor (TSHR), resulting in the activation of Thyroglobulin gene which codes for thyroid hormonestriiodothyronine (T3) and thyroxine (T4) [4,5].

These antibodies are thyroid-stimulating antibodies (TSAbs) [6], - often called TSH receptor autoantibody (TSHRAb or TRAb) [7] that replicates the effect of the TSHR’s natural ligand Thyrotropin (TSH) by stimulating cyclic adenosine monophosphate (cAMP)–dependent signal transduction, thereby activating the Thyroglobulin gene present within the thyroid gland [6,8]. Thyroglobulin gene is naturally activated by hypothalamus upon the binding of TSH to TSHR. Since the activation of Thyroglobulin gene by the binding of TSIs to THSR is not natural, there is no feedback system leading to the excess production of thyroid hormones: triiodothyronine (T3) and thyroxine (T4) that is catalyzed by Thyroid Peroxidase (TPO) [6,9]. This activation also triggers follicular hypertrophy and hyperplasia, thyroid enlargement [4].

The synthesis of triiodothyronine T3 and thyroxine T4 hormones from Thyroglobulin gene is catalyzed by Thyroid Peroxidase (TPO) upon the iodination and coupling of tyrosyl residue in Thyroglobulin [10] and is hence the target site of anti-thyroid drugs for the treatment of GD. TPO_201-500 is a 933 residues transmembrane homodimer [11-13], affixed to the apical membrane of thyroid follicle cells by its TMD region on its C-terminal end (residues 847–871) [14,15]. The Xray crystallised structure of TPO was determined 16 years ago, but of low resolution [16,17]. TPO shares 47%, homology with eosinophil peroxidase (EPO), 48% with Lactose Peroxidase (LPO), and 47% with myeloperoxidase (MPO) [18,19]. Based on these homology studies, it was found that Asp238, His239, Glu399, and His494 were important for proximal heme linkage. The proximal histidine residues present in heme peroxidases are crucial in the redox properties of heme iron for catalysis [20,21]. The homology modeling of TPO was done recently, reporting that its cis and trans conformations are involved equally in its antigenicity [22].

The antithyroid drug (ATD), Propyl- thiouracil (6-propyl-2- thiouracil) (PTU) and Methimazole (1-methyl-2-mercaptoimidazole, Tapazole) (MMZ), are prescribed exclusively for controlling the hyperthyroidism in GD [23-26]. Both ATDs obstructs the iodination of tyrosyl residues in Thyroglobulin catalyzed by Thyroid Peroxidase, thereby arresting the production of thyroid hormones [27-29]. MMZ is preferred to PTU, due to its better adherence and bioavailability [23,30-32]; however, PTU is also recommended for pregnant women [33,34]. Irradiation studies too have been done to treat hyperthyroidism [35].

The drug binding mechanism of PTU/MMZ to TPO have only been done in-vivo [36-39], however the molecular insights are yet to be understood. Due to the close homology, LPO has been used as surrogates to study interaction of TPO with anti-Thyroid drug MMZ. MMZ-LPO catalyzed oxidation of 2,2- azio-bis-3-ethylbenthiazoline- 6-sulfonic acid (ABTS) and H₂O₂ was carried out where the turnover was estimated to be ~7.0 lM. [40]. The enzyme inhibition study between LPO- MMZ showed that inhibition occurred at the heme group and side chains Arg255, Glu258, and Leu262 of LPO. Gln105 and His109 of the LPO too were involved in the enzyme inhibition [41].

In this study, we have taken up PTU and MMZ as our drug of choice considering the disease. The protein structure of TPO_201-500 was constructed first, and its stability was asserted. The complex structures of TPO_201-500 with both the drugs were constructed using PatchDock server [42,43] which gave the results based on the maximum surface area and minimum atomic contact energy. LigPlot+ software v.1.4.5 [44] was employed to investigate the various bonding and non-bonding interactions between the drugs and TPO_201-500. Since the residues- Asp238, His239 and His494 belonging to the heme linkage region of the protein TPO were seen to be involved in the bond interaction with both drugs, these sites will now no longer be available for the iodination and coupling of tyrosyl residues by TPO for the synthesis of thyroid hormones. R491 was also found to play crucial role in the drug-protein communication by forming hydrogen bond. Binding energy calculation further confirmed that both drugs possessed high affinity towards TPO_201-500.

The antechamber [45] module of AMBER12 [46] software package was used for prepping the ATDs prior to the simulation. Each of the TPO_201-500-ATD complexes was subjected to Molecular Dynamics simulation to study its stabilizing properties. In addition, Poisson– Boltzmann and generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA) [47-51] techniques were endorsed to calculate binding free energies of MMZ and PTU to TPO_201-500. PBTOT for TPO_201-500-MMZ and TPO_201-500-PTU was -11.05 and -14.02 kcal mol-1 while GBTOT for TPO_201-500- MMZ and TPO_201-500-PTU was found -9.43 and -15.47 kcal mol-1 respectively, which put forth the inhibitors to have high affinity to TPO_201-500 in the dynamic system. Our findings affirm that both the drugs are equally effective as an anti-thyroid drug. Moreover, our results may further provide guidance to developing more potent inhibitors of TPO and preventing thyroid hormone production in GD.

Computational model

The amino acid sequences were obtained from UniProt [52], ID No. P07202 by querying the TPO_201-500 protein sequence against RCSB PDB web server [53,54] domain of Thyroid Peroxidase protein sequence (residues 201-500) was submitted to Iterative Threading ASSembly Refinement server (ITASSER) [55-57] for the 3D structure prediction. From I-TASSER we obtained 5 best models ranked according to their C-scores out of which best model is chosen based on C-score, Template Modelling (TM) score [58] and Root Mean Square Deviation (RMSD) score. Using RAMPAGE [59] and VERIFY 3D server [60,61] we validated the structure of TPO_201-500.

Based upon the 3D structural validation checks, Model 1 showed consistently good results among the rest, hence chosen as an ideal structure of TPO_201-500 for further studies. SOPMA (Self-optimized prediction method alignment) server [62] was adopted in studying the secondary structure of TPO_201-500. Anti-thyroid drug: Propylthiouracil (PTU) and Methimazole (MMZ/MMI) were chosen for this study. 3D conformer of both the ATDs were procured from PubChem server [63] in SDF format, which was later converted to pdb file using open Babel software v.2.3.2 [64].

Docking

The initially constructed 3D structures of TPO_201-500 with PTU and MMZ were processed separately for the rigid molecular docking using the PatchDock [42,43] web server. PatchDock works on the principle of geometric surface docking algorithm. The clustering score of 1.5 Å RMSD was given to remove the unwanted models. The ten complexes of TPO_201-500-PTU and TPO_201-500-MMZ were obtained from the server. Based on highest surface area and lowest atomic contact energy, one complex each of TPO_201-500-MMZ and TPO_201-500-PTU were selected.

Determining interacting residues and bond types

The best structure of TPO_201-500 in complex with PTU and MMZ processed by PatchDock was selected among the rest considering the minimum atomic contact energy and maximum surface area. The best structure was visualized in LigPlot+ software v.1.4.5 [44] for evaluating the interacting residues of TPO_201-500 interacting with the drugs. Interface residues are defined as the residues with a contact distance less than 6 Å from the interacting partner.

Molecular dynamics simulations of thyroid peroxidase (TPO_201-500) with anti-thyroid drugs (ATD) - methimazole and propylthiouracil

The molecular dynamics study was performed using AMBER12 [46] software package where ff99SB force field [65] parameters were used for protein (TPO_201-500) and generalized AMBER force field (GAFF) [66] parameters of the anti-thyroid drugs−Methimazole and Propylthiouracil were established with the Antechamber [45] program, wherein AM1-BCC [67,68] was used to fix the proportionate partial atomic charges. Using xLEaP and antechamber module in AMBER12, we prepared the initial coordinate and the topology files for two the complex systems.

Xleap also added all the missing parameters to the complexes. All the canonical systems were then neutralized with Na+ ions, followed by the hydration with TIP3P [69] water molecules in a box with a buer size of 10 from the solute in x, y, and z axis. The initial minimization of the solvated structures was done imposing the constraints over the solute for the first 500 steps using the steepest descent algorithm followed by another 500 steps with a conjugate gradient method. Second minimization was carried out freely for 500 steps. 8 Å was applied as the cut off for both the minimization steps. Gradual heating was done to both the complexes with the successive increase in the temperature from 0-300K. The generated ensemble was then equilibrated for 100ps using the NPT ensemble at 300 K and 1 atmospheric pressure. Full MD production run was carried out till 15ns.

Periodic boundary conditions along with all the electrostatic interactions were maintained throughout the simulations using the particle mesh Ewald (PME) method [70,71] under isothermal and isobaric conditions respectively. The temperature was governed using Berendsen thermostat [72]. The production run was attained with the NPT ensemble for all the formerly mentioned Thyroid Peroxidase-Anti thyroid drug complexes for 15 ns with 2 fs time step. Furthermore, the shake algorithm [73] was used to restrain all the bonds including the hydrogen atoms.

MD analyses

To check the stability of the obtained MD trajectories, each trajectory were clustered to acquire comparative Root Mean Square Deviation (RMSD), Radius of Gyration (Rg), B-factor and Solvent Accessible Surface Area (SASA) using the CPPTRAJ algorithm [74] of Amber12 package. We also proceeded further to analyze hydrogen bond properties based on all possible hydrogen donors (HD) and acceptors (HA) in the system, the percentage occupancy of hydrogen during each step of MD trajectory and considering bond distances and bond angle made by HA–H–HD atoms. Here, the hydrogen bonds having consistency for more than 5% of the simulation period were taken as an output result. Molecular graphics and analyses of the trajectories obtained were performed using UCSF Chimera package [75] and VMD v.1.9.3 [76].

Binding Free Energy (BFE) calculations

The binding free energy between the protein and the drugs were assessed by considering the thermodynamics pattern of Figure 1. The binding free energies of two anti-thyroid drugs- Methimazole and Propylthiouracil to Thyroid peroxidase were computed using MMGBSA (Molecular Mechanics Generalized Born Surface Area) and MM-PBSA (Molecular Mechanics Poisson–Boltzmann Surface Area approach [47-51] in AMBER program. We took 100 conformational snapshots extracted from the last 1 ns of MD trajectories. Independent MM-GBSA/PBSA analysis was done for (i) Thyroid peroxidase (TPO_201-500) (ii) anti-thyroid drugs (ATD) - MMZ/PTU and (iii) the complex- TPO_201-500-MMZ/PTU. The binding free energies of TPO_201-500 with both the drugs were calculated using equation 1.

ΔG_binding = ΔG_complex- [ΔG_TPO+ ΔG_ATD] --------------->(1)

Figure 1: Schematic representation of the thermodynamic cycle used for the calculation of binding free energies (BFE) of Thyroid peroxidase (TPO) and Anti-Thyroid drugs (ATD).

Where G-binding corresponds to the total binding free energy, Gcomplex, GTPO and GATD are the solvation free energies of the complex, thyroid peroxidase and anti-thyroid drugs, respectively. BFE as stated by second law of thermodynamics can also be summed as the difference between enthalpy ΔH and the conformational entropy -TΔS [48,51] of the system which can be represented by the equation 2.

ΔG_binding = ΔH - TΔS ---------------> (2)

ΔH can be explained as sum of the changes of the gas-phase molecular mechanics energies (ΔE_MM), polar and non-polar solvation energy (ΔG_solv); see equation 3. Also, ΔE_MM is the sum total, of the differences in the bond, angle and dihedral energies (ΔE_internal), van der Waals (ΔE_vdW), and electrostatic energies (ΔE_ele), as shown in the equation 4. Accordingly, equation 5 divides ΔGsolv into polar (ΔG_PB/GB) and non-polar solvation free energy (ΔG_surf);

ΔH = ΔE_MM + ΔG_solv ---------------> (3)

ΔE_MM = ΔE_internal + ΔE_vdW + ΔE_ele ---------------> (4)

ΔG_solv = ΔG_PB/GB + ΔG_surf ---------------> (5)

The terms E_int, E_vdW, and E_ele are internal energy, van der Waals interaction, and electrostatic interaction of the systems respectively. The terms ΔG_PB/GB and ΔG_surf are polar and non-polar solvation free energies. The modified GB model of Onufriev [77] was considered in the determination of the ΔG_GB. A value of 80 was considered for the internal dielectric constant of the solvent (water) while, the external dielectric constant of the solute was set to 1. The non-polar solvation free energy (ΔG_surf) for both the systems was estimated equation 6.

ΔG_surf = γ × SASA +β ---------------> (6)

Where SASA is the abbreviation of solvent-accessible surface-area. In case of PB method, γ was to 0.00542 kcal (mol^-1 Å^-2) and β = 0.92 kcal mol^-1, while for the GB method, γ was customized to 0.0072 kcal (mol^-1Å^-2) and β = 0 kcal mol^-1 respectively [78]. 1.4 Å was the probe radius of the solvent used for these calculations.

Structural modelling and validation of thyroid peroxidase

We modelled the structure of Thyroid Peroxidase (TPO_201-500) from I-TASSER (Figure 2). The server provided us with five models using the threading algorithm. The ideal model was selected based upon the C-score computed from the comparative structural significance. The C-score for TPO_201-500 was 1.20. The overall quality of model was achieved by various bioinformatics tools. TM scored TPO_201-500 in the rage of 0.88 ± 0.07 and RMSD estimated was 3.8 ± 2.6 Å. The structural justification of the ITASSER yielded structure was done by RAMPAGE server. We observed 81.5% of residues were within most favored region, 11.4% in allowed region and 7.0% in disallowed region. Verify-3D server ascertains the accuracy of the predicted secondary structure (3D) model with its respective residues (1D) by assigning a structural class based on its location and environment. 81.6% of residues had an averaged 3D-1D score >= 0.2 which confirms the accuracy of the predicted model to near correct to the actual structure. SOPMA server identified 39.33% α-helices, 47.67% random coil, 7.33% or extended strand and 5.67% of β-strands respectively in the modeled structure of Thyroid Peroxidase.

Figure 2: Modelled structure of Thyroid Peroxidase (TPO) obtained from ITASSER server with the C-score of 1.20.

Protein-drug interaction study

It is an evident fact that Methimazole and Propylthiouracil interferes in the synthesis of thyroid hormones by regulating the enzyme inhibition of thyroid peroxidase. The in-vivo/in-vitro binding mechanisms of PTU/MMZ to TPO_201-500 have also been done exclusively, yet the interactions at the molecular level are not known. So as to understand the drug-protein interaction at the molecular stage, we docked Thyroid peroxidase with both the drugs using PatchDock, server. This particular server works on the principle of geometric hashing which is implied as a best molecular shape complementarity between the docked structures considering the wide interface area and lesser steric hindrances. The docking gave us array of docked structures for TPO_201-500-MMZ/PTU complex. These complexes were graded basis of geometric shape, interface area and atomic contact energy as shown in Figures 3 and 4. The best ranked structure in case of TPO_201-500-MMZ complex had a geometric shape complementarity score of 2300, approximate interface area of 260.2 Ų and atomic contact energy of (ACE) of -152.3 kcal mol^-1 (Solution Structure 1 from Figure 3). Likewise, the top scored structure for TPO_201-500-PTU was found to have the geometric score of 2912, interface area of about of 328.8 Ų and ACE of -146.6kcal mol^-1 (Solution Structure 1 from Figure 4).

Figure 3: Docked structures of thyroid peroxidase-methimazole (MMZ) obtained from patchdock server.

Figure 4: Docked structures of thyroid peroxidase and propylthiouracil (PTU) obtained from patchdock server.

Next, we examined the hydrogen and hydrophobic interaction occurring in protein-drug complex using the Ligplot+ tool as shown in Figures 5A and 5B. Figure 5A details the interaction between Thyroid Peroxidase and Methimazole. Arg491 formed a hydrogen bond with the thiol atom i.e., Sulphur at the distance of 3.25 Å. Hydrophobic bonds are formed between the protein and the drug by residues Asp238, His239, Phe243, Thr487 and His494. Figure 5B showcases the molecular interaction between TPO_201-500 and Propylthiouracil. We can clearly observe the hydrogen bond formed between residue Arg491 and Sulphur atom of the drug PTU at the distance of 3.04 Å. Amino acids Asp238, Gln235, His239, Phe243, Glu399, Thr487 and His494 formed the hydrophobic bonds with the drug. In both the scenarios, drug is seen to have bonded with proximal heme binding residues of the protein- Asp238, His239 and His494.

Figure 5: LIGPLOT representation of Anti-Thyroid drugs-Methimazole (MMZ) and Propylthiouracil (PTU) binding to Thyroid Peroxidase (TPO). Ligand, protein, and hydrogen bonds are in thick purple, thick brown and green dash lines, respectively. Hydrophobic contacts forming residues are shown in sky blue colour as an arc with spikes towards ligand atoms (A) Thyroid Peroxidase (TPO) bound to Methimazole (MMZ). (B) Thyroid Peroxidase (TPO) bound to Propylthiouracil (PTU).

Since these residues are bound to the drug, which renders the active site of the protein unavailable for its ritual catalyzing activity thereby inhibiting the production of the thyroid hormones. Another homolog study of TPO was done in LPO by Singh et al., where they found that the heme moiety residues held MMZ tightly on one side while hydrophobic amino acids - Arg255, Glu258, and Leu262 held it on the other end. Many extensive studies on LPO as TPO surrogate has supported the importance of heme binding and catalytic site residues mainly Hisidines and Aspartic acid in mammalian peroxidases [79-81]. A review on MPO highlighted its ligand binding site to be present within the distal heme pocket which upon the binding of PTU, isoniazid could lead to irreversible inactivation of the enzyme [82]. The homology studies done by Taurog in 1999 and Furtmüller et al. have pointed out the essentiality of the proximal histidine residues for the prime functioning of TPO_201-500 to catalyse the iodination of tyroysl residues in thyroglobulin. And we have noted these histidine residues His239 and His494 to have been engaged by anti-thyroid drugs due to hydrophobic linkage, also relieving the protein of its catalytic role. Comparative study on the effective dosage of Methimazole and Propylthiouracil in Patients with GD by Nakamura et al. [83] revealed that 15 mg/d of MMI is prescribed for mild and moderate GD, 30 mg/d of PTU is recommended for severe forms of GD. A team of scientists from Spain and Argentina [84] determined the pharmacological activities of a copper based propylthiouracil ([Cu(PTU)₂]₂, and concluded that this combination of PTU would be good candidate due to its antioxidant, antimicrobial and alkaline phosphatase activity.

MD simulation study on the TPO_201-500-MMZ/PTU complex

The molecular dynamics simulations were administered to probe the efficacy of both the anti-thyroid drugs- Methimazole and Proplythiouracil as an inhibitor for thyroid peroxidase. Extracting the trajectory output, we analyzed the stability of TPO_201-500-MMZ/PTU complexes in an explicit environment. The conformational dynamics of TPO_201-500 complexed with the drugs-MMZ and PTU were estimated by plotting the Root Mean Square Deviation (RMSD) by taking their respective equilibrated structures as a reference. The comparative RMSD graph for the Cα atoms of the TPO_201-500- MMZ/PTU complexes as a function of time can been seen in Figure 6A. MMZ bound TPO_201-500 started to settle at the RMSD value of around 5 Å while MMZ bound XPA was exhibiting fluctuations even at the RMSD value of ~ 7.5 Å. The stability of PTU over MMZ in the biological system can also be explained from the study done by Yoshihara et al., where they found that over exposure of the in utero foetus to MMZ increased the risks of congenital malformations and least with PTU prescriptions.

Figure 6: Comparative MD analyses of TPO-MMZ and TPO-PTU. (A) RMSD values for TPO-MMZ and TPO-PTU complexes relative to the starting structure during MD simulation. (B) Solvent accessible surface area (SASA) values for TPO-MMZ and TPOPTU complexes. (C) Radius of Gyration (Rg) values of TPO-MMZ and TPO-PTU. (D) B-factor distribution for TPO-MMZ and TPOPTU.

Figure 6B shows the solvent accessible surface area (SASA) for the both the protein-drug complexes. Both the systems, be it TPO_201-500- MMZ or TPO_201-500-PTU had proper access to the solvent provided in the water box throughout the simulation. TPO_201-500- MMZ presented the SASA value of within 20000-22000 Ų while PTUshowed SASA value in the range of 21000-22500 Ų. TPO_201-500- PTU complex required much larger surface area to access the solvent in the system.

We determined the comparative radius of gyration (Rg) for both the drugs complexed with TPO_201-500. Rg quantifies the compactness and the spatial dispersal of the atoms, of a particular component molecule, around an axis. Figure 6C shows the radius of gyration for TPO_201-500 bound to Methimzole and Propylthiouracil as a function of time. Methimazole bound TPO_201-500 showed oscillations from 6 ns onwards, while TPO_201-500-PTU was settled till 9 ns, then it started oscillating.

Next, we analyzed B-factor graph for the backbone Cα-atoms for the TPO_201-500 in complex with both the drugs to observe the variations in their structural dynamics from their average position. Figure 6D shows that both the systems presented more or less coherent atomic fluctuations, and the only difference between them being the higher peaks as seen in the TPO_201-500-PTU complex, whereas TPO_201-500-MMZ had shorter peaks meaning less atomic fluctuations.

Hydrogen bond analysis as seen in Figure 7 showed intermolecular hydrogen bonds formed between the enzyme and the drug during the simulation period. As per Figure 7A, we saw two hydrogen bonds formation in TPO_201-500-PTU complex form, while only single hydrogen bond formed in case of TPO_201-500-MMZ (Figure 7B) during simulation. We further calculated the inter-molecular hydrogen bond stability of both the systems.

Figure 7: The number of hydrogen bonds within the two systems during MD simulation. (A) two H-bonds formed between TPOPTU. (B) Single H-bonds formed between TPO-MMZ.

Table 1 shows hydrogen bond occupancy of the inter-molecular hydrogen bonds formation as the function of time period. We found the important amino acids Arg491 had high fraction of occupancy in the complex of TPO_201-500-MMZ while Asp238 and Phe485 had greater fractions for hydrogen bond formation in case of TPO_201-500- PTU. Higher fraction of hydrogen bond definitely suggested the stronger interaction between the protein TPO_201-500 and the ATDs. These occupancies align exactly with the Figures 7A and 7B for the total number of the hydrogen bonds formed within the systems. We simultaneously took comparative snapshots of both the TPO_201-500- MMZ and TPO_201-500-PTU systems at different intervals of simulation period and are presented in Figure 8.

Figure 8: Comparative snaphots for TPO-MMZ and TPO-PTU at 3 ns, 6 ns, 9 ns and 15 ns respectively.

Table 1: H-bond occupancy among all MD simulations.

Both the systems presented the β-sheets formations. PTU bound TPO_201-500 had loss of the α-helices over the simulation period.

Calculations of binding free energies

To measure the efficacy of the anti-thyroid drugs- Methimazole and Propylthiouracil in inhibiting the enzyme activity of Thyroid peroxidase in Graves’ disease, the binding free energies were evaluated using MM-PBSA/GBSA method as provided in AMBER program. The calculated values are summarized in Table 2. According to Table 2, the binding free energy of the TPO_201-500-MMZ complex using MM-PB/ GBSA methods is found to be -11.05 (PBTOT) and -9.43 (GBTOT) kcal mol^-1. TPO_201-500-PTU complex had the binding free energy of -14.02 (PBTOT) and -15.47 (GBTOT) kcal mol^-1. As per the data obtained, the drug-protein binding was considerably enhanced by the contributions of the van der Waals interactions (ΔEvdW), non-polar solvation energy (ΔGsurf) and electrostatic interactions (ΔGele).

Table 2: Binding free energies (kcal mol^-1) between Thyroid peroxidase (TPO) and Anti-Thyroid drugs- Methimazole (MMZ) and Propylthiouracil (PTU).

The PBTOT or GBTOT values of both the complexes showed not much difference in their overall binding energy. For PBTOT, the differences were 14.02-11.05 = 2.97 kcal mol^-1, while their GBTOT difference was however quite large i.e., 15.47-9.43 = 6.04 kcal mol^-1. Based on these results, Propylthiouracil showed better results in enzyme inhibition than Methimazole. Nonetheless, both the drugs are equally qualified in target binding of Thyroid peroxidase and thereby preventing the production of thyroid hormones. A population study conducted by Azizi and his team in 2014 for the management of hyperthyroidism in Asia too showed the preference of PTU over MMZ by the clinician’s due to its better outcome. A density functional theory method was opted by Tappa and his group [85] to evaluate the inhibitory effect of the MMZ tautomers-1-methyl-1H-imidazole-2 (3H)-thione (M1) and 1-methyl-1H-imidazole-2-thiol (M2) where the effect of M1 was superior to M2 tautomer. Gupta and Kumar [86] did metabolic modelliing and tested the drug activity of PTU against TPO, TSHR and sodium iodide symporter (NIS) using systems biology approach and deduced that PTU is effective on TPO. They also observed TSHR to be potent target for therapeutic study. A study done on PTU, KClO4, or TSH [87] revealed that these compounds reduced the oxidative stress due to high concentrations of iodide in the thyroids of both Metallothioneins(MTs) MT-I/II KO and WT mice. Bhabak and Mugesh in 2010 [88], and Manna, Roy and Mugesh in 2013 [89], reported Selenium analogues of ATDs to have better inhibitory effect compared to Sulphur based thinamide drugs on TPO suggesting that this selenium based drugs may be the next potential ATDs.

In this computational study, the protein structure of thyroid peroxidase was modelled first, followed by the independent docking of the protein with MMZ and PTU. We have identified the residues that are involved in the drug-protein interaction. Both the drugs engaged the heme linkage active site residues- Asp238, His239 and His494, depleting the protein of its function to catalyze the iodination of tyrosyl residues on thyroglobulin for T3 and T4 hormone production. The hydrophobic interaction with residues Asp238, His239, Phe243, Thr487 and His494 was seen in both the cases. Proximal histidine residues His239 and His494 were also found to be interacting with the drug making the enzyme unfit to carry out its biological role. Arg491 was also found to be binding with the drug by hydrogen bonding. RMSD showed the stability of MMZ to be much higher than PTU. TPO_201-500-PTU had higher values for SASA and B-factor graphs. Hydrogen bond analyses too substantiated the stronger and stable interaction shared by these anti-thyroid drugs to the enzyme. MMPBSA/ GBSA binding free energy results showed both drugs to be equipped in binding to TPO_201-500 and in inhibiting the catalytic activity of TPO_201-500 in thyroid hormone synthesis as seen in the literatures. Propylthiouracil showed comparatively higher binding efficiency than methimazole as per GBTOT analysis, yet there was not much difference between their binding energies considering the PBTOT results. After studying all the obtained evidences, we deduce that both the anti-thyroid drugs are competent enough to inhibit the enzyme activity of Thyroid peroxidase, and in process keeping the unwanted production of thyroid hormone at bay.

The authors extend their deepest gratitude towards Tezpur University for the start-up grant. This research was also supported by Department of Biotechnology (DBT) funded Bioinformatics Infrastructure facility (BIF) in the Department of Molecular Biology and Biotechnology at Tezpur University.