The molecular docking and molecular dynamics study of flavonol synthase and flavonoid 3’-monooxygenase enzymes involved for the enrichment of kaempferol

Abstract Kaempferol is a natural flavonol that shows many pharmacological properties including anti-inflammatory, antioxidant, anticancer, antidiabetic activities etc. It has been reported in many vegetables, fruits, herbs and medicinal plants. The enzyme flavonol synthase (FLS, EC 1.14.20.6) catalyses the conversion of dihydroflavonols to flavonols. Whereas flavonoid 3’-monooxygenase (F3’H, EC 1.14.14.82) catalyses the hydroxylation of dihydroflavonol, and flavonol. FLS is involved in the synthesis of the kaempferol whereas F3’H causes degradation of kaempferol. The present study aimed to analyse the binding affinity, stability and activating activity of enzyme FLS as well as inhibitory activity of enzyme F3’H involved in the enrichment of the kaempferol using the in-silico approaches. Computational study for physico-chemical properties, conserved domain identification, 3-D structure prediction and its validation, conservation analysis, molecular docking followed by molecular dynamics analysis of FLS and F3’H, protein-activator (FLS-LIG Complex) and protein-inhibitor (F3’H-LIG Complex) complexes have been performed. Other structural analyses like root mean square fluctuation (RMSF), root mean square deviation (RMSD), surface area solvent accessibility (SASA), radius of gyration (Rg), hydrogen bond analysis, principal component analysis (PCA), Poisson-Boltzmann analysis (MM_PBSA) and the dynamic cross correlation map (DCCM) analysis to explore the structural, functional and thermodynamic stability of the proteins and the complexes were also studied. The molecular docking result showed that FLS binds strongly with the activator ascorbate (CID _54670067) while F3’H binds with the inhibitor ketoconazole (CID_456201). The most powerful inhibitor (ketoconazole for F3’H) and activator (ascorbate for FLS) is determined by computing the thermodynamic binding free energy through MM_PBSA analysis. The current work provides wide-ranging structural and functional information about FLS and F3’H enzymes showing detailed molecular mechanism of kaempferol biosynthesis and its degradation and hence kaempferol enrichment. Finding of the present work opens up new possibilities for future research towards enrichment of kaempferol by using activator (ascorbate) for FLS and inhibitor (ketoconazole) for F3’H as well as for its large-scale production using in vitro approaches. Communicated by Ramaswamy H. Sarma


