3D-QSAR, molecular docking and molecular dynamics studies of a series of RORγt inhibitors

The discovery of clinically relevant inhibitors of retinoic acid receptor-related orphan receptor-gamma-t (RORγt) for autoimmune diseases therapy has proven to be a challenging task. In the present work, to find out the structural features required for the inhibitory activity, we show for the first time a three-dimensional quantitative structure–activity relationship (3D-QSAR), molecular docking and molecular dynamics (MD) simulations for a series of novel thiazole/thiophene ketone amides with inhibitory activity at the RORγt receptor. The optimum CoMFA and CoMSIA models, derived from ligand-based superimposition I, exhibit leave-one-out cross-validated correlation coefficient (R2cv) of .859 and .805, respectively. Furthermore, the external predictive abilities of the models were evaluated by a test set, producing the predicted correlation coefficient (R2pred) of .7317 and .7097, respectively. In addition, molecular docking analysis was applied to explore the binding modes between the inhibitors and the receptor. MD simulation and MM/PBSA method were also employed to study the stability and rationality of the derived conformations, and the binding free energies in detail. The QSAR models and the results of molecular docking, MD simulation, binding free energies corroborate well with each other and further provide insights regarding the development of novel RORγt inhibitors with better activity.


Introduction
The nuclear hormone receptor (NHR) includes receptors for thyroid, steroid hormones, retinoids, vitamin D and different 'orphan' receptors, which are composed of DNA-binding domains (DBDs) and ligand-binding domains (LBDs) (Aranda & Pascual, 2001). These receptors display diverse functions in reproduction, development, and homeostasis by regulating cell growth, differentiation, and apoptosis, thus resulting in inflammatory diseases and cancer to endocrine and metabolic diseases (Kastner, Mark, & Chambon, 1995). Recent studies have shown that Th17 cells play important roles in immune response, inflammation, and cancer diseases (Ivanov et al., 2006;Serody & Hill, 2012;Tosolini et al., 2011), which is a linkage of CD4 + T helper cells, expressing several cytokines, IL-17A and IL-17F. The differentiation and function of Th17 cells are controlled by the transcription factor retinoic acid receptor-related orphan nuclear receptor-gamma-t (RORγt).
Unlike other NHRs, RORγt expression is sufficient to induce transcriptional activation of a reporter, which suggests that RORγt is either constitutively active or its activating ligands are ubiquitously present. Three ROR homologues are delineated in mammals: RORα, RORβ, and RORγ. RORβ is not expressed in immune cells; RORα and RORγt play partially redundant roles in mouse Th17-cell differentiation (Yang et al., 2008), and RORγt alone can block human Th17-cell differentiation in vitro. Animal experiment has found that RORγt is expressed at the double-positive stage of T-cell development (Eberl & Littman, 2004). Additionally, affymetrix gene chip analysis on CD4 + T cells cultured under conditions favoring Th17 and favoring Th1 differentiation, results demonstrated that RORγt mRNA was elevated and had central role in Th17 differentiation. Different from other transcription factors, RORγt consisting of a ligand-binding pocket is hence to be an excellent target for pharmacologic intervention in inflammatory and immune diseases. Thus, the RORγt inhibitors have potential utility in controlling the activity of Th17 cells and can be developed as therapeutic agents for the treatment of Th17-related diseases.
Studies have identified numerous structural diverse inhibitors of RORγt targeting the ligand-binding pocket. A small-scale small-molecule screening was performed, and digoxin was identified as a RORγt inhibitor by inhibiting murine Th17-cell differentiation without affecting other T-cell lineages; however, digoxin exhibited side effect (toxic for human cells) at 300 nM (Huh et al., 2011). Using the T0901317 scaffold as a lead compound, Griffin and Burris developed a derivative, SR1001, a new class of compound, which is specific to both RORα and RORγt (Solt et al., 2011). By screening a small chemical library, ursolic acid, a small molecule present in herbal medicine, was identified, which selectively and effectively inhibited RORγt, resulting in greatly decreased IL-17 expression (Xu et al., 2011). Solt et al. (2012) have reported a screening of SR1555 (IC 50 ≈ 1.5 μM) in a GAL4-NR chimeric co-transfection assay, which target both suppression of TH17 and stimulation of T regulatory cells. Using a modular chemistry approach, further chemical modification to the SR1001 scaffold has developed a RORγt specific inhibitor SR2211 . In addition, a small molecule library was screened at the NIH Chemical Genomics Center, a series of diphenylpropanamide compounds including compound ML209-a selective RORγt inhibitor were found (Huh et al., 2013). However, all these RORγt inhibitors are not suitable for oral dosing, which leads to a continued investigation of RORγt inhibitors.
More recently, novel thiazole/thiophene ketone amides as potent RORγt inhibitors have been discovered and synthesized, the binding activity was analyzed by incubating the RORγt-LBD protein with various concentrations of compounds in the assay buffer, the IC 50 value was determined using the Prism software with nonlinear regression. And some compounds showed in vivo efficacy in both mouse experimental autoimmune encephalomyelitis and collagen-induced arthritis models via oral administration (Wang et al., 2014).
To facilitate the drug discovery process, in silico modeling approaches (Hao et al., 2011;Wang, Li, Ma, Wang, & Wang, 2012;Wang et al., 2011) as a productive and cost-effective technology have been found to be valuable in further optimization and development of novel compounds. Quantitative structure-activity relationship (QSAR) is the most significant applications of chemometrics giving useful information for the design of novel inhibitors interacting with a specific target. Threedimensional quantitative structure-activity relationship (3D-QSAR) is a useful tool which connects the steric, electrostatic, hydrophobic, hydrogen-bond donor and hydrogen-bond acceptor structural requirements of compounds, and the biological activity. Molecular docking and molecular dynamics (MD) simulations, which can elucidate the interactions between the ligands and the receptor, and the results, can be analyzed to support the conclusions of 3D-QSAR models (Kirkpatrick, 2004).
In this study, mathematical models which are unique for this series of RORγt inhibitors were obtained through 3D-QSAR, molecular docking, MD simulations, and molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA). To gain the most reliable 3D-QSAR models, two different alignment rules were used. The reliability and robustness of the models were estimated by the internal and external validation. The 3D contour maps and molecular docking and MD simulations can validate each other, and the MM/PBSA method was employed to determine the major contributions to the binding free energy. Until recently, to our knowledge, this is the first 3D-QSAR study on this class of RORγt inhibitors, we hope the developed models can provide useful information for further design of novel inhibitors with improved inhibitory activity.

Methods and materials 2.1. Data sets and biological activity
The structures and the inhibitory activities of thiazole/thiophene ketone amides were taken from the previous literature (Wang et al., 2014), the structure and inhibitory activities of the RORγt inhibitors were used as the data set for molecular modeling in this study (Table S1). The corresponding pIC 50 values were used as dependent variable in CoMFA and CoMSIA model. The whole compounds were divided into a training set (including 37 compounds) for model generation and a test set (including 15 compounds) for model validation, respectively. The test set was selected by considering that the pIC 50 values are similar to that of the training set, and the structures cover as large diversity as the data set.
All molecular modeling and 3D-QSAR studies were performed using SYBYL package (Tripos Associates, St. Louis, MO). The Gasteiger-Hückel charges were added to the inhibitors. And the energy minimization and conformational search were performed by Tripos force field (Clark & Cramer, 1989) with the Powell conjugate gradient minimization algorithm. The energy minimization was employed to guarantee a low energy conformation with suitable bond distances and angles. The energy gradient convergence criterion was set to .05 kcal/mol Å in order to obtain the most stable conformation, and the maximum iterations were set to 1000.

Molecular alignment
Molecular alignment was considered as the most important step for 3D-QSAR (Thaimattam, Daga, Banerjee, & Iqbal, 2005). In the present work, two different alignment rules were employed to develop reliable CoMFA and CoMSIA models. The first alignment rule is ligand-based alignment (Superimposition I), during the model generation, the most potent inhibitor (compound 50) was selected as a template; all the inhibitors were then aligned to the common substructure of the template using the 'align database' function, which is shown in Figure 1(A) (marked in blue), and the result of superimposition is depicted in Figure 1(B); the second alignment rule is receptor-based alignment (Superimposition II) which is shown in Figure S1, the conformation of each compound obtained from molecular docking was aligned together to compound 50.

CoMFA and CoMSIA field calculation
The aligned compounds were placed in a 3D grid box with a grid spacing of 2 Å and extending 4 Å units beyond the aligned molecules in x, y, and z directions. Two descriptors: steric (Lennard-Jones potential) and electrostatic (Coulomb potential) field energies were calculated to build the CoMFA model , both steric and electrostatic interactions at each grid point were computed by applying the sp 3 carbon probe atom with a van der Waals radius of 1.52 Å and a charge of +1.0, the automatically generated CoMFA fields were scaled by the CoMFA-STD method.
CoMSIA method has several advantages over CoMFA such as greater robustness regarding both region shifts and small shifts within the alignments (Somayeh & Jahan, 2010). In CoMSIA, five different similarity fields: steric (S), electrostatic (E), hydrophobic (H), H-bond donor (D), and H-bond acceptor (A) were calculated. And the descriptors were calculated using a probe atom with radius 1.0 Å, +1.0 charge, +1 hydrophobicity, +1 H-bond donor, and +1 H-bond acceptor properties. A Gaussian method was applied to evaluate the distance between the probe atom and each molecule atom. In general, similarity indices, A F,K between the compounds of interest were computed by placing a probe atom at the intersections of the lattice points and using Equation (1):

Àar
(1) where q represents a grid point, i is the summation index over all atoms of the molecule j under computation, k represents five physicochemical properties (S, E, H, D, and A), W ik is the actual value of physicochemical property k of atom i, and W probe,k is the value of the probe atom, α is the attenuation factor with the default value of .3 used.
Partial least squares (PLS) method ) was used to analyze the training set. First, a molecular spreadsheet and PLS methods were employed to generate and evaluate the 3D-QSAR models, the cross-validation analysis was performed using the leave-one-out method in which one compound was removed from the data set, and its activity was predicted using the model derived from the rest of the data set. PLS was conjunct with cross-validation option to determine the optimum number of components and the lowest standard error of prediction, which was further employed to calculate the cross-validated correlated correlation coefficient R 2 cv . Second, the non-cross-validated analysis was performed to generate non-cross-validated correlation coefficient R 2 ncv . Additionally, the CoMFA and CoMSIA results were graphically represented by field contour maps using the field type 'StDev*Coeff'. In order to assess the robustness and statistical significance of the CoMFA and CoMSIA models, the test set which was not included in the training set was predicted. The predictive correlation coefficient (r pred 2 ) was defined as follows: where PRESS represents the sum of the squared deviations between predict and actual pIC 50 values for the test set compounds, SD represents the sum of the squared deviations between the biological activities of the test set molecules and the mean activity of the training set compounds.

Molecular docking
The crystal structure of RORγt (PDB entry code 4NIE) was extracted from Brookhaven protein database (http:// www.rscb.org/pdb). The binding mode of all inhibitors was obtained by a grid-based docking program Auto-Dock3.05 (Morris et al., 1998). The interaction energy between ligand and receptor was evaluated using atom affinity potentials calculated on a grid similar to that described by Goodford (1985). All ligands and water molecules were removed from the crystal structure, and polar hydrogen and united atom Kollman charges were assigned for the receptor. During the docking process, the relative ligand (compound NBH) in the crystal structure against 4NIE was used as a standard docking model. For the molecular docking, AutoGrid was used for the preparation of the grid map using a grid box with a number of points in xyz of 60 × 60 × 60 Å, the box spacing was .375 Å. And the grid center was designated at dimensions (xyz): 16.481 × 3.038 × −12.827. Finally, AutoDock was run using maximum number of 10,000 retries and 27,000 generations, respectively. The docking possibilities were calculated by the genetic algorithm with local search. For each inhibitor, the simulation was composed of 100 docking runs using the standard AutoDock parameters. During the docking simulation, inhibitors but not receptors were allowed to move. Then, the conformations with the lowest binding energy were extracted and aligned together for further QSAR analysis.

MD simulations
In order to identify a functionally validated complex from molecular docking, we also performed a 12 ns MD simulations to examine the best inhibitor on the receptor and to investigate the conformational changes. The MD simulations were performed with NAMD2.9 MD package (Phillips et al., 2005). Partial charges for the ligand were calculated with AM1 (Jakalian, Jack, & Bayly, 2002) in MOPAC7 and fitted into the AM1-bcc type of charge. The force field parameters for the ligand were assigned by the GAFF force field (Wang, Wolf, Caldwell, Kollman, & Case, 2004) using antechamber (Wang, Wang, Kollman, & Case, 2001) in AmberTools 13 (Case et al., 2012). Hydrogen atoms were added to the protein by use of the xleap module, and the protein was defined by the newly published and tested force filed, FF12SB. Two chloride counter ions were added to neutralize the system. The system was then solvated with a rectangular box of TIP3P (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983) water molecules. The minimum distance from the surface of each complex to the faces of the box was set to 12 Å. The solvated system contains approximately 41,256 atoms, 37,184 water molecules, and the final box dimension is 79.6 × 89.6 × 71.7 Å 3 . The SHAKE algorithm (Essmann, Ciccotti, & Berendsen, 1977) was applied to fix all bond lengths involving hydrogen bonds, permitting a 2 fs time step. The particle-mesh Ewald method (Essmann et al., 1995) was used to handle Coulombic interactions, and a 10 Å cutoff was applied on all van der Waals interactions. After 15,000 steps of minimization, an equilibration protocol started with heating the system to 300 K in 100 ps gradually, and then Langevin dynamics was switched on using constant pressure (1 atm) and fixed temperature (300 K) in another 100 ps. In the production phase, complex conformations were collected every 10 ps for 12 ns NPT MD simulation with constant pressure (1 atm) and isotropic molecule-based scaling. Finally, 200 snapshots from MD trajectories were sampled from the stable fluctuation region for postprocessing analysis and free energy calculations.

MM/PBSA binding free energy
The free energy of ligand binding in the binding pocket can be calculated by the MM/PBSA method (Kollman et al., 2000;Wang, Morin, Wang, & Kollman, 2001), which is summarized as: DG com , DG rec , and DG lig , stand for the free energies of the complex, the receptor, and the ligand, respectively. Each of them (DG com=rec=lig ) is equal to the enthalpy (DH) minus the change of conformational entropy (T DS), which was not considered here as there is no good way of estimating binding entropy without running extensive simulations using alchemical or similar approaches; thus, it is not necessary to calculate it especially when comparing the binding affinities of the similar compounds (Gohlke & Case, 2004) (Equation (4)). The enthalpy was calculated by aggregating the internal energy in gas phase (E gas ) and the solvation free energy (DG sol ) in Equation (5). E gas was the standard gas phase energy, including internal energy (E int ), van der Waals interactions (E vdw ), and electrostatic energies (E ele ) (Equation (6)). Because we performed a single molecular dynamics trajectory protocol, the contribution of internal energy term (E int ) to the binding energy is equal to zero. The solvation free energy, DG sol , was calculated with a PB/SA model, which decomposes solvation free energy by sum of a non-polar section (DG NP ) and an electrostatic section (DG PB=GB ) (Equation (7)). For the MM/ PBSA calculation, the DG PB part was calculated using the Poisson-Boltzmann function with the default cavity radii from the AMBER prmtop file. The dielectric constant was set to 1 for the interior solute and 80 for the surrounding solvent. While in terms of the MM/GBSA method, the alternativeDG PB part was replaced byDG GB , the Hawkins, Cramer, and Truhlar pairwise generalized Born model (Hawkins, Cramer, & Truhlar, 1996;Hou, Chen, McLaughlin, Lu, & Wang, 2006), and the parameters are described by Tsui and Case (Onufriev, Bashford, & Case, 2000). The LCPO approach (Weiser, Shenkin, & Still, 1999) was used to calculate the solvent accessible surface area (SASA) for the estimation of the nonpolar solvation free energy (DG NP ) using Equation (8) with g = .00542 kcal/mol Å 2 and b = .92 kcal/mol (Bringmann & Rummey, 2003).

QSAR model
To judge whether the QSAR model is robust and has a credible fitting for prediction of unknown inhibitors, a number of statistical parameters should be considered, that is, the cross-validated correlation coefficient (R 2 cv ), non-cross-validated correlation coefficient (R 2 ncv ), standard error of estimate (SEE), and F-statistic values (F), etc. CoMFA model calculated steric and electrostatic properties. In addition to steric and electrostatic fields, hydrophobic, hydrogen-bond donor and acceptor descriptors were calculated in CoMSIA models. It has been argued that the five different descriptor fields may not be totally independent of each other, and such dependencies of the individual fields usually decrease the statistical significance of the results (Bohm, Stürzebecher, & Klebe, 1999;Bringmann & Rummey, 2003), thus in the present work, all combinations of the five fields were calculated to choose the suitable CoMSIA models (shown in Table S2 and Table S3). The statistic results of the CoMFA and CoMSIA models applying two different alignments are summarized in Table 1.
Since the QSAR models are sensitive to the alignment rules, differences in the statistical values are observed with different superimpositions. As shown in Table 1, the results of both CoMFA and CoMSIA models based on superimposition I are better than those obtained from superimposition II. Thus, in the present work, superimposition I model will be discussed in the following sections. The actual and predicted pIC 50 of these compounds are shown in Table S4 and Table S5.
For CoMFA model, the PLS analysis gives a relatively high R 2 cv value of .859, R 2 ncv value of .993, F value of 689.960, and a SEE value of .120 at 6 components, suggesting that this model is a useful tool for predicting the related inhibitory activities. The CoMFA model includes 75.7% of the steric field descriptor, but only 24.3% of the electrostatic filed, thus, the steric field has greater influence than the electrostatic field.
For CoMSIA model, 31 possible descriptors' combinations were calculated (shown in Table S2). Finally, the model using electrostatic field and H-bond donor field is superior to the other models. The best CoMSIA model (ED) gives an R 2 cv of .805 with an optimized component of 4, which suggests that this model should be considerably reliable to predict the IC 50 values. The non-crossvalidated PLS analysis produces an R 2 ncv of .920, F value of 92.548, and a low SEE of .386. Contributions of E and D fields are 89.1% and 10.9%, respectively, which indicates that the D field makes less contribution to the ligand-binding affinity.
A test set of 15 inhibitors excluded from the model derivation was employed to test the predictive ability of the model. The predictive correlation coefficient R 2 pred is .7317 and .7097 for CoMFA and CoMSIA models, respectively. The graph of actual versus predicted pIC 50 values of the training set and test set is illustrated in Figure 2, the plots represent a uniform distribution around the regression line with respective slope and intercept very close to one and zero, indicating the satisfactory predictive capability and accuracy of the model.
Additionally, compound 39 is detected as outlier, outliers are inhibitors those do not fit the derived model Notes: R 2 cv = cross-validated correlation coefficient using the leaveone-out methods. R 2 ncv = non-cross-validated correlation coefficient; SEE = standard error of estimate; F = ratio of R 2 ncv explained to unexplained = R 2 ncv / (1 − R 2 ncv ). R 2 pred = predicted correlation coefficient for the test set of compounds; SEP = standard error of prediction; N C = optimal number of principal components; S = steric, E = electrostatic, H = hydrophobic, D = H-bond donor, A = H-bond acceptor. Superimposition Ⅰ, ligand-based alignment; Superimposition II, receptor-based alignment. or are poorly predicted by it due to the low inhibitory activity (Egan & Morgan, 1998). By removal of the outlier, stronger, and more significant CoMFA and CoMSIA models are produced (R 2 pred of .9198 and .9026 for CoMFA and CoMSIA, respectively).

3D-QSAR contour maps
CoMFA and CoMSIA contour maps were generated to visualize the 3D-QSAR model and rationalize the 3D spatial regions around the inhibitors, which show regions where differences in molecular fields are associated with differences in the inhibitory activity (Mao, Li, Hao, Zhang, & Ai, 2012). To aid in visualization, the most potent inhibitor 50 is shown in the contour maps as depicted in Figures 3 and 4.

CoMFA contour maps
The CoMFA steric and electrostatic field contour maps with the most potent inhibitor 50 as a reference are shown in Figure 3. In the CoMFA steric contour (Figure 3(A)), green maps represent areas where bulky groups are favored, while the yellow contour maps indicate areas where bulky groups are disfavored. A large green contour map is found over the 4-position of ring A, suggesting that inhibitors with bulky groups at this position should be more active than those with no or smaller groups. The fact that compound 42 (having -Cl at this position) is more potent in activity than compound 47 (with -H at the same location) is a good example. Around the 3-Cl of ring B, it is noted that there is a small yellow-colored map indicating a small group is advantageous to the inhibitory activity. This is in line with the fact that compound 38 with an -CH 2 NMe 2 substituent shows less potency than its counterpart compound 37 with a -CN group. Additionally, a small yellow contour map is appeared at the 2-F of ring A, indicating that its preference for the steric moderate and less crowed groups in this area. This may explain why compound 50 (-F) exhibits higher inhibitory activity than compound 48 (-Cl). It should be noted that some yellow contours are located around Region A (shown in Figure 5), suggesting that substituent with the bulky steric is detrimental to the inhibitory activity. By comparing compounds 10, 11 and 12 and compounds 8 and 9, it is found that the presence of a larger group is harmful to the activity.
For the electrostatic field, blue contours indicate regions where electropositive groups will enhance the inhibitory activity, and red contours represent regions where electronegative groups are favorable. A blue contour map is located above ring A, implying electropositive groups at this position will increase the inhibitory activity. This may explain why the activity of compound 31  In addition, if a benzene ring is connected to other groups such as -Cl, the density of the electron cloud near -Cl will be increased and will be impaired around benzene ring, this could be shown by the order for the activity of compounds, such as activity of compound 50 (-Cl) > 46 (-H). Another large blue contour map around ring C indicates that the electron-donating substituent can enhance the activity. This is in good correlation with the experimental activity distribution of the inhibitors. By comparing the inhibitory activities, the rank would be: 7 > 6 > 5 (the -SO 2 can be in resonance with the benzene ring, thus the electron-withdrawing tendency of -SO 2 is higher than -CN and -Cl). On the other hand, most of molecules have electronegative substituent at Region A, suggesting that there is room for modification in this area. Two red contours located around the linker between ring A and ring D indicate that electronegative groups in this region would benefit the potency. This could be explained by the fact that the activity of compound 28 (possessing -O at this position) is higher than compound 27 (possessing -C at this position).

CoMSIA contour maps
For the CoMSIA, an electrostatic field contour map is depicted in Figure 4(A). Similar results are obtained for the electrostatic contours of CoMSIA as those of the CoMFA model, with a red contour located at the linker between ring A ring D, and a blue contour around ring C. However, a different large blue contour map is depicted around ring A, indicating that more positively charged group situated in this area will enhance the biological activity, which can be explained by the finding that the inhibitory activity of compound 50 (X=C) is higher than compound 43 (X=N). In fact, this finding is fully consistent with observations obtained from ref 14, which systematically compared the effect of changing the substitutes in this position by experiments and found that if the -C was changed to -N, the inhibitory activity will be reduced. The H-bond donor maps are shown in cyan and purple (Figure 4(B)), indicating areas where H-bond donor groups would increase and decrease the inhibitory activity. A cyan contour is found near the linker between ring C and ring D which suggests that H-bond donor groups at this position would enhance the inhibitory activity by forming hydrogen bonds with the receptor. In fact, all of the molecules bear -NH groups at this position. Additionally, a purple contour near the cyan contour indicates that the presence of strong H-bond acceptor groups extending into this contour is beneficial to the biological activity. In addition, another large purple contour at the end of Region A depicts the unfavorable interactions with the H-bond donor groups. In the most potent compound 50, the -SO 2 group served as the H-bond acceptor is favorable for the activity and its replacement of them would lead to lower activity such as compounds 1-6.

Docking analysis and comparison with 3D contour maps
Molecular docking helps us to study the ligand-receptor interaction and obtain the best geometry of ligand-receptor complex so that the energy of interaction between ligand and receptor is minimum (Bhadoriya, Sharma, Jain, Raut, & Rananaware, 2013). All inhibitors in the data set were docked into the active site of RORγt. To illustrate the interaction mechanism between the inhibitor and the receptor, attention was focused on the most active compound 50, which is shown in Figure 5. To check whether the parameters used in the procedure of docking are reasonable, the docking mode of RORγt -50 is superimposed on the crystal structure of 4NIE, as shown in Figure S2, which reveals that the two compounds locate in the same binding site and have the similar orientation. The active site mainly consists of hydrophobic residues (Leu287, Ala321, Ala327, Ala368, Phe377, Phe378, Phe388, Val361, Met365, and Ile400) and electropositive residues (Arg364 and Arg367) (Figure 6(A)). To evaluate the docking results, H-bond as the most widely used parameters was calculated (Figure 6(B)). The -NH located between ring C and ring D acts as a donor to form hydrogen bond with the side chain of Phe377 (-O···HN, 2.02 Å, 162.2°) (H-1), thus, the H-bond interactions at this position may favor the binding activity which is shown by the cyan contour at the H-bond contour maps. Another H-bond interactions between -SO 2 at Region A and the backbone of Arg367 is observed (-O···HN, 1.81 Å, 142.9°) (H-2), (-O···HN, 2.29 Å, 115.2°) (H-3). In addition, we can find that -SO 2 also forms a hydrogen bond with the -NH on the residue Leu287 (-O···HN, 1.89 Å, 147.9°) (H-4), which is consistent with the purple plot indicating activity-favorable H-bond acceptor area at the corresponding position on the contour maps of the CoMFA model and the CoMSIA model.
The binding pocket of RORγt can be divided into two sub-pockets which are marked as P1 and P2 (shown in Figure 7). The P1 pocket is formed by residues Cys285, Gln286, Leu287, Arg364, Arg367, and Ala368, thus is positively charged, analysis indicates that electropositive residues Arg364 and Arg367 are significant for the inhibitor binding. Docking results demonstrate that Region A is placed into P1 pocket where the volume is small, thus suggesting that too large groups may introduce severe steric clash. This finding is also inferred from the steric yellow maps of the CoMFA model. A large hydrophobic cavity (P2 sub-pocket) is formed among Cys320, Ala321, Phe388, and Ile400, which accommodates ring A and ring B, thus a large substituent, such as the benzene ring, in this position will increase the activity. This is fully supported by the CoMFA steric maps results.
The comparison of the docking results and QSAR models indicates that the QSAR models developed are reasonable and can offer constructive suggestions to further modification of RORγt inhibitors.

Comparison of 3D-QSAR, molecular docking and molecular dynamics simulations
In order to further validate the results of 3D-QSAR and molecular docking, MD simulations were applied to predict more reliable ligand-receptor interaction mechanisms. The 12 ns MD simulations in aqueous explicit solution were performed to observe the dynamic picture of the conformational changes for 4NIE-50. To explore the dynamic stability of the complexes, the root-meansquare deviations (RMSDs) of complex and ligand of the MD trajectories were calculated respectively (shown in Figure 8 complex sharply rise to 2.2 Å in the first .8 ns and maintain stable at the range of 1.8-2.1 Å ever since then, which indicates that the whole system is in equilibrating state from .8 to 12 ns. As expected, the RMSD of the ligand stabilizes at lower window (.75 Å), however, it vibrates in a larger range in .5-1.5 Å. This is because we took all atoms (including hydrogen) in the RMSD calculation. The ligand RMSD also witness a drastic conformational variation at 5 ns, the jump is then stabilized further between 7 ns and 12 ns, which is confirmed by the superposition of coordinates of the protein after MD simulation onto the initial conformation without significant fluctuation (shown in Figure 8(C)). In addition, the RMSD for residues at the ligand-receptor interface is stabilized in 1.5 Å, analysis indicates that these residues include Cys285, Gln286,Leu287,Leu292,Cys320,Ala321,His323,Leu324,Met358,Arg364,Met365,Arg367,Ala368,Val376,Phe377,Phe378,Phe388,Leu391,Cys393,Leu396,Ile397,Ile398,Phe399,and His479. The RMSD plots indicate that the tendency of the RMSD of residues in the binding pocket copes well with the movement of the ligand; furthermore, the RMSD plot of residues within 5 Å of the ligand is lower than that of the overall backbone RMSD, this means that the ligand has taken up the entire pocket space and banded firmly to the receptor.
In order to validate the convergence of the dynamic attribution of the system, we also compared the average B-factor obtained from the all-atom MD simulation and from the crystal structure (4NIE.pdb), shown in Figure 8(B). Most regions of the protein show great agreement with the X-ray structure. Interestingly, the simulated results yield several spectacular peaks observed in the crystallographic structure around residues 264-276 and residues 469-508, all these residues with higher fluctuation values are found belonging to the highly mobile solvent-exposed amino acids. This result suggests that these residues play significant roles in the protein movements, which are hard to obtain in the low temperature X-ray experiment. Moreover, we find those residues in the ligand-binding domain possess less fluctuation and low degree of flexibility, which can also be illustrated by the RMSD of residues in the ligandbinding pocket. Therefore, we speculate that the ligand interacts with the receptor steadily.
From Figure 6(B) and Figure 8(D), we can find that the complex has undergone several changes during MD simulation, for example, the distance between the -NH group (located between ring A and ring C) and Phe377, the -SO 2 at Region A and Arg367, the -SO 2 and the backbone of Leu287 is changed from 2.02 to 2.05 Å, 1.81 to 2.33 Å, 2.29 to 3.03 Å and 1.89 to 3.17 Å after MD simulation; however, the binding pocket and the conformation of the ligand are still stable, suggesting the rationality and validity of the molecular docking model. The binding affinity of compound 50 with RORγt was further calculated using MM/PBSA method. The results are shown in Table 2. An inspection of the free energy components of the complexes reveals that the van der Waals and the electrostatic interactions play vital roles in the ligand binding for the complex. The non-polar solvation interactions, composing of the van der Waals interaction and the non-polar solvation effect, are the dominant driving forces for the burial of inhibitors' hydrophobic groups on binding. The polar interaction is −20.38 and 42.38 kcal/mol for DE ele and DG PB , respectively, suggesting that electrostatic interaction (favorable) cannot overcome the larger polar energy barrier. For the results of molecular docking and MD simulation, we can see that the binding pocket is mainly hydrophobic, which is also suggested by the MM/PBSA calculation, where the value of the DE vdw (−66 kcal/mol) is greater than that of the DE ele (−20.38 kcal/mol).

Conclusions
In the present work, a series of RORγt inhibitors have been estimated by developing 3D-QSAR models based on ligand-and receptor-based superimpositions, the 3D-QSAR models based on ligand exhibit better statistical results than those based on receptor does, considering several evaluation criteria, such as R 2 cv , R 2 ncv , SEE, F and R 2 pred . The 3D contour plots provide useful information about the intermolecular interactions of the inhibitors with the surrounding environment. In addition, the 3D-QSAR models, molecular docking and MD simulations exhibit a good consistency indicating that the derived models are reliable and robust. The MD simulation results indicate that compound 50 uses one hydrogen atom in the amino group as hydrogen-bond donor and three oxygen atoms in the sulfur dioxide group as hydrogen-bond acceptors in interactions with Phe377, Arg367, and Leu287. The binding free energy calculated by MM/PBSA method has a good correlation with the ligand-receptor interaction mode, and suggesting that the major contributions to the binding process are van der Waals contacts and electrostatic interactions.
According to the detailed contour analyses of the models, some useful information on the structural requirements for the observed inhibitory activity is obtained as follows: (1) The minor and electropositive groups around ring A have a major impact on the inhibitory activity.
(2) The electronegative and H-bond acceptor groups located between ring A and ring D might increase the activity. (3) Substituents with electropositive groups above ring C are favored. (4) The introduction of large substituents at the 3-Cl of ring B would enhance the ligand binding. (5) Minor and H-bond acceptor groups near the end of Region A could result in larger binding affinity.
Overall, these findings will aid in a better structural understanding of inhibitors targeting RORγt and provide an insight into some instructions for further design of novel RORγt inhibitors.