A Mathematical Modeling and Molecular Dynamic Simulations in the Investigation of Novel Type I Diabetes Treatment

In the present analysis, 30 compounds were selected as novel type 1 diabetes inhibitors to combat diabetes mellitus resistance and were subjected to 2D-QSAR and 3D-QSAR studies using CoMFA (ffdsel) and CoMFA (uvepls). In the 2D-QSAR study, multiple linear regression was used. A decent relationship was found between the experimental activity and those got by model 1 and model 2 such as (r2 = 0.7616 and q2 = 0.6327) and (r2 = 0.7721 and q2 = 0.5834), respectively. Model 1 with a low standard error of estimate and high valve of q2 seen promising than model 2. In the 3D-QSAR study, docking-based alignment was used. The results showed that CoMFA (uvepls) (q2 = 0.6897; r2 = 0.9999) have good stability and predictability. The internal validation indicated that CoMFA (uvepls) MIFs possess good predictive power than COMFA (ffdsel). The molecular docking study showed that conventional hydrogen bonding, hydrophobic (pi-Alkyl) and van der Waals interactions are key residues such as Glu50, Ser271, Arg 266, Pro270 and Tyr66 play an important role in ligand- receptor binding. MD simulation analysis showed that compound 17 formed hydrophobic (pi-Alkyl) and conventional hydrogen bonds with the residues (Glu50, Ser271, Arg 266, Pro270 and Tyr66), and binds more closely to type 1 diabetes. These results provide useful insights into the discovery and design of a new generation of type 1 diabetes inhibitors.


Introduction
Diabetes, otherwise known as blood sugar is a stipulation that vitiates the body's power to process blood glucose [1]. There are three (3) main types of diabetes: type 1, type 2 diabetes and gestational diabetes. Other less prevalent are monogenic diabetes and cystic fibrosis-related diabetes. Type 1 diabetes, also known as juvenile diabetes, occurs because the body does not produce insulin [2]. Individuals with type I diabetes are insulin-subordinate, which means that they can take false insulin every day to stay alive has been one of the serious issues jeopardizing human life since it was accounted for in 1988 [3,4]. World health organization (WHO) anticipated that 347 million individuals have blood sugar (diabetes), which is anticipated to be the seventh driving reason of deaths by as found out by Boukarai, et al. [7], utilized use of principal components analysis (PCA), multiple linear regression (MLR), multiple non-linear regression (MNLR) and the artificial neural network (ANN) to plan 2D QSAR model of Flavonoid subordinates as in vitro inhibitors operators of Aldose Reductase (ALR2) protein for diabetic difficulties. They revealed good internal and external predictions. The measurable outcomes demonstrate that the model is factually critical and shows excellent soundness towards information variety in leave-one-out (LOO) cross-validation.
Vitthal, et al. [19] who estimate the 3D QSAR and pharmacophore modeling on substituted cyanopyrrolidines as future dipeptidyl peptidase-IV inhibitors of type II antidiabetic agents. They reported that the inhibition provides useful information on the structural characteristics that provide the cyanopyrrolidine inhibitory potency. Pharmacophore modeling was used in addition to their work to identify the structural features of cyanopyrrolidines which are essential to the structure's biological activity. In this work, we perform 2D-QSAR, 3D-QSAR (CoMFA), docking, and molecular dynamics (MD) modeling study of 30 compounds as LYP inhibitors (type 1 diabetes). The model can be successfully employed to design selective diabetes type 1 inhibitors and to discover novel frameworks through screening of compound databases.

Methodology
The dataset used in this study was taken from the PubChem [https://pubchem.ncbi.nlm.nih.gov/bioassay/435024].The studied materials of this research work are 30 compounds. These compounds' structure and their PubChem-CID, experimental, and predicted binding values (-log(IC50 * 10 ^ -6) are available in Table 1. Before anything else, the structures of all the compounds have been optimized and then encrypted the structural features (descriptors) of the molecules. We used Spantan's version 14 program with a semi-empirical PM3 basis set for geometrical optimization. Finally, chemical descriptors were revealing for each molecule that was obtained from the PaDEL-Descriptor [20] output. Molecular descriptors calculated with PaDEL-Descriptor [20] were subsequently subjected to variable reduction using a Data Pre-Treatment Tool (V-WSP Tool) [21]. The data set was then separated into two groups: training and test sets using the [http:/dtclab.webs. com/software-tools] technique of Kennard Stone. For the model generation, the training set, consisting of 18 molecules, was used and the test set, consisting of 6 molecules, was used to take care of the overtraining and evaluate the predictive power using multiple linear regression (MLR) of the produced model.  In the development of CoMFA (FFDSEL and UVEPLS) studies, alignment of the chemical structure is of crucial importance. The precision of the CoMFA model forecasts and the stability of the contour models are highly dependent on the structural orientation of the molecules [22]. As a prototype for superimposition, the most active stable compound was used, assuming its conformation is the most bioactive conformation at the level of the active receptor site.

Molecular Alignment
In Open3DALIGN [23], the fragment was chosen from the docking-base alignment method and the majority of the molecules had been matched with it. In Figure 1, the compatible compounds are shown. The 3D descriptors were used as independent variables and the pIC50 as the dependent variable to derive 3D-QSAR models. Open3DQSAR, a free software, is used today for highthroughput MIFs analysis for pharmacophore exploration [24].
This Open3DQSAR method is useful for the pharmacophore evaluation and drug formulation based on ligand [24]. Recent studies show that the findings obtained using this technique are comparable to CoMFA / CoMSIA [22]. Molecular interaction fields (MIFs) for van-der-Waals and electrostatic field interactions for Open3DQSAR analysis were measured using a grid box around the aligned structures with a phase size of 1.5Å and a distance of 5.0Å [22], similar to the CoMFA report.

Drugs
The 3D Structure of Lyp with inhibitor I-C11 receptor protein (PDB code: 2QCT) crystal structure [25]

Molecular Dynamic Simulations
Molecular dynamics calculations are performed using VMD [27] and NAMD [28] and the CHARMM force field [29]. The interaction parameters were computed using the CHARMM27 force field.
The periodic boundary condition was used and the system was immersed in a cubic water box of extended simple point charge water molecules to calculate the solvated protein. The protein was solvated with explicit water with 0.1M NaCl salt concentration for neutralization. To maximize the initial structure of the protein, the minimization was performed. After that, the system temperature was steadily heated up by 100ps from 0 K to 310 K. Finally, at 310 K for 100ps with the NVT ensemble, the system was balanced. The system was simulated for 1ns (50,000 steps) of chemical time. The MD simulation and results from the analysis were performed on the DELL INSPIRON; Pentium® Dual-Core CPU T4500 @ 2.30GHz and 3GB of RAM, 64-bit-Operating System, x64-based processor.
A full summary of the input parameters is given in the NAMD documentation (www.ks.uiuc.edu/Research/namd/).

Results and Discussion
To discover the structure-activity relationship (SAR) of 24 different compounds acting as antidiabetics, a 2D-QSAR study was     The observed and predicted pIC50 is shown in Figure 2, a list of the descriptors that led to the current analysis and the predicted behavior is shown in Table 2. By measuring their variance inflation factors (VIF), which can be interpreted mathematically as   [32,33].
The mean effect indicates the contribution of each descriptor to the model and its values are shown in Table 3 for Model 1 and Table 4 for Model 2 respectively. Its sign shows the way of variance in task values as a function of the rise or decrease in descriptor values. As seen in Tables 3 & 4, ATS3e and bpol descriptors have high mean effect values relative to other descriptors, which indicates that these descriptors have a significant impact on the pIC50 of the compounds studied.

Results of the 3D-QSAR Models
Interpretation of 3D-QSAR Polyhedrons: The plot of actual and predicted activity obtained from CoMFA (ffdsel) and CoMFA (uvepls) analyses are shown in Figure 4. The best statistical analysis obtained by the CoMFA (ffdsel) and CoMFA (uvepls) are listed in Table   7, which indicated that the obtained 3D-QSAR model is reliable and able to predict binding affinities of the new compound. In CoMFA (ffdsel), the van der Waals (steric) and Coulomb (electrostatic) contributions were found to be 52% and 48%, respectively. In CoMFA (uvepls), the van der Waals and Coulomb contributions were found to be 54% and 46%, respectively. Thus, the steric field had a larger effect than the electrostatic field, indicating that the steric interactions of the molecules with the receptor could be an important factor for anti-diabetes activity. The van der Waals and Coulomb contributions were also evaluated by LOO, LTO, and LMO cross-validation and test set evaluation. These models showed good correlative and predictive values which indicated that the CoMFA (ffdsel) models were less in need of the fields engaged. The statistical importance of the CoMFA (uvepls) 3D-QSAR models is far greater than that of the QSAR models developed using the 2D method, as shown by stronger q2 values and smaller standard error estimates [37]. Furthermore, the QSAR models developed using 3D methods offer an opportunity to describe the observed behavior variance in terms of contour maps as discussed below. The contour maps created during docking-based Open3DQSAR analyses are depicted in Figure 5-8 respectively.     where an increase or decrease of steric bulk, respectively, are predicted to enhance binding affinity. In the CoMFA (ffdsel) steric contour map of compound 4, this steric bulk (red color) was found making interaction with the compound base on the auxochrome that is attached to cyclohexane and benzene rings ( Figure 5A).
The steric contour map of compound 7, the large red coloration scattering around the compound indicating increased steric bulk located around compound 7 ( Figure 6A) will enhance the activity.
These contours are represented in Figure 5B.  Figures 8E-8H). This means that the substituents around these regions should be bulky, that is why compound 4 has the highest binding affinity than compound 7, another green zone is found close around compound 4, indicating that electron-donating character is favored and would exhibit high diabetic inhibitory activity. An increase in the electron density/ electropositive (green contour maps) in this region as shown in Figure 8G and decrease in electron density (yellow contour maps) in this region will increase the potency of compound 7.

Interpretation of Molecular Docking Results
The derivatives were docked to the receptor (PDB: 2QCT) and the energy values were computed using PyRx (Autodock vina residues as shown in Figure 9. The inhibitory assay and molecular docking study reveal that compound 4 inhibits insulin epitope drives type 1 diabetes may be because of its hydrogen bond interactions with Gln63, Tyr7, and Tyr171 amino acid residues (Table 8).

Interpretation of MDs Simulation
By examining the results of 2D-QSAR, 3D-QSAR, and molecular docking, compound 4 was utilized as the lead compound, and the predominant changed area dictated by the above outcomes was almost altered. The root means square deviation (RMSD) plot, temperature, kinetic energy, and potential energy plot of protein during the MD simulations are shown in Figure 10. The simulation result reveals that the RMSD tends to be steady and wavered approximately at 2.0 Å, the kinetic and potential energy is fluctuated approximately at 28,896.58 and -104,501 kcal/mol, respectively, after running 50000 steps (1ns), which implies that the structure of the complex is basically in a stable state. Molecular docking analysis of compound 4 was performed on the A-chain of simulated protein (PDB ID: 2QCT) as shown in Figure 10 using the 2016 discovery studio client, this did not result in an improvement in the number of hydrogen bond interactions ( Figure 11). This was fundamentally similar to the underlying docking systems, thereby showing the efficiency of the docking. It also proposed that compound 4 formed a stable protein complex when protein was simulated to or beyond

Conclusion
In this study, a molecular modeling study was performed on 30 compounds as type 1 diabetes inhibitors using 2D-QSAR, 3D-QSAR, molecular docking, and MD simulations. The 2D-QSAR, CoMFA, and CoMSIA approaches were used to conduct the QSAR analysis and the final models were adequate based on the findings of the statistical validation. The contour maps obtained described the structural activity relationship (SAR) of the studied compounds and were capable of predicting the activity of the compound training set reliably ( Figure 10D). The molecular docking analysis was then used to describe the effects of 3D-QSAR and to investigate the interactions between compound 4 and ligands; Tyr7, Gln63, and Tyr171 were identified as essential residues for inhibitory activity.
However, under realistic conditions, the findings of molecular docking are refined by MD simulations to show the consistency of the association between protein and ligand. The proteinligand interactions described in the molecular docking phase are almost stable under dynamic conditions. We thus had a greater understanding of the SAR and the binding properties of compound 4 as a type 1 diabetes regulator using molecular simulation. We assume that these findings may be useful in the production of new type 1 diabetes inhibitors.