Introduction
Medicinal plants have key role in traditional medicine as well as in modern pharmaceutical drugs where different plant parts are used to prepare herbal medicines as well as curing the various diseases. The use of medicinal plants in treating the variety of diseases is a natural protective strategy (Dianat et al., 2014). Millettia pinnata plant also has many medicinal and pharmacological properties such as antimicrobial, antioxidant, anti-inflammatory, anti-cancerous and anti-diabetic (Chopade et al., 2008). Many pharmacologically important phytoconstituents in M. pinnata and a variety of flavones and its derivatives including kaempferol has been isolated (Muqarrabun et al., 2013).
Kaempferol is a characteristic flavonol (a class of flavonoid) and reported from various plants including broccoli, tea, grapes, ginkgo biloba leaves, tomatoes, and medicinal plants . This flavonoid shows numerous pharmacological properties such as cancer prevention agent, antidiabetic, anticancer and antimicrobial (Kim & Park, 2020). The anticancer properties of kaempferol have been accounted for in malignancy cells from various organs including, blood tumors, breast, lungs, ovarian, gastric and pancreatic disease (Sak, 2014;Zhang et al., 2008). Kaempferol interferes with a number of cancerous cells that rely on signalling pathways .
Flavonoid 3 0 -monoxygenase enzyme catalyses the chemical reaction as shown below.
Recent advancement in molecular biology techniques has paved a new platform for researchers to study developmental and inducible regulation of key enzymes for various metabolic pathways in plants. Enzyme inhibitor and enhancer may play significant role in enriching the quantity of a medicinally important plant compounds. Previously, an in-silico study towards inhibition of geraniol dehydrogenase activity was reported using virtual screening of ligand database, homology modelling, and molecular docking (Sahu et al., 2011). Another multistage virtual screening method has been reported for screening of the Cyclin-dependent kinase 2 (CDK2) inhibitors belonging to the family of Ser/Thr protein kinases (Liang et al., 2020). Also, in vitro grown control plantlets could biosynthesize a high quantity of different bioactive cardenolides, which can be enriched using various biotic and abiotic elicitors (Singh et al., 2019). M. pinnata is an important medicinal plant with a wide range of pharmacological properties and have been employed as medicinal agents for a range of diseases (Bhalerao & Sharma, 2014) Many studies have described the beneficial effects of dietary kaempferol in reducing the risk of chronic diseases, especially cancer . Furthermore, due to their interactions with numerous biomolecules, flavanoids are strong antioxidants that have been shown to alter the activity of various enzymes (Dakshayani et al., 2002). Kaempferol, quercetin, karanjin, kanjone, pongaglabrone, gammatin, pongaglabol, kanugin, and other bioflavonoids have been identified in M. pinnata (Satyavati et al., 1987). Kaempferol has more anticancer, antioxidant, and anti-inflammatory properties than other antioxidant present in this plant (Wang et al., 2018). It can be used to treat a variety of acute and chronic inflammation-related disorders, such as colitis and intervertebral disc degeneration, as well as acute lung injury and postmenopausal bone loss (Ren et al., 2019). Although, kaempferol from M. pinnata has several significances but no detailed investigation of this phytoconstituent has been done in this plant. Therefore, use of computational approaches to study the enrichment of kaempferol from M. pinnata is our great of interest.
There are many recent in-silico techniques involving the structure-based discovery of drugs for the ligand-target recognition study from the molecular docking the static one to molecular dynamic enhanced strategies (Salmaso & Moro, 2018). Expanding and underscoring the in-silico approach in terms of the reduction, refinement and replacement, will lead to effective and more pr ecised research in case of cancer (Jean-Quartier et al., 2018). Therefore, it will be more interesting to know whether in-silico-based retrieval and subsequent validation of some elicitors and inhibitors (chemical enhancers such as metal ions and organic compounds) has any effect on enrichment of important components like kaempferol in M. pinnata. Here a detailed study towards the enrichment of kaempferol by using activator of FLS and inhibitor of F3'H towards enhancing and inhibiting the activity of synthesizing and degrading enzyme using computational approach was performed. Many researchers have employed binding free energies to determine ligand-protein interactions (Rahul et al., 2020). This study will provide comprehensive information about metabolic pathways as well as the strategies to use in M. pinnata or other medicinal plants towards enriching its important compounds and use them in traditional medicine. The main goals of the study were to use the molecular docking approach to determine protein-activator, protein-inhibitor interactions and inhibitory potential in the form of binding affinity. Further, molecular dynamics simulations, validation of the molecular docking results and comprehensive analysis of protein-ligand complexes were also studied. Finally, Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) analysis was performed to determine the most powerful inhibitor and activator by computing the thermodynamic binding free energy.

Material and methods
In-silico screening of FLS and F3'H The protein sequences of flavonol synthase (FLS, NP_001237419.1) and flavonoid 3 0 -monooxygenase (F3'H, XP_003525447.1) of Glycine max were retrieved from NCBI. These sequences were chosen depending on the overall quality specifications in FASTA. The conserve domain was identified by NCBI conserve domain identification database. The expasy-ProtParam tool (http://web.expasy.org/protparam) was used to evaluate physiochemical parameters viz. extinction coefficient, theoretical pI, molecular weight, number of amino acids and composition, instability index, aliphatic index, grand average of hydrophobicity (GRAVY) of the FLS and F3'H (Dunker & Obradovic, 2001). According to the recent 7 th CASP experiment rank, I-TASSER is the best approach for the prediction of 3 D structure of the proteins . Therefore, we have used I-TASSER server (Roy et al.,2010) for predicting the 3 D structure of FLS and F3'H proteins. The best model was selected on the based on highest C-score value and the quality of the models was validated by the various tools like Qmean server available at http://swissmodel.expasy.org/qmean (Benkert et al., 2008 and2009), SAVeS Server available at https://services.mbi.ucla.edu/ SAVES and ERRAT (Colovos & Yeates, 1993). For FLS, F3'H proteins, the Consurf algorithm (http://consurfdb.tau.ac.il/) was the tool used for divulging functional regions in macromolecules analysing evolutionary dynamics of atoms in proteins (Celniker et al., 2013).

Docking
Docking is the molecular modelling technique which helps to predict interaction of protein with ligands. As per BRENDA database, the synthesizing enzyme FLS was docked with all the activating compounds that activates the activity of the FLS, while F3'H was docked with all the inhibitors that inhibits the activity of F3'H. The docking study was done by software Maestro-Schrodinger (Bowers et al., 2006) that involves pre-processing of the protein followed by generation of receptor grid of 20 A 0 and then ligand preparation (like 3 D optimization, energy minimization) was done using Ligprep and confgen options available in Maestro and the file was saved as. maegz format. Molecular docking was done for the study of the interaction as well as their binding affinities between protein and ligand using Glide at SP algorithm.

Molecular dynamics simulation
Molecular dynamics (MD) simulation study of proteins FLS and F3'H as well as docked complexes FLS-LIG Complex and F3'H-LIG Complex was performed using GROMACS version 5.1.1(GROningen Machine for Chemical Simulations) (Abraham et al., 2015). The protein parameters were generated using Charm36-Jul2020 forcefield (Vanommeslaeghe et al., 2010), while the ligand parameters were generated using CGenFF server (Yu et al., 2012) for the same forcefield for FLS-LIG Complex and Acpype server (Da Silva & Vranken, 2012) for F3'H-LIG Complex. The gmx solvate tool helps to perform solvation with SPC water model to solvate the specified structures, water model molecules were poured in a cubic simulation box and using the "gmx genion" script and the counter ions were added into the system to electrically stabilise it. GROMACS energy minimization procedures were also followed to stabilise the entire simulation box. After the neutralization process, the energy minimization was performed to remove the steric clashes as well as optimise the structures. Energy minimization was performed that includes 50,000 steps using steepest descent algorithm. The MD balance was carried out in two steps: (i) NVT (No. of molecules, Volume and Temperature) equilibration, 300 K to stabilise the temperature of the system, (ii) NPT (No. of molecules, Pressure and Temperature) equilibration was done to stabilise the density and pressure of the system. LINC (Linear Constraint Solver) algorithm helps to conserve the bond lengths (Hess et al., 1997). Also, the particle mesh Ewald (PME) helps in handling/maintaining the long rang interactions (Kohnke et al., 2020). The SHAKE algorithm was used to limit hydrogen bond lengths and the coulombic and van der Waal interactions were cut to 1.2 nm allowing for a time step of around 2fs. Finally, the MD simulation of this well equilibrated system was performed at desired temperature and pressure for 100 ns. Also, various dynamic analysis parameters like RMSD, RMSF, ROG and others was carried out by gromacs tool (Grant et al., 2006).

Trajectory analysis and protein-ligand interaction energy analysis
The analysis of gromacs trajectory was carried out by using various gromacs tools like the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) using gmx rms and gmx rmsf tools, respectively. Whereas, radius of gyration (Rg) and solvent accessible surface area (SASA) was done by gmx gyrate and gmx sasa, respectively. The gmx do_dssp tool helps in the estimation of secondary structure, while gmx h_bond analysed the hydrogen bonds. The Pymol (DeLano, 2002) and VMD (Humphrey et al., 1996) was used for the visualisation process and the graphs were plotted using OriginPro 9.0 software.

Principle component analysis
The PCA is the generally utilized insightful procedure to display slow and functional movements of the biomolecules (Bahar et al., 1998). PCA was performed on the trajectory to observed the essential motion of the C-alpha atoms present in the protein during the simulation (Rout et al., 2018a;Rout et al., 2018b). The direction of the motion is represented by the eigenvectors whereas magnitude of motion along with direction is represented by eigenvalues (Chen et al., 2017;Yan et al., 2018). The PC components (PC1 and PC2) were obtained by solving and diagonalizing the eigenvectors and eigenvalues for the covariance matrix. PC1, the 1st principal component gives the highly significant information. The porcupine plot was generated using the first eigenvector with the help of pymol and represented in Figure 2.
The PCA analysis was plotted based on following equation: Where, ri is the cartesian coordinate of the ith atom and rj is the cartesian coordinate of the jth atom and t is the average time period over the whole trajectory (Amadei et al., 1993).

MM_PBSA analysis
The Gromacs utilities tool g_mmpbsa was used to calculate the interaction energies of trajectory files. The MM_GBSA approach comprise of three steps calculations to give detail information of the total changes occurring in the binding energy during the simulation. They are (a) Calculation of potential energy (MM) in vaccum that is formed on the parameters of Molecular Mechanics (MM) and it comprises both non-bonded and bonded parameters. The bonded parameters include the torsions, bond, and energies angles, while non-bonded parameters include Electrostatic and Vander waal interactions (Gouda et al., 2020). The Potential energy (E MM ) is represented by the below mentioned equation: And the other two (b) calculation of Polar solvation energy and (c) non-polar solvation energy is account for the calculation of free energy of solvation using an indirect solvation model (Kumari et al., 2014). This is represented by the equation shown below:

DCCM analysis
A prominent method for examining the trajectories of molecular dynamics (MD) simulations is dynamic cross correlation (DCCM) analysis. It is feasible to calculate the dynamic correlation between all atoms within the molecule or the degree to which they move together by analysing the system's trajectories, as after attaining the equilibration every Ca atoms shows the correlative motions (Kormos et al., 2006). And dynamic cross-correlation was performed to check the differences and similarities in the internal dynamics of protein and the complex ). The Dynamic Cross Correlation Map was constructed and analysed using the time evolution of the MD Simulation trajectory ranging from 0 to 100 ns. During the fluctuation/ displacement of the atom, the DCCM plot returns a matrix of all atom-wise cross-correlations that are correlated or anticorrelated with each other. The covariance matrix was built, and the cross-correlation coefficient (M xy) was determined by the motions of two atoms, as shown below: M xy ¼ hDr x: Dr y i=ͱ ½hDr x i2: hDr y i 2 Here, x and y are atoms and Dr x and Dr y are their respective movement vectors.
The N Ã N heatmap is generated by this dynamic cross-correlation tool using Bio3d package (Grant et al., 2006) in "R", where N denotes the no. of C alpha atoms. The values of correlation range from À1 to þ1, where þ1 denote perfect correlation, 0 denote no correlation, and À1 denote complete anti-correlation.

Results and discussion
In-ilico screening of FLS and F3'H The M. pinnata show maximum resembles with Glycine max Wegrzyn et al., 2016). Therefore, G. max protein sequences of flavonol synthase (FLS, NP_001237419.1) and flavonoid 3 0 -monooxygenase (F3'H, XP_003525447.1) containing 334 and 513 amino acids, respectively were selected and retrieved from NCBI for the subsequent study. The Expasy protparam tool was used to find the various chemical and physical properties of the proteins as shown in Table 1.
To quantify the accuracy of the I-TASSER predictions, a scoring function (C-score) based on the relative clustering structure density and the consensus significance score of several threading templates is introduced. With a correlation coefficient of 0.91, a large-scale benchmark test reveals a good correlation between the C-score and the TM-score. Both false positive and false negative rates are less than 0.1 when using a C-score threshold of > À1.5 for models with correct topology. The accuracy of the I-TASSER models could be predicted using C-score and protein length, with an average error of 0.08 for TM-score .
The Ramachandran plot analysis of predicted 3 D model of FLS and F3'H showed 16.1%, 19.4% in additional allowed regions, 81.1% and 76.6% of the amino acid residues in most favoured regions, respectively (Supplementary Figure 1). The models were also validated by ERRAT tool and noticed that FLS and F3'H models have overall quality factor of 86.503 and 94.455, respectively. The Z-score for FLS and F3'H proteins are À4.99 and À5.59, respectively. From the above validation results, it can be predicted that the modelled protein FLS and F3'H are stable, good, compatible, reliable and highly acceptable. The ConSurf algorithm tool predicts the conservation pattern in the protein models FLS and F3'H, the Bordeaux color indicates the highly conserved amino acids, intermediately conserved positions are indicated in white color however, the turquoise color indicates the variable amino acids and yellow color show no conservation score due to the insufficiency of data. In both FLS and F3'H, in one angle variables regions are more than the conserved regions while in other angle conserved regions are more as shown below (Figure 4).

Analysis of protein-ligand interactions
In this investigation, we generated and computed the binding capability of a set of potential activators and inhibitors with the FLS and F3'H, respectively. All of the activators with FLS, while inhibitors with F3'H were docked in their respective binding sites. For the evaluation of molecular docking interactions, the binding affinity, glide score, gscore, glide emodel and glide RMSD to input were the essential parameters. Based on these parameters, we selected activator and inhibitor for FLS and F3'H as listed in Table 2. Importantly, the best binding affinity is shown by activator ascorbate (CID_54670067) with FLS and Inhibitor ketoconazole (CID_456201) with F3'H, having the maximum negative score considering the whole protein as receptor. The molecular docking helps in predicting activity of a molecule bound to a receptor having its best orientation and affinity (Hakes et al., 2007). The 3 D and 2 D interaction of the docked complexes were generated using Discovery studio and are shown in Figure 5.
Present finding infer that the selected activator and inhibitor not only form the interactions found in a standard molecule, but they also formed some significant interactions that could assist the enhancement of the stability of the ligands in the receptor's binding pocket. The key residues involved in H-bonds formation in FLS-LIG Complex are ASP (51), GLU (52), GLY (53) while residues involved in H-bonds includes GLY (447), VAL (300), GLY (303), ALA (302), VAL (117). The receptor-ligand interaction is dynamic therefore the stable binding mode of the ligands could be determined by docking results, but they cannot be regarded unequivocal because of its dynamic nature John et al., 2017;Sharma et al., 2019;Tanwar & Purohit, 2019).
As a result, MD simulations on docked complexes were run for 100 ns, followed by PCA analysis, DCCM and MM-PBSA computations. MD simulations helps to understand the patterns of interaction between protein-ligand complexes as well as their functions Sharma et al., 2020). It helps in assessing the stability of the docked protein-ligand complex, thus validates the docking findings (Rajendran et al., 2018).

Molecular dynamics simulations
Molecular dynamics simulation is a potent tool for studying protein and other biopolymer dynamics. MD simulation is a reliable method for analysing the effect of mutations on protein dynamics, protein folding, target identification and drug development (Rajith & George Priya Doss, 2011;Wu et al., 2008;Xu & Wang, 2006;Bhardwaj et al., 2021, Duan et al., 2010 and many other applications. The docking data are static poses of protein-ligand interactions, thus they can't tell us about the structural changes that occur as a result of ligand binding. The stability of proteins FLS, F3'H and complexes FLS-LIG Complex, F3'H-LIG Complex were monitored by calculating several parameters like RMSD, RMSF, ROG and others for a period of 100 ns of MD Simulation.

Conformational stability of structures
The binding of a ligand to the active site of a protein can alter the overall stability of a protein (Teilum et al., 2011). RMSD is a validated tool for analysing the conformational dynamics and stability of biomolecular structures. The RMSD measurements are sufficient to characterize the ligand's conformational stability in the binding pocket . Significant RMSD fluctuations of all the complexes ranging from 0 to 100 ns are observed as shown in Figure 5. The mean distance between the backbone C-alpha atoms of the superimposed proteins is represented by the RMSD. With an average RMSD of 0.15 nm, the FLS-LIG Complex showed the least deviation from the starting structure. The F3'H and F3'H-LIG Complex complexes deviated slightly higher as compared to FLS-LIG Complex. The average RMSD of F3'H and F3'H-LIG Complex was $ 0.35 nm and 0.4 nm. Beyond 55 ns, the trajectories for all the complexes reach a stable state with average RMSD values ranging 0.28 nm to 0.40 nm. The least rmsd variation is in FLS and FLS-LIG Complex as compared to F3'H and F3'H-LIG Complex.  The minor variances in RMSD values and low levels of fluctuations revealed that the complexes had stable trajectories after simulations and it lays the groundwork for further research. The ligand's stability inside the binding pocket is indicated by the stabilized RMSD trajectories (Rajendran et al., 2016).

Root mean square fluctuation
The RMSF was calculated to determine the reaction of a residue that affects the conformational stability of the protein and the complexes. As shown in Figure 6(B), the binding site RMSF values of proteins FLS and F3'H vary from 0.4 to 1.6 nm, whereas the rmsf of complexes FLS-LIG Complex and F3'H-LIG Complex range from 0.4 to 0.5 nm. The root-mean-square fluctuation (RMSF) of complexes evaluates a conformational stability in the complexes. When a simulation is equilibrated and the structure fluctuates significantly with a stable symmetry, the concept of identifying the fluctuations of each residue of the protein and the complex relative to the simulation usually appears. The RMSF value of residues with remarkable flexibility was higher, as expected in dynamic molecule-protein binding. The RMSF value of all the protein and the complexes demonstrated great adaptability and relativistic vibrational mobility throughout the simulations for every residue marked in the active site, compared to the complex with conventional molecules.

Radius of gyration
The compactness of all four molecules is measured by radius of gyration (Rg) shown in the Figure 6(C). Over a time period of 100 ns at 300 K temperature, all four molecules including proteins FLS, F3'H and the complexes FLS-LIG Complex, F3'H-LIG Complex in the simulation remained externally stable in their compact form. The radius of gyration (Rg) is a measurement of how well regular secondary structures fold into three-dimensional protein structures. Rg denotes a change in the compactness and overall size of the protein structure. The effect of activator and inhibitor binding on the Rg value of the protein was calculated and compared to proteins that did not include any activator or inhibitor. And it is concluded that the Rg value has higher fluctuation in case of F3'H and F3'H-LIG Complex whereas FLS and FLS-LIG Complex shows less fluctuation.

Analysis of inter-molecular hydrogen bonds
Hydrogen bonds play a significant role in protein-ligand interactions (Kollman, 1977). Hydrogen bonds fulfil the specific geometrical requirements for protein-ligand interactions (Yesylevskyy et al., 2006), and hence it requires thorough analysis. As shown in Figure 6(E), we computed the number of hydrogen bonds between the selected protein-ligand complexes. The standard activator (ascorbate) formed less hydrogen bonds as compared to standard inhibitor (ketoconazole) during the simulation in the protein-ligand complex and was more strongly linked to the active site than the F3'H-LIG Complex ( Figure 6E). Hydrogen bonds are responsible for the stiffness and stability of proteins and complexes.
The biomolecular surface which is accessible for the solvent molecules are determined by the parameter called Solvent Accessibility Surface Area (SASA). The SASA value of FLS-LIG Complex and F3'H-LIG Complex is higher as shown in shown in Figure 6(D). However, it helps in the determination of the changes in protein before and after docking as FLS and F3'H SASA value increases after docking.

Principle component analysis
The Principal Component Analysis helps in the determination of overall movement of c-alpha atoms of both the structures proteins and complex during the simulation and it also concludes the behaviour of conformational changes with gained principal components i.e., eigenvectors (which details the direction of maximum variance of the structures) and eigenvalue (describes the total percentage of variance of positional fluctuations of atoms present in each dimension). There is periodic jump of conformers throughout the simulation and this is indicated by -blue to -white to -red, color scale. The maximum information is stored in PC1 (Mart ınez-Archundia et al., 2019) and contribute the maximum percentage of variance in the structure i.e. (60-80%). The protein and complex projection are shown as scattered plot PC1 v/s PC2 (Figure 7). The PCA utilises the results of MD simulation trajectory for demonstrating the motion of the c-alpha atoms present in the protein and the complex. The Porcupine plot (Figure 2) shows the overall movement of the Ca atoms present in the protein and the complex. This movement gives on to the opening of the active sites and thus facilitates the ligand binding and then the other amino acid moves from the loop side chain and thus closes the active sites.
To evaluate the link between different conformations in the samples during the journey, PCA was used on each protein as well as complexes. The obtained principal components (orthogonal eigenvectors) describing conformational differences can be used to determine the nature of conformational differences the structure's axes of maximum variance, the positional fluctuations captured in each dimension is the eigenvalue of that dimension. The colour scale, which runs from blue to red, represents the periodic jump between conformers along the trajectory. Our result demonstrates that the first few PCs account for about 60-80% of structural variance FLS ( Figure  7A). For all structures, the PC1 and PC2 describe the information concerning trajectory jumping of phases. However, in Figure 7(C) and (D), no difference in conformations is observed before and after docking of F3'H. The Protein and its complexes were projected along PC1 vs PC2, and the results were displayed in a scatter plot. Four separate groupings of points along PC1 and PC2 could be depicted in the dispersed plot (Grant et al., 2006).  The loop region showed the largest variation, indicating structural and conformational changes taking place in this area to accommodate substrates of different sizes and shapes. Ensemble based PCA investigates three PCs that account for roughly 90% of structural variance and are divided into four groups, each of which has one subgroup. The maximum fluctuations are mainly due to loop region indicating that the conformational and structural changes occur in this region to accommodate substrates of various shapes and sizes.

MM-PBSA analysis
Using g_mmpbsa gromacs tool (that works on the principle of Poisson-Boltzmann) we have calculated free binding energy and electrostatic interaction of the complexes to see the detailed distribution of electric potential in the solution and also to interpret the binding mechanisms namely van der Waals interaction, electrostatic interactions, polar and non-polar solvation energy. Scripts of selected complexes were extracted to probe the computational study. The summary of MM-PBSA analysis of FLS-LIG Complex and F3'H-LIG Complex as listed in the Table 3. The binding free energy study revealed that the F3'H-LIG Complex (-213.104 þ/-14.909) has the highest DE binding energy (kJ/mol) as compared to FLS-LIG COMPLEX (-177.768 þ/-0.917). Our result (Table 3) shows the binding free energy of the FLS-LIG Complex, F3'H-LIG Complex separated into polar solvation, SASA, van der Waal and electrostatic interactions. Electrostatic energy is the primary controller and supporter of protein-ligand interactions while Van der Waal energy has a minor role. In the activation of the loop, negative electrostatic interactions are thought to maintain the stability of the loop and hence enhance the affinity. The SASA played an advantageous role in complex formation whereas polar solvation energy has more unfavourable role (positive value). The stability of complexes is achieved with a reasonable binding free energy as per the interpretation of the protein trajectory. F3'H-LIG Complex has higher binding energy than FLS-LIG Complex inferring that it binds to the receptor more tightly. Furthermore, these findings laid the groundwork for further comprehensive research into the development of novel activators and inhibitors for FLS and F3'H, respectively.
The summary of MM-PBSA analysis of FLS-activator and F3'H-inhibitor is listed in the table below.

Dynamic cross correlation map analysis
During atomic fluctuation/displacement a DCCM plot is associated with a matrix of all-atom wise cross-correlation (i. e. either correlated or anticorrelated). Each built-up structure as well as the protein-ligand complex structure was subjected to DCCM analysis (Figure 8). Our result shows the stability of the proteins and the complexes. The correlations range from À1 to þ1, or negative to positive, and are depicted in pink and cyan, respectively. FLS-LIG Complex has maximum residues in positive correlation, while F3'H residues has a slightly weaker but still positive correlation (Figure 8). FLS has the most residues in negative correlation while F3'H-LIG Complex has the most residues in negative correlation but also has some positive correlation (Figure 8). The overall correlation between the binding site and the other amino acids has increased as observed in the DCCM plot of FLS-LIG Complex ( Figure 8). In FLS-LIG Complex and F3'H-LIG Complex, the correlation between binding areas increases. The DCCM plot validates all the structures displaying the correlation map of each residue. From result it is concluded that F3'H-LIG Complex shows a strong open structure ensemble in the protein that relates to the flexibility of different residues in this state. However, FLS-LIG complex shows both open and closed conformations.

Conclusion
The present work includes the extensive bioinformatics analysis of the enzyme FLS and F3'H involved in the biosynthesis and degradation of kaempferol. The conserved domain of FLS belongs to PLN02704 super family (cl31547), while conserved domain of F3'H resembles to p450 super family (cl12078). The 3 D structure prediction revealed good, reliable, stable, compatible and highly acceptable structure as validated by Ramachandran plot, saves server, Errat, verify 3 D, Consurf algorithm tool. Molecular dynamics simulation of 100 ns of proteins (FLS and F3'H) as well as protein-ligand complexes exhibited the thermodynamics stability. Molecular docking and molecular dynamic simulation study of proteins and complexes (FLS-LIG Complex, F3'H-LIG Complex) analyses (like RMSD, Rg, RMSF, SASA, Hydrogen bond analysis, PCA, DCCM and MM_PBSA) has been done to check the binding affinity, functional and structural stability. The Vdw interactions were responsible for binding free energy and helps in defining the binding affinity of ligands. The F3'H-LIG Complex has more binding free energy than FLS-LIG Complex. The newly selected activator (ascorbate for FLS) and inhibitor (ketoconazole for F3'H) was found to be best towards a strong inhibitory capability of Ketoconazole for F3'H and strong activator capability of ascorbate for FLS. Also, this is the first and a comprehensive in-silico analysis of the flavonol synthase and flavonoid 3 0 -monooxygenase enzymes in which we have demonstrated the sequence analysis to MD simulation using different computational approaches. As these enzymes, flavonoids and medicinally important components are widely distributed in many monocots, dicot and medicinal plants, the experiment can be repeated for their enrichment analysis in other plants. Moreover, validation of in silico findings need to be carried out in future using in vitro based experiment.