Molecular Modeling Study for the Design of Novel Peroxisome Proliferator-Activated Receptor Gamma Agonists Using 3D-QSAR and Molecular Docking

Type 2 diabetes is becoming a global pandemic disease. As an important target for the generation and development of diabetes mellitus, peroxisome proliferator-activated receptor γ (PPARγ) has been widely studied. PPARγ agonists have been designed as potential anti-diabetic agents. The advanced development of PPARγ agonists represents a valuable research tool for diabetes therapy. To explore the structural requirements of PPARγ agonists, three-dimensional quantitative structure–activity relationship (3D-QSAR) and molecular docking studies were performed on a series of N-benzylbenzamide derivatives employing comparative molecular field analysis (CoMFA), comparative molecular similarity indices analysis (CoMSIA), and surflex-dock techniques. The generated models of CoMFA and CoMSIA exhibited a high cross-validation coefficient (q2) of 0.75 and 0.551, and a non-cross-validation coefficient (r2) of 0.958 and 0.912, respectively. The predictive ability of the models was validated using external validation with predictive factor (r2pred) of 0.722 and 0.682, respectively. These results indicate that the model has high statistical reliability and good predictive power. The probable binding modes of the best active compounds with PPARγ active site were analyzed, and the residues His323, Tyr473, Ser289 and Ser342 were found to have hydrogen bond interactions. Based on the analysis of molecular docking results, and the 3D contour maps generated from CoMFA and CoMSIA models, the key structural features of PPARγ agonists responsible for biological activity could be determined, and several new molecules, with potentially higher predicted activity, were designed thereafter. This work may provide valuable information in further optimization of N-benzylbenzamide derivatives as PPARγ agonists.


Introduction
Type 2 diabetes (T2D) is a disease that is generally characterized by relative insulin deficiency caused by insulin resistance in target organs, and pancreatic β-cell dysfunction [1].In 2014, there were 422-million people with diabetes, with more than 90% estimated to have T2D, worldwide.Unfortunately, this number will increase to approximately 552-million by the year 2030 [2].Accordingly, T2D is generating a significant socioeconomic burden, as a pandemic disease with a high and increasing fatality [3,4].
The peroxisome proliferator-activated receptor γ (PPARγ) is generally regarded as a molecular target for the thiazolidinedione class of anti-diabetic drugs [5,6], as it plays a key role in the generation and development of diabetes mellitus [7][8][9].Recent studies have shown that PPARγ agonists, including rosiglitazone and pioglitazone [10], may be used as insulin sensitizers in target tissues to lower glucose, as well as fatty acid levels in T2D patients.
However, both rosiglitazone and pioglitazone have been withdrawn from the market because of significant hepatotoxicity and cancer development concerns [11].Hence, there is an urgent need for the development of safer PPARγ modulating drugs.One severe side-effect of known PPARγ agonists, involves sodium and water retention, which may be dangerous for patients suffering from congestive heart conditions [12].Recently, various new N-benzylbenzamide compounds have been shown to act as PPARγ agonists that, not only lowered blood pressure and reduced systemic glucose, triglycerides, and free fatty acid levels, but have also been shown to maintain water and electrolyte homeostasis [13].Therefore, a variety of N-benzylbenzamide compounds have since been identified as safer PPARγ modulators for the treatment of T2D.
Based on CoMFA [14], along with CoMSIA [15], methods involving 3D-QSAR determinations allow for the structure-activity relationship of N-benzylbenzamide compounds to be studied.Molecular docking was also applied to reveal the most likely binding modes between the compounds and PPARγ.On the basis of 3D-QSAR and molecular docking results, valuable information can be retrieved for further structured-based drug design, with higher activity.Finally, a series of new potent molecules with a higher predicted activity than the template compound, the latter exhibiting the best activity reported in the literature, have been designed.Our study will potentially provide guidance for the future design of selective and potent PPARγ agonists.

CoMFA and CoMSIA Results
The 3D-QSAR models were obtained using a training set of 27 compounds, and a test set of six compounds.The statistical parameters associated with CoMFA and CoMSIA can be found in Table 1.In general, various alignment strategies can lead to different statistical values in the constructed QSAR models.The best CoMFA and CoMSIA models were generated employing a partial least square (PLS) analysis, which produced cross-validated coefficients (q 2 ).When a cross-validation coefficient, q 2 > 0.5, was used, the QSAR model demonstrated statistical significance.
As shown in Table 1, two descriptor fields in CoMFA form all three possible combination models, including steric (S), electrostatic (E) and SE models.The CoMSIA models, with a combination of five descriptor fields, including S, E, hydrophobic (H), hydrogen bond donor (D) and acceptor (A), were developed to generate the optimal 3D-QSAR model.However, some models with a low q 2 value did not meet the criterion (q 2 > 0.5), indicating an unacceptable 3D-QSAR model.Still, overfitting seemed to occur for some models (those with a large number of components).From Table 1, we can see that the best established models (CoMFA and CoMSIA) exhibited high q 2 (0.75 and 0.551), r 2 (0.958 and 0.912), and F-values (76.113 and 43.388), along with a low standard error of estimate (SEE) value (0.097 and 0.138), and a suitable number of components (6 and 5), which indicated good statistical correlation of the models.Moreover, the predictive capabilities of the generated models were assessed by calculating their predictive correlation coefficient (r 2 pred ) involving their corresponding test set molecules.The generated CoMFA and CoMSIA models with maximum external predictive ability (r 2 pred 0.722 and 0.682), were considered the best models.The distribution of actual predicted pEC 50 values of the training and test sets for CoMFA and CoMSIA are shown in Figure 1.The CoMFA and CoMSIA models show a good fit along the diagonal line.Both models also exhibited satisfactory predictive ability throughout the training and test sets.

CoMFA Contour Map Analysis
The steric and electrostatic fields of the CoMFA model are presented as contour maps in Figure 2. Finally, compound 24c was selected as the template molecule.The green contours represent regions indicating favorable steric fields, while the yellow contours represent the regions indicating unfavorable steric fields.Moreover, the blue and red contours highlight the positions where

CoMFA Contour Map Analysis
The steric and electrostatic fields of the CoMFA model are presented as contour maps in Figure 2. Finally, compound 24c was selected as the template molecule.The green contours represent regions indicating favorable steric fields, while the yellow contours represent the regions indicating unfavorable steric fields.Moreover, the blue and red contours highlight the positions where electropositive groups and electronegative groups would be favorable, respectively.

Steric Contour Map
The steric contour map in CoMFA (Figure 2A) has a yellow contour near the ortho position of the benzene ring, which indicates that the presence of steric substituents in this region is unfavorable.Furthermore, the yellow contour explains why a -CF 3 substituent in ortho position of the benzene ring in compound 2b is more potent than in compound 10b, which bears a -OCF 3 substituent.Likewise, a small yellow contour map appeared in the para position of the benzene ring, which indicates that the large size of the substituent was not preferred in this area.Moreover, a -OCH 3 group in this position in compound 18b, could be found within the steric field, which led to decreased biological activity.Finally, the large yellow region on the R 2 substitutes may explain why compound 32b, bearing a phenyl group, was less active than compound 30b bearing a propyl group.

CoMFA Contour Map Analysis
The steric and electrostatic fields of the CoMFA model are presented as contour maps in Figure 2. Finally, compound 24c was selected as the template molecule.The green contours represent regions indicating favorable steric fields, while the yellow contours represent the regions indicating unfavorable steric fields.Moreover, the blue and red contours highlight the positions where electropositive groups and electronegative groups would be favorable, respectively.

Steric Contour Map
The steric contour map in CoMFA (Figure 2A) has a yellow contour near the ortho position of the benzene ring, which indicates that the presence of steric substituents in this region is unfavorable.Furthermore, the yellow contour explains why a -CF3 substituent in ortho position of the benzene ring in compound 2b is more potent than in compound 10b, which bears a -OCF3 substituent.Likewise, a small yellow contour map appeared in the para position of the benzene ring, which indicates that the large size of the substituent was not preferred in this area.Moreover, a -OCH3 group in this position in compound 18b, could be found within the steric field, which led to decreased biological activity.Finally, the large yellow region on the R2 substitutes may explain why compound 32b, bearing a phenyl group, was less active than compound 30b bearing a propyl group.

Electrostatic Contour Map
A large blue contour area near the para position of the benzene ring indicates that the presence of an electropositive group may increase activity (Figure 2B).This assumption becomes even more significant in the case of compounds 20b and 12b, as these compounds contain a -Cl and -F substituent, respectively.However, due to the presence of different electron-donating groups,

Electrostatic Contour Map
A large blue contour area near the para position of the benzene ring indicates that the presence of an electropositive group may increase activity (Figure 2B).This assumption becomes even more significant in the case of compounds 20b and 12b, as these compounds contain a -Cl and -F substituent, respectively.However, due to the presence of different electron-donating groups, compound 12b was found to be less biologically active than compound 20b.The red contours present on the ortho position of the benzene ring suggest that an electron negative group would be favorable in this area, an assumption that proves to be true for compounds 2b and 6b, which contain a -CF 3 group and a -CH 3 group, respectively.However, since a -CF 3 group proves to be more electron-withdrawing than a -CH 3 group, compound 2b was determined to be more biologically active than compound 6b.

CoMSIA Contour Map Analysis
The CoMSIA steric and electrostatic contour maps were both similar to the CoMFA contour maps discussed above (Figure 3A,B).Thus, only hydrophobic, hydrogen bond donor, as well as hydrogen bond acceptor fields of CoMSIA, were analyzed in this section.The CoMSIA steric, electrostatic, hydrophobic, hydrogen bond donor along with hydrogen bond acceptor contour maps are shown in Figure 3, respectively.Compound 24c was selected as the corresponding reference molecule.
The CoMSIA steric and electrostatic contour maps were both similar to the CoMFA contour maps discussed above (Figure 3A,B).Thus, only hydrophobic, hydrogen bond donor, as well as hydrogen bond acceptor fields of CoMSIA, were analyzed in this section.The CoMSIA steric, electrostatic, hydrophobic, hydrogen bond donor along with hydrogen bond acceptor contour maps are shown in Figure 3, respectively.Compound 24c was selected as the corresponding reference molecule.

Hydrophobic Contour Map
In the hydrophobic contour map (Figure 3C), the yellow contours show favorable hydrophobic regions, while white contours represent unfavorable hydrophobic regions.For the hydrophobic map,

Hydrophobic Contour Map
In the hydrophobic contour map (Figure 3C), the yellow contours show favorable hydrophobic regions, while white contours represent unfavorable hydrophobic regions.For the hydrophobic map, one white unfavorable region could be found around the R 2 substitutes, indicating that the addition of hydrophobic substituents in this region would lead to decrease in activity.Further evidence for this notion can be obtained from compound 30b bearing a propyl group, which is considerably less hydrophobic than a phenyl group.Therefore, compound 30b proves to be more active than the biologically less active compound 32b, bearing a phenyl group.The other white contour area observed in the para position of the benzene ring, indicates that hydrophobic substituents were not preferred in this region.This finding can be further explained by the fact that compound 21c contains a -Cl substituent that generally leads to a higher potency compared to a methoxy group present in compound 19c.

Hydrogen Bond Donor Map
The contour map for the hydrogen bond donor field is shown in Figure 3D.Cyan and purple contours represent a hydrogen bond donor field favorable region and hydrogen bond donor unfavorable region, respectively.For the hydrogen bond donor map, a cyan contour appeared around the hydroxyl group.This suggests that a hydrogen bond interaction, with the hydrogen atom of the hydroxyl group acting as a hydrogen bond donor, is favorable for increased activity.

Hydrogen Bond Acceptor Map
In the hydrogen bond acceptor contour map (Figure 3E), the magenta contours represent a favorable hydrogen bond acceptor field, while the red contours represent an unfavorable hydrogen bond acceptor field.For the hydrogen bond acceptor map, one favorable polyhedral surface (magenta) is found around the carboxyl group, which suggests that hydrogen bond interactions between the oxygen atom of the carbonyl group, acting as a hydrogen bond acceptor, and a hydrogen atom of the group, lead to an increase in activity.

Design of More Potent Compounds
Based on CoMFA and CoMSIA models obtained in the present study, the structure-activity relationships of PPARγ agonists could be determined, and several new potent molecules could be designed.The chemical structures of the newly designed compounds, as well as their activity characteristics on PPARγ, were predicted by the CoMFA and CoMSIA models, as seen in Table 2.The predicted activities of the newly designed compounds on PPARγ were all significant.A set of the molecules demonstrated an even better activity than the most active agonist previously reported, further validating the superiority of the models, and indicates that the structure-activity relationships in the work reported herein, may potentially be used in structural modification and optimization.a -Cl substituent that generally leads to a higher potency compared to a methoxy group present in compound 19c.

Hydrogen Bond Donor Map
The contour map for the hydrogen bond donor field is shown in Figure 3D.Cyan and purple contours represent a hydrogen bond donor field favorable region and hydrogen bond donor unfavorable region, respectively.For the hydrogen bond donor map, a cyan contour appeared around the hydroxyl group.This suggests that a hydrogen bond interaction, with the hydrogen atom of the hydroxyl group acting as a hydrogen bond donor, is favorable for increased activity.

Hydrogen Bond Acceptor Map
In the hydrogen bond acceptor contour map (Figure 3E), the magenta contours represent a favorable hydrogen bond acceptor field, while the red contours represent an unfavorable hydrogen bond acceptor field.For the hydrogen bond acceptor map, one favorable polyhedral surface (magenta) is found around the carboxyl group, which suggests that hydrogen bond interactions between the oxygen atom of the carbonyl group, acting as a hydrogen bond acceptor, and a hydrogen atom of the group, lead to an increase in activity.

Design of More Potent Compounds
Based on CoMFA and CoMSIA models obtained in the present study, the structure-activity relationships of PPARγ agonists could be determined, and several new potent molecules could be designed.The chemical structures of the newly designed compounds, as well as their activity characteristics on PPARγ, were predicted by the CoMFA and CoMSIA models, as seen in Table 2.The predicted activities of the newly designed compounds on PPARγ were all significant.A set of the molecules demonstrated an even better activity than the most active agonist previously reported, further validating the superiority of the models, and indicates that the structure-activity relationships in the work reported herein, may potentially be used in structural modification and optimization.

Docking Analysis
In order to obtain the probable binding conformations between the molecules and the protein, Surflex-dock was carried out to dock the compounds to the binding site of PPARγ.In this study, compound 24c (template) and a newly designed compound, N1, N9, and N12, were placed in the corresponding binding sites, respectively.The docking score of compound 24c was 8.913.Meanwhile, the docking scores of compounds N1, N19, and N12 were 11.0573, 11.010, and 11.690, respectively.The docking scores of compounds N1, N9, and N12 are higher than compound 24c, which has the highest activity in the training set.This result is in good agreement with corresponding predicted activities of CoMFA and CoMSIA models.
Figure 4 shows the surface of the binding site of PPARγ, the binding modes between compounds 24c, N1, N9, and N12, and the binding site of the protein.Figure 4A-D, illustrate the surface of the binding site and the conformations of compounds 24c, N1, N9, and N12 (yellow), and the original ligand (purple), as well as the key residues (white) at the binding site.High resemblance between these molecules is observed and they occupied nearly the same binding pocket as PPARγ.It is representative of the active conformation to dock selected compounds.Here, compounds were positioned in the pocket, surrounded by His323, Tyr473, Ser289, Leu453, Ile341, Cys285, etc.As seen in Figure 4E-H, the carbonyl group of the ligand that acts as a hydrogen bond acceptor, formed a hydrogen bond with the backbone O-H of Ser289.The backbone N-H of His323, and O-H of Tyr473, which act as hydrogen bond acceptors and formed hydrogen bonds with the hydroxyl group of the ligand, respectively.Thus, these three hydrogen bond interactions played a major role in the combination of these drugs and the receptor.

Docking Analysis
In order to obtain the probable binding conformations between the molecules and the protein, Surflex-dock was carried out to dock the compounds to the binding site of PPARγ.In this study, compound 24c (template) and a newly designed compound, N1, N9, and N12, were placed in the corresponding binding sites, respectively.The docking score of compound 24c was 8.913.Meanwhile, the docking scores of compounds N1, N19, and N12 were 11.0573, 11.010, and 11.690, respectively.The docking scores of compounds N1, N9, and N12 are higher than compound 24c, which has the highest activity in the training set.This result is in good agreement with corresponding predicted activities of CoMFA and CoMSIA models.
Figure 4 shows the surface of the binding site of PPARγ, the binding modes between compounds 24c, N1, N9, and N12, and the binding site of the protein.Figure 4A-D, illustrate the surface of the binding site and the conformations of compounds 24c, N1, N9, and N12 (yellow), and the original ligand (purple), as well as the key residues (white) at the binding site.High resemblance between these molecules is observed and they occupied nearly the same binding pocket as PPARγ.It is representative of the active conformation to dock selected compounds.Here, compounds were positioned in the pocket, surrounded by His323, Tyr473, Ser289, Leu453, Ile341, Cys285, etc.As seen in Figure 4E-H, the carbonyl group of the ligand that acts as a hydrogen bond acceptor, formed a hydrogen bond with the backbone O-H of Ser289.The backbone N-H of His323, and O-H of Tyr473, which act as hydrogen bond acceptors and formed hydrogen bonds with the hydroxyl group of the ligand, respectively.Thus, these three hydrogen bond interactions played a major role in the combination of these drugs and the receptor.

Data Set
A set of 33 N-benzylbenzamide derivatives and the corresponding activity data were collected from the work of René Blocher et al. [13] (Table 3).The data set was randomly divided into the training set of 27 compounds (82%) to generate the 3D-QSAR model, and the test set of 6 compounds (18%) to verify the predictive ability of the model.The bioactivities of the compounds were expressed as pEC50 (-logEC50), which was used as a dependent variable in further investigations.To avoid possible issues during the external validation, the selection of the training and the test set was carried out such that both sets included structurally diverse compounds and all types of activity [16][17][18].

Data Set
A set of 33 N-benzylbenzamide derivatives and the corresponding activity data were collected from the work of René Blöcher et al. [13] (Table 3).The data set was randomly divided into the training set of 27 compounds (82%) to generate the 3D-QSAR model, and the test set of 6 compounds (18%) to verify the predictive ability of the model.The bioactivities of the compounds were expressed as pEC 50 (-logEC 50 ), which was used as a dependent variable in further investigations.To avoid possible issues during the external validation, the selection of the training and the test set was carried out such that both sets included structurally diverse compounds and all types of activity [16][17][18].

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with 1000 iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the lowest energy conformation of each molecule was determined for the definitive QSAR studies.Molecular alignment is one of the most essential steps for the generation of the best CoMFA and CoMSIA models [21].Thus, molecular alignment was performed using the Distill alignment technique, a user-defined common core of the Sybyl tools [22].Compound 24c, exhibiting the highest activity in the complete data set, was selected as the template molecule.The remaining compounds in the Mol2 database were aligned by their corresponding maximum common substructures, as shown in Figure 5A.The rigid body alignment of the molecules is shown in Figure 5B.

Molecular Modeling and Alignment
To obtain the best conformers for each molecule, the Sybyl X-2.1.1 software package was used for all compound modeling and optimization parameters.All structures of the compound series were subjected to preliminary geometry optimization using the Tripos force field with 1000 iterations [19].Partial atomic charges were calculated by the Gasteiger-Hückel scheme, with an energy gradient convergence criterion of 0.05 kcal/mol Å [20].Based on the analysis method described above, the lowest energy conformation of each molecule was determined for the definitive QSAR studies.Molecular alignment is one of the most essential steps for the generation of the best CoMFA and CoMSIA models [21].Thus, molecular alignment was performed using the Distill alignment technique, a user-defined common core of the Sybyl tools [22].Compound 24c, exhibiting the highest activity in the complete data set, was selected as the template molecule.The remaining compounds in the Mol2 database were aligned by their corresponding maximum common substructures, as shown in Figure 5A.The rigid body alignment of the molecules is shown in Figure 5B.

CoMFA Method
The CoMFA method is often used to describe steric and electrostatic fields.Lennard-Jones and Coulomb potentials were employed to calculate two fields.A 3D cubic lattice, with grid spacing of 2.0 Å, was generated to surround the aligned molecules in all directions.These grid points were generated using the Tripos force field, a sp 3 carbon atom probe with a Van der Waals radius of 1.52 Å, and a charge of +1.00 (default probe atom in Sybyl).Based on the CoMFA method, steric and

CoMFA Method
The CoMFA method is often used to describe steric and electrostatic fields.Lennard-Jones and Coulomb potentials were employed to calculate two fields.A 3D cubic lattice, with grid spacing of 2.0 Å, was generated to surround the aligned molecules in all directions.These grid points were generated using the Tripos force field, a sp 3 carbon atom probe with a Van der Waals radius of 1.52 Å, and a charge of +1.00 (default probe atom in Sybyl).Based on the CoMFA method, steric and electrostatic fields were scaled with a default energy cut off of 30 kcal/mol, the latter being the optimal parameter for this model [23].

CoMSIA Method
The CoMSIA analysis is similar to CoMFA, in regard to the descriptors around the aligned molecules.Three other fields, (i.e., hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields), were calculated together with the same standard settings used in the CoMFA calculations.More importantly, the distance dependence between the probe atom and each molecule atom was measured by a Gaussian function [24].

Internal Validation and Partial Least Squares (PLS) Analysis
Partial least square (PLS) regression analysis was performed on the training set to construct the correlation between the QSAR model and activity values [25].To evaluate the reliability of the models generated from PLS analysis, cross-validation analysis was performed through the leave-one-out (LOO) method, which determines the square of the cross-validation coefficient (q 2 ) and the optimal number of components (ONC).To obtain the non-cross-validation coefficient (r 2 ), a final non-cross-validation analysis was performed using the ONC derived from cross validation analysis and the corresponding standard error of estimate (SEE).The value for q 2 , a measure of the internal quality of the models, was evaluated as follows: ∑ y obs − y pre 2 ∑(y obs − y mean ) 2 where y obs , y pre , and y mean are observed, predicted, and mean activity in the training set, respectively.

External Validation of the QSAR Model
To evaluate the predictive ability of CoMFA and CoMSIA models on the test set, the predictive power of the models generated by the CoMFA or CoMSIA analyses with the training set was assessed by calculating the predictive factor r 2 (r 2 pred ) [26], and measuring the predictive performance of the PLS model.The factor r 2 was calculated as follows: the sum of squared deviation, (SD) between the biological activity of molecules in the test set and the mean biological activity of the training set molecules; the sum of squared deviations between actual and predicted activity values (PRESS), for every molecule in the test set.Coefficients and QSAR results in the contour maps were produced with the field type "STDEV*COEFF".

Molecular Docking
In an effort to explore the interaction mechanism and investigate suitable binding modes, a molecular docking study was performed using the Sybyl package [27].The crystal structure of PPARγ was retrieved from the RCSB (Research Collaboratory for Structural Bioinformatics) Protein Data Bank (PDB ID: 5TWO) [28].In the protein preparation phase, the A-chain was used for the docking study.Crystallized ligands and water molecules of the B-chain were deleted and the hydrogen atoms along with the united atom Gasteiger charges were assigned for the receptor [29].Based on a protomol generation with a threshold parameter of 0.5 and a bloat parameter of 1 Å, the intended active sites where putative ligands could align to and generate potential interactions, were created using the Sybyl package.Binding affinities were presented by Surflex-Dock total scores.In general, conformations of each ligand were ranked by total scores of docking, with the best conformation of the ligand taken into consideration for the corresponding binding interactions.In this study, compound 24c (template) and the newly designed compound N1, were selected and docked to the binding pocket, using the parameters optimized previously.

Conclusions
In this paper, 3D-QSAR and molecular docking studies were utilized to investigate the structural requirements for improving the potency of N-benzylbenzamide derivatives as PPARγ agonists.The established CoMFA and CoMSIA models were both statistically significant, with high external prediction characteristics, indicating that the models could be used to successfully predict compound activity.Surflex-Dock analysis also demonstrated the binding interactions of the template compound with amino acids.Using the model parameter analysis and contour maps, the corresponding structure-activity relationships were determined (Figure 6).Based on the information derived from the different contour maps, several new compounds with improved activities, were designed, further validating the ability of the generated model.We surmise that this will be helpful for the future development of new PPARγ agonists, in the design and screening of new high-activity compounds.
with amino acids.Using the model parameter analysis and contour maps, the corresponding structure-activity relationships were determined (Figure 6).Based on the information derived from the different contour maps, several new compounds with improved activities, were designed, further validating the ability of the generated model.We surmise that this will be helpful for the future development of new PPARγ agonists, in the design and screening of new high-activity compounds.

Figure 1 .
Figure 1.Plots of Actual versus predicted pEC50 values, for the training set and test set compounds, for CoMFA (A) and CoMSIA (B) models.

Figure 1 .
Figure 1.Plots of Actual versus predicted pEC 50 values, for the training set and test set compounds, for CoMFA (A) and CoMSIA (B) models.

Figure 1 .
Figure 1.Plots of Actual versus predicted pEC50 values, for the training set and test set compounds, for CoMFA (A) and CoMSIA (B) models.

Figure 4 .
Figure 4. Docking results.(A) The surface of the binding site, and the conformation comparison of compound 24c (yellow), the original ligand (purple), and the key residues (white) at the binding site; (B) The surface of the binding site and the comparison of the conformation of, compound N1 (yellow), the original ligand (purple), and the key residues (white), at the binding site; (C) The surface of the binding site, and the comparison of the conformation of compound N9 (yellow), the original ligand (purple), and the key residues (white), at the binding site; (D) The surface of the binding site, and the comparison of the conformation of compound N12 (yellow), the original ligand (purple), and the key residues (white), at the binding site; (E) Interaction between compound 24c (yellow) and residues (white); (F) Interaction between compound N1 (yellow) and residues (white); (G) Interaction between compound N9 (yellow) and residues (white); (H) Interaction between compound N12 (yellow) and residues (white).

Figure 4 .
Figure 4. Docking results.(A) The surface of the binding site, and the conformation comparison of compound 24c (yellow), the original ligand (purple), and the key residues (white) at the binding site; (B) The surface of the binding site and the comparison of the conformation of, compound N1 (yellow), the original ligand (purple), and the key residues (white), at the binding site; (C) The surface of the binding site, and the comparison of the conformation of compound N9 (yellow), the original ligand (purple), and the key residues (white), at the binding site; (D) The surface of the binding site, and the comparison of the conformation of compound N12 (yellow), the original ligand (purple), and the key residues (white), at the binding site; (E) Interaction between compound 24c (yellow) and residues (white); (F) Interaction between compound N1 (yellow) and residues (white); (G) Interaction between compound N9 (yellow) and residues (white); (H) Interaction between compound N12 (yellow) and residues (white).

⁎
Test set molecules.

Figure 5 .
Figure 5. Molecular alignment.(A) Common structure retrieved from compound 24; (B) Alignment of the compounds in the training set.

Figure 5 .
Figure 5. Molecular alignment.(A) Common structure retrieved from compound 24; (B) Alignment of the compounds in the training set.

Figure 6 .
Figure 6.Diagram of structure-activity relationship based on core structure of template compound 24c.

Table 1 .
Statistical parameters of the CoMFA and CoMSIA models.

Table 2 .
Newly designed compounds and predictive activity.

Table 2 .
Newly designed compounds and predictive activity.

Table 3 .
Actual and predicted pEC 50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.

Table 3 .
Actual and predicted pEC50 values of PPARγ agonists.
⁎ Test set molecules.