An in-silico receptor-pharmacophore based multistep molecular docking and simulation study to evaluate the inhibitory potentials against NS1 of DENV-2

Abstract DENV-2 strain is the most fatal and infectious of the five dengue virus serotypes. The non-structural protein NS1 encoded by its genome is the most significant protein required for viral pathogenesis and replication inside the host body. Thus, targeting the NS1 protein and designing an inhibitor to limit its stability and secretion is a propitious attempt in our fight against dengue. Four novel inhibitors are designed to target the conserved cysteine residues (C55, C313, C316, and C329) and glycosylation sites (N130 and N207) of the NS1 protein in an attempt to halt the spread of the dengue infection in the host body altogether. Numerous computer-aided drug designing techniques including molecular docking, molecular dynamics simulation, virtual screening, principal component analysis, and dynamic cross-correlation matrix were employed to determine the structural and functional activity of the NS1-inhibitor complexes. From our analysis, it was evident that the extent of structural and atomic level fluctuations of the ligand-bound protein exhibited a declining trend in contrast to unbound protein which was prominently noticeable through the RMSD, RMSF, Rg, and SASA graphs. The ADMET analysis of the four ligands revealed a promising pharmacokinetics and pharmacodynamic profile, along with good bioavailability and toxicity properties. The proposed drugs when bound to the targeted cavities resulted in stable conformations in comparison to their unbound state, implying they have good affinity promising effective drug action. Thus, they can be tested in vitro and used as potential anti-dengue drugs. Communicated by Ramaswamy H. Sarma


Introduction
The incidence of dengue and the proliferation of its viral mosquito vectors (most commonly Aedes aegypti and Aedes albopictus), has significantly increased in recent years exposing the fatality of the disease.Dengue as 'a spectrum of disease' is manifested as dengue fever (DF), dengue hemorrhagic fever (DHF), and dengue shock syndrome (DSS).Dengue is caused by an arthropod-borne virus called dengue virus (DENV), belonging to the family Flaviridae and genus Flavivirus (Murugesan & Manoharan, 2020).World Health Organization (WHO) has conferred the disease an 'endemic' status in recent times, spreading over 100 countries, affecting up to 400 million people annually worldwide causing around 40,000 fatalities out of which 25,000 cases were classified as severe.A high number of cases were reported in Bangladesh (101,000), Malaysia (131,000) Philippines (420,000), and Vietnam (320,000) in Asia during the 2019 outbreak.The endemic regions of dengue as per the WHO report are the US, South-East Asia, and Western Pacific, with Asia, carrying about 70% of the global disease burden.The statistics suggest dengue is now even spreading to new areas including Europe, typically France and Croatia.In 2021-2022, a total of 164,103 dengue cases have been reported in India alone (https://economictimes.indiatimes.com/news/india/india-reported-1-64-lakh-dengue-cases-during-2021-against-2-05-lakh-cases-in-2019-govt-to-rajya-sabha/articleshow/88009894.cms).Among the different climatic zones, five dengueendemic states (Punjab, Haryana, Rajasthan, Gujarat, and Kerala) have been the hotspots of dengue in the country, contributing to more than 50% of the country's annual infections (Mutheneni et al., 2017).Fever, rash, headache, eye discomfort, nausea, and vomiting are frequent clinical manifestations of dengue fever.Patients with severe dengue suffer from plasma leakage, fluid buildup, respiratory difficulties, severe bleeding, or organ dysfunction, and yet there has been no successful discovery of an antiviral drug in the market to target dengue, till date (Srikiatkhachorn, 2009).Dengvaxia, DENVax, and TV003 are few dengue vaccines available (Torres-Flores et al., 2022).But these vaccines show major side effects, such as headache, malaise, low-grade rash, and transient local reactogenicity (Torres-Flores et al., 2022).Implementing the vaccine in the field will be difficult due to vaccine cost and stability issues, long-term monitoring for potential adverse events, and compatibility with current vaccination schedules (McArthur et al., 2013).An anti-malarial drug called chloroquine has shown some of its effects as an antiviral drug.However, its application found limitations when chloroquine was unable to decrease the overall duration for which the disease lasted (Borges et al., 2013).
There are four genetically distinct dengue virus serotypes (DENV-1, -2, -3, and -4) with a fifth strain namely, DENV-5, isolated back in October 2015, having a milder impact but a higher secondary infection rate (Mustafa et al., 2015).DENV-2 strain of dengue is observed to be linked with a drastic decrease in platelet count, which results in fatal consequences in the infected individuals, making it the most life-threatening strain.The genome of the dengue virus type 2 (DENV-2) is a positive sense single-stranded RNA and is �11 kb in length, processed inside the host cell.It encodes three structural proteins (capsid [C], membrane [M], and envelope [E]) that make up the virion's components, and seven non-structural proteins (NS1, NS2 A/B, NS3, NS4 A/B, and NS5) (Perera & Kuhn, 2008).Among all the NS proteins, NS1 is the solitary protein being capable of being both, anchored on the surface of infected cells in the form of a membrane-associated homodimer and being diverged from infected cells into the blood stream as a hexamer (Rastogi et al., 2016).The DENV NS1 glycoprotein is a highly conserved 46 kD glycoprotein consisting of 12 conserved cysteine residues and two N-linked glycosylation sites (Flamand et al., 1999).A unique feature of this protein is to repeatedly polymerize and exist in various monomeric, dimeric, and barrel-shaped hexameric forms (Akey et al., 2015).The polymerization and the distribution of NS1 which are important for its functionality, depend on its glycosylation site.This polymerization ensues in the endoplasmic reticulum (ER) lumen.NS1 is first expressed in the infected cells as a monomer.It translocates co-transcriptionally to the ER lumen, where it undergoes modification and produces homodimers that adhere to the cell membrane and organelle membranes (Alcal� a et al., 2018).Presumably, dimerization may be the cause of the change of NS1 into a membrane-associated protein with enhanced hydrophobicity.NS1 dimers effectively localize with cholesterol-associated molecules on membranes of cells and also contribute to the modulation of cellular lipid dynamics during flavivirus infection (Mackenzie et al., 1996).The third and most recurrent form of NS1 is the hexameric structure formed with glycosylated combinations with lipids.Since hexameric NS1 and high-density lipoprotein are similar, it has been suggested that NS1 might interfere with endogenous lipoprotein particle interaction or biogenesis and hence alter the coagulation cascade (Chen et al., 2018).The progression of dimer to hexamer formation of NS1 protein has been experimentally demonstrated by scientists like Gutsche et al. (2011).
Considering an experimental demonstration by Fan et al. (2014) using DENV transfected cells it has been derived that NS1 mutants, with alanine replacement C55A was incapable of causing CPE and generating any visible plaques in the transfected cells (Płaszczyca et al., 2019).This inability might trace back to the fact that mutation at C55 rendered NS1 dysfunctional in giving rise to any viral particle (Fan et al., 2014;Płaszczyca et al., 2019).In another research conducted by Pryor and Wright (Pryor & Wright, 1993) it was observed that alanine scanning at C313, C316, C329 produced only monomeric versions of the proteins, indicating towards their inability to dimerize or form dimers of reduced stability.Thus, it can be hypothesized that mutation of cysteine residues impairs the correct formation of disulphide bonds within the system of NS1 in its entirety.Notably, mutation C55A in NS1 hampers its functionality in inducing viral replication.Evidently, the cysteine residues at the carboxyl end of NS1 have elucidated to play a crucial role in functionality of the protein mediated by dimer formation.
The maturation of the NS1 protein is primarily dependent on glycosylation, which further facilitates the virus's secretion, pathogenesis, and replication (Flamand et al., 1999).First, the NS1 monomer is glycosylated with high mannose glycans at the residues namely, N130 and N207.Again, the result of an experiment conducted by Yap et al. (2017) established that reduced yields of infectious DENV strains were obtained when N130 and N207 conserved glycan sites were deglycosylated.The viral protein's oligomer is rendered unstable due to disruption in its folding, and solubility patterns by the lack of glycosylation (Yap et al., 2017).Inference derived from a mutational study administered by Somnuke et al. (2011) established that N130Q mutation induced monomer formation.Inspecting and analyzing all the previous experimental outcomes, it can be proposed that N130 stimulates the oligomeric assembly of secreted NS1 by affecting its generation efficiency along with the stabilization of the secreted hexamer (Somnuke et al., 2011).On the other hand, glycosylation of residue N207 is involved in hyperpermeability and endothelial glycocalyx-like layer (EGL) degradation.A detailed analysis demonstrates the N-linked glycan at residue 207 facilitates secretion and extracellular protein stability, responsible for NS1's affinity for binding to the complement host cells and is crucial for NS1 endosomal transport during clathrin-mediated endocytosis, a typical mechanism in the pathogenesis of NS1 (Dahlb€ ack et al., 1983).Through interactions with the plasma complement regulator C4BP, NS1 binds to the cell surface through N207.Evidently, NS1 roleplays as a protectant of the virus against complement-mediated neutralization and averts the death of infected cells (Avirutnan et al., 2011;Leung et al., 2006;Rawal et al., 2009).Hence, N-inked glycan at position 130 is experimentally observed to be vital for the stability of the hexamer of NS1 protein whereas the N-inked glycan at 207 th residue is crucial for secretory purposes and extracellular protein stability.
Our research approach aims to design four de-novo inhibitors targeting four specific sites: the first inhibitor (D1) is to block the terminal cysteine residues (C313, C316, and C329) that mediate homodimer formation, subsequent polymerization, and structure sustenance.The second inhibitor (D2) would be designed targeting and blocking the N130 site to disrupt hexametric stability of the NS1 protein as it is one of the pivotal glycosylation sites, vital to their stability, contagiousness, and antigenicity.The third inhibitor (D3) is meant to block the C55 residue position that is engaged in forming disulphide linkages to strengthen the protein monomeric folds, thus in turn stopping the CPE.Lastly, the fourth inhibitor (D4) would be targeting and blocking another active glycosylation site at N207 that is concerned with extracellular protein stability and cell surface attachment and thus is hypothesized to halt the infection in its primary stage (Figure 1).

Structure extraction of NS1 protein
The NS1 protein structure having entry ID 4O6B was retrieved from RCSB PDB.Discovery Studio Visualizer (DSV) was employed for the conversion of dimeric form to its monomeric form by eliminating the recurring sequence of the protein.

Cavity detection and pharmacophoric features extraction using CavityPlus
CavityPlus (Xu et al., 2018) facilitates the computational detection of protein cavities for successful drug designing, and discovery endeavor including target selection and identification.CAVITY module identifies the binding sites on the surface of the designated protein and classifies them based on their druggability and ligandibility.The CavPharmer module initiates pharmocophore feature modelling within the cavities.CavPharmer generates six pharmacophore class features which are hydrogen bond donor, hydrogen bond acceptor, positive ionizable center, negative ionizable center, hydrophobic interactions, and aromatic rings (Xu et al., 2018).From the CavityPlus result obtained, four favorable cavities were selected (D1_CAV, D2_CAV, D3_CAV, and D4_CAV orderly) based on the occurrence of active residues (C313-C316-C329, N130, C55, and N207) which were then chosen and submitted to CavPharmer for generation of pharmacophoric features of the cavities.

Pharmacophore-based database screening and construction of initial ligand libraries
This step includes the screening of chemical libraries for the recognition of compounds that exhibit specific pharmacophoric features and as a result expresses similar biological activity towards the target molecule.ZINCPharmer identifies ligands with similar pharmacophoric features, and therapeutic and interaction profiles indispensable for the activity, affinity, interaction, and predicted inhibitory action of the ligands on the target domain regions (Koes & Camacho, 2012).These features could be enlisted as hydrogen bond acceptor, hydrogen bond donor, hydrophobic center, and negative or positive ion.Various essential filters were applied for the higher precision of the search engine.The RMSD values for D1_CAV, D2_CAV, D3_CAV, and D4_CAV were set at 0.5, 0.31, 1, and 0.9 Å, respectively.For all the aforementioned cavities, the molecular weight of the drug and druglike small molecules was set between 180 and 480 and the number of rotatable bonds was 0-9 as per Lipinski's rule.Consequently, the relative hits of all the cavities were generated.ZINCPharmer generated a library of 43,298; 38,028; 40,565; and 51,916 compounds against cavities D1_ CAV, D2_CAV, D3_CAV, and D4_CAV, respectively.The libraries were saved in sdf format for further study.This ligand library was subjected to redundancy removal and energy minimization using Open Babel.For energy minimization, MMFF94 forcefield was utilized.Steepest descent and successively conjugate gradient algorithm with 2500 default steps were performed.

Construction of final ligand library
The ZINCPharmer generated a library of 43,298; 38,028; 40,565; and 51,916 compounds against cavities D1_CAV, D2_ CAV, D3_CAV, and D4_CAV.After redundancy removal and energy minimization, a library of 13,000; 20,401; 13,257; and 11,520 was generated against cavities D1_CAV, D2_CAV, D3_ CAV, and D4_CAV, respectively.This final library will be used for docking onto the target receptor NS1.

Structure-based virtual screening
AutoDock Vina (Eberhardt et al., 2021) was utilized to conduct structure-based virtual screening.Vina scoring function utilizes a Monte-Carlo (MC) iterated search combined with the BFGS gradient-based optimizer.Firstly, a configuration file was created for each of the four cavities by recognizing their respective dimensions with respect to the 3D space.The grid box was adjusted manually for the active cavities (D1_CAV, D2_CAV, D3_CAV, and D4_CAV).CAV7 was set with the coordinates X: 32.307, Y: 53.415, Z: 45.081, and dimensions (Å) X: 72, Y: 40, Z: 62.For CAV10 coordinates were set at X: 22.873, Y: 46.21, Z: 5.865 and dimensions (Å) X: 40, Y: 46, Z: 42.CAV14 was set at coordinates X: 30.987,Y: 47.455, Z: 6.434 and dimensions (Å) X: 40, Y: 40, Z: 40.For CAV15 coordinates were set at X: 38.005, Y: 71.822, Z: 32.38 and dimensions (Å) X: 40, Y: 34, Z: 40.For all the four cavities, exhaustiveness was set at 8, maximum number of binding modes generation was set to 1 and energy range (kcal/mol) was set to 4. Consequently, Instadock was used to screen the generated output library to obtain the top 100 ligands from each ligand library in accordance with the best binding affinity to their respective cavities and was saved in csv format.

ADMET evaluation
ADME evaluation is the most pivotal step in drug discovery and designing.It optimizes the balance of different drug properties, i.e. absorption, distribution, metabolism, and excretion which are necessary to convert the leads into effective medicines and reduce the toxicity level and increase the druggability (Wang et al., 2015).For this evaluation, SwissADME (Daina et al., 2017) was used.The top 100 ligands were provided as input an only the leads generated by the server which completely fell under the bio sensitivity radar plot, were selected for further analysis.The top 50 leads were selected for further progression.For the (3) NS1 protein with the target residues highlighted in purple (C55), yellow (N130), pink (N207) and green color (C313, C316, C329).( 4), ( 5), ( 7) monomeric (4), dimeric (5), and hexameric forms (7).NS1 is first expressed in infected cells as a monomer.It translocates co-transcriptionally to the ER lumen.( 6) Further oligomerization leads to the formation of NS1 hexameric structure formed with glycosylated combinations with lipids.(4a, 4b) D3, blocking the C55 residue position that is involved in the formation of disulphide linkages, thus in turn stopping the CPE.(5a) D1 will block C313, C316, and C329 that is majorly involved in homodimer formation.The dimers further modify cellular lipid dynamics during CPE and/or lead to the recruitment of cholesterol and triglycerides to the viral-replication complex.(7a) Through interactions with the plasma complement regulator C4BP, NS1 binds to the cell surface through N207.D4 would be targeting and blocking the active glycosylation site at N207. (7b) N130 stimulates the oligomeric assembly of secreted NS1 along with the stabilization of the secreted hexamer.D3 would be targeting and blocking the N130 site to disrupt the hexametric stability of the NS1 protein vital to viral stability, contagiousness, and antigenicity.
prediction of various toxicity endpoints ProTox-II (Banerjee et al., 2018) was utilized.The inhibition constant (K i ) was also calculated to measure the predicted inhibitor's inhibition potency.It is the concentration required to induce half the maximum inhibition of the drug.It was obtained from the binding energy (DG) using the formula: , where R is the universal gas constant (1.985 � 10 À 3 kcal/mol K) and T is the temperature (298.15K) (Ortiz et al., 2019).

Final molecular docking run
The final molecular docking run was carried out among the top 50 screened library that was extracted from the SwissADME website.The protein preparation as well as the ligand preparation was carried out in a manner similar to what has been discussed in section 2.5.After the final docking, the selection of the top three leads were assessed and extracted.The top lead was chosen as the final ligand based on its high binding affinity of the ligand to the protein and a better ADMET profile on account of its pharmacodynamics and pharmacokinetics.

Molecular dynamics simulation of unbound NS1 protein
MD simulation of the unbound protein structure was carried out via GROMACS v.2019.2(Abraham et al., 2015;Spoel et al., 2005).
The input parameters chosen during the submission of the unbound protein included the force field, utilized as CHARMM27 (Bjelkmar et al., 2010).The water model chosen was SPC, cubic box type was opted and salt type was NaCl.Neutralization was required along with the addition of 0.15 M salt.The energy minimization parameter included the steepest descent with 1000 steps (Chen et al., 2013).Equilibration and MD run parameters covered equilibration type as NVT/NPT, temperature as 310 K (physiological temperature of the human body to maximally mimic the natural environment of inhibitor action), pressure was set to 1.0 bar, the algorithm for MD integration was leap-frog (Van Gunsteren & Berendsen, 1988).The total simulation time for the system was 100 nanoseconds.

MD simulation and data analysis of bound proteinligand complex
We employed GROMACS simulation to perform fully solvated molecular dynamics simulations and trajectory analysis.The complex proteins along with the topology files prepared using Prodrg (Sch€ uttelkopf & van Aalten, 2004), were uploaded.The input parameters chosen during submission of the each complex protein included the force field, utilized as GROMOS43a1 with the rest of the parameters being similar to that of unbound NS1 protein simulation.

MD data analysis
The results obtainable from MD simulation are RMSD, RMSF, R g , and SASA.The structural fluctuations are visualized in the form of the Root-Mean-Square-Deviation (RMSD) which is a measure of the average instantaneous displacement of the atoms superimposed to the reference structure (Kufareva & Abagyan, 2012).Protein residual fluctuations given by RMSF calculate the displacement of one particular atom, or a group of atoms, associated with the reference protein structure, averaged over all the atoms (Farmer et al., 2017).Estimation of the radius of gyration (R g ) is an indicator to predict the structural activity of a protein, i.e. its stability and flexibility (Lobanov et al., 2008).Solvent Accessible Surface Area or SASA of a protein is defined as the characteristic surface around a hypothetical center of a solvent (represented as a sphere) with the van der Waal's contact surface of the protein molecule (Ali et al., 2014).A hydrogen bond or Hbond is a strong non-covalent interaction that helps maintain the structural integrity of the protein (Chen et al., 2016).Two vital short-range interactions: Lennard Jones short-range potential and Coulombic short-range potential were determined to evaluate the binding affinity of the ligand to the protein (Asthagiri et al., 1999).Hydrogen and hydrophobic bond analysis are crucial in drug designing because of their strong influence on the binding specificity of the ligand towards the target molecule (Pantsar & Poso, 2018).Ligplot þ was utilized to analyze hydrogen bonding and hydrophobic interaction between the two entities.

Conformational analysis during the simulation
The conformations corresponding to 25, 50, 75, and 100 ns were extracted in both the bound and unbound forms of NS1.The structures were then analyzed based on their respective secondary structure composition with an aim to predict the conformational alterations induced due to the binding of each of the ligands to the NS1 receptor protein.The webservers that aided this purpose were 2Struc (Klose et al., 2010) and STRIDE (Heinig & Frishman, 2004).

Binding affinity evaluation of the complexes throughout the simulation
The first 20 minima were extracted from each trajectory of the four complexes.PRODIGY-LIG was employed to calculate the binding affinity of the minimia conformations (Vangone & Bonvin, 2017).The resulting binding affinity predictor models DG score (Equation 1) and DG prediction (Equation 2) were employed for ranking ligands and predicting the affinity: (Vangone & Bonvin, 2017) (Vangone & Bonvin, 2017) where AC CC , AC NN , AC OO , and AC XX stands for the ACs between Carbon-Carbon, Nitrogen-Nitrogen, Oxygen-Oxygen, and between all other atoms and polar hydrogens, respectively.

Dynamical fluctuations-correlation and communication network analysis
Proteins' functional qualities are determined not only by their rather inflexible overall structures but also, perhaps more critically, by their dynamic properties.The residual fluctuations are inevitably one of the most crucial variables influencing the dynamic characteristics of a protein.The dynamic cross-correlation matrix (DCCM) technique has been widely used to assess and quantify the correlation coefficients of movements between atoms (Tong et al., 2021).Furthermore, to observe whether these correlations affected protein flexibility, PCA of the unbound and bound complexes was conducted.
To designate accessible motions over a broad range of time scales and spatial scales, protein conformations are best characterized by a vector space that extends to a large number of dimensions.PCA is a linear transform that extracts the most vital elements in the data using a covariance matrix or a correlation matrix constructed from atomic coordinates that define the accessible degrees of freedom of the protein, such as the Cartesian coordinates that outline atomic displacements in each conformation comprising a trajectory (Sittel et al., 2014).PCA was carried out with respect to the C-alpha atoms of the protein in ligand-bound conditions.GROMACS inbuilt modules were employed through Galaxy webserver (Senapathi et al., 2019) to calculate PCA and DCCM of the complexes using the solvated.groand xtc files of the protein-ligand complexes.

Total energy profile analysis
The total energy of a protein can be defined as the sum total of its potential and kinetic energy (Guan et al., 2011).The total energy of a protein primarily corresponds to its spontaneity of folding.For a 3D protein structure, the various components of total energy are electrostatic energy, covalent bonding energy, and Van der Waals energy.Conformations that deviate maximally from their native state have higher energy values.The conformation with the lowest energy represents its proximity to a supposed-native state.The total energy of a protein conformation s i is given in Equation (3): (Taylor & Jaffe, 2018) where, x i 2 R n represents the collection of the energy terms for s i , and w 2 R n stands for the weight coefficients.Ideally, the total energy of a protein decreases on complexation.Total energy analysis was made using OriginLab (v2023) using R g values to represent the total energy surface and RMSD and R g values for plotting the total energy landscape.Initially, we used Galaxy webserver (Senapathi et al., 2019) to tabulate the total energy of each conformation at a given instant of time.Extraction of the lowest energy conformation of the unbound NS1 protein was performed which was further utilized for comparison of the lowest energy conformations of the four respective complexes.
The entire methodology applied in the study has been drawn in Figure 2.

NS1_DENV2 structure extraction and analysis
The DENV NS1 protein is a highly conserved 46 kDa, 352 amino acid length glycoprotein consisting of 12 conserved cysteine residues and 2 N-linked glycosylation sites (Flamand et al., 1999).The unique feature of the said protein is to repeatedly polymerize and exists in various monomeric, dimeric, and hexameric forms (Akey et al., 2015).The polymerization and distribution of NS1 which is important for its functionality, depends on its glycosylation state, which ensues in the ER lumen.It has been determined that the three particular cysteine residues at the carboxy terminus of NS1(namely C313, C316, and C329) are essential for the protein's functioning, which is mediated by dimer formation.Moreover, N-linked glycan at the residues 130th position of the NS1 protein is critical for hexameric structural stability, impacting the binding affinity to the complement host cells, while the glycosylation at the 207th position is important for secretory purposes and mediating NS1 endosomal transport during clathrin-mediated endocytosis.

Protein cavity, pharmacophore detection, and virtual screening
Sixteen ligandable cavities were detected by the Cavity module run.Four cavities namely Cav 7, Cav 10, Cav 14, and Cav 15 are selected in correspondence to the presence of residues with following significances: Cav7 contains cysteine residues at 313th (C10), 316th (C11), and 329th (C12) position of the NS1 protein, which is responsible for intramolecular disulphide linkages, inducing stability to the molecule during homodimer formation, Cav14 contains the C55 residue that stabilizes the monomeric forms, whereas, Cav10 and Cav15 carries the conserved glycosylation sites of N130 and N207, respectively.The cavities prediction corresponding to the 7th pocket, 10th pocket, 14th pocket, and 15th pocket exhibited maximum pK d values of 6.99, 6.24, 5.81, and 5.77, respectively.A suitable good binding site is characterized by pK d with values around and >6 (Williams et al., 2003).Thus, all four cavities exhibited the properties of good binding site (Figure 3).The details of selected cavities are entered into CavPharmer and pharmacophore features like hydrophobic center, hydrogen bond donor, hydrogen bond acceptor, and charge center relevant to each cavity were generated.The results generated by CavPharmer are listed in Supplementary Table 1.
All these characteristics are crucial for the activity, affinity, contact, and binding efficacy of these ligands on our targeted cavities.Chemical library databases like the ZINC database are screened using the generated pharmacophore features to identify novel chemicals that have comparable therapeutic and interaction profiles matching our targets.Following the application of the filters, the ZINCPharmer search engine generated four ligand libraries for four cavities.These ligands obtained are likely to display a similar interaction with the selected cavities as they contain all the pharmacophore characteristics of the chemicals utilized to create the pharmacophore.These were finally docked onto the NS1 receptor.After the virtual screening, the top 100 ligands from each library were selected and validated as per their ADME properties.The top 50 validated ligands from each library were subjugated to the final molecular run and the top ligands were selected.The top ligands along with their binding affinities are tabulated in Table 1.

Root mean square deviation (RMSD) analysis
RMSD is accountable for characterizing global changes all along the protein chain (Kufareva & Abagyan, 2012).The average protein backbone of superimposed protein structures is measured by root-mean-square deviation (RMSD) to compare the structural deviance (Kufareva & Abagyan, 2012).
A constant RMSD curve is indicative of the fact that the protein has reached an equilibrium state (Guo et al., 2015).
From Figure 4(a), the conformational variations of the NS1 protein backbone with D1 throughout a 100 ns MD trajectory were plotted to determine the stability of the protein-ligand complex.The RMSD in both the bound and unbound conditions was calculated using the backbone atoms of the protein, and comparisons were made between the two curves.The high amount of fluctuations in the unbound protein conformations is owing to its instability.In the unbound state, the curve contains points of numerous crests and peaks.It fails to attain a steady and equilibrated conformation, lacking structural stability.During the 100 ns simulation, the unbound protein fluctuates twice: once from 0 to 40 ns and again from 70 to 100 ns.Significant RMSD variation was observed in the curve during the 0-40 ns interval, however, the RMSD variations between 70 and 100 ns interval were less pronounced than those between 0 and 40 ns, but they were nonetheless consequential.Such fluctuations are suggestive of a protein that lacks conformational stability.The reason might be the high content of destabilizing random coil percentage and 3 10 helix composition in the protein during the simulation.
The protein in the bound condition (D1 ligand and NS1), furthermore, exhibited inconspicuous RMSD deviations throughout the simulation.Such non-fluctuating nature indicates a stable protein-ligand complex, which means that the association of the D1 ligand with the NS1 protein, decreased the overall flexibility of the protein by inducing a remarkable amount of steadiness and compactness in the structure of the complex, therefore, giving the protein an equilibrated state.
The RMSD curves for unbound_NS1 and D2Complex were superimposed (Figure 4(b)) to elucidate the difference in stability between the two structures.The RMSD plot curve for unbound protein emphatically fluctuates towards the initial course of simulation in comparison to D2Complex which is almost stable throughout the simulation runtime.This data is suggestive of lower structural rigidity in unbound_NS1 with respect to its D2complex.It can be concluded that after the binding of NS1 to D2, there occurs an overall flexibility in the protein structure thereby reducing the extent of fluctuations due to excessive mobility of atoms.The principle behind this is that the stability of interaction rises on ligandtarget protein binding.
Superimposed RMSD curves of unbound_NS1 and D3Complex (Figure 4(c)) show that the binding of a ligand to the target protein results in the decrease of the RMSD fluctuations in the protein.The D3Complex exhibits minimal fluctuation in the RMSD which is a symbol of more stability and steadiness compared to the RMSD of the unbound one.The reduced amount of RMSD fluctuations depicts the gradual increase in the overall stability of the protein when it binds with D3 ligand.This is induced due to the strong interaction between the ligand and the protein which results in the formation of a very stable and compact complex to attain its equilibrium state and remain consistent throughout the simulation.
Comparing unbound_NS1 and D4Complex (Figure 4(d)), it can be perceived that the RMSD curve of D4Complex fluctuates less, maintaining an almost steady and plateau-like RMSD curve throughout the simulation.The minimum variations represent the steady and stable nature of the complex.Upon ligand binding the protein becomes more stable, compact, and less flexible.This compact nature of the protein could be a result of strong, stable interaction between the entities of the complex resulting in a well equilibrate conformation.

Radius of gyration (R g ) analysis
The radius of gyration (R g ) indicates the mass weight analysis of a collection of atoms at their centre of mass as well as the protein's compactness.It depicts the equilibrium conformation of the total system (Lobanov et al., 2008).A variation in R g curve implies improper protein structure folding (Kony et al., 2007).The R g value of the NS1 protein in its both bound and unbound state was extracted during 100 ns MD simulation run.Lower the value of R g , constricted is the atomic packing of the protein residues.The average R g of unbound protein is 2.2023 nm, whereas for D1Complex, D2Complex, D3Complex, and D4complex the average R g amounted to 2.0929, 2.1075, 2.1270 and 2.0577nm, respectively.
The radius of gyration of the unbound_NS1 decreases by 0.1094 nm upon binding of the ligand D1 (Figure 5(a)), resulting in an increase in the overall compactness of the protein, and a decrease in flexibility.The R g curve of the unbound protein shows noticeable fluctuations during the initial 40 ns of the simulation; on the other hand, the R g curve of the bound protein exhibited minimal fluctuations only during the initial 20 ns of the simulation.Thus, it is suggested that D1 has the potential to strongly bind to the protein by enhancing the protein's compactness owing to conformational changes.Superimposition of R g curves of unbound_NS1 and D2Complex (Figure 5(b)) explains the decline in the structural flexibility of the bound complex due to fall in the R g value by 0.0948 nm with respect to its unbound state.A drastic difference in fluctuations of R g curve between unbound_NS1 and D2Complex is observed from the R g plots of the simulated structures.Evidently, D2 would serve as a proficient drug to inhibit viral replication mediated through our target protein NS1.
Association of the D3 ligand with the protein consequently caused a higher degree of compactness in NS1 protein structure which is characterized by the significant decline in the average R g value by 0.0753 nm with respect to its unbound state.Figure 5(c) projects the superimposed curves of the R g values generated for both the unbound and bound complex.The curve for D3Complex fluctuates lesser in contrast to the curve for the unbound protein, depicting the gain in conformational stability due to the steady binding and strong interactions between the ligand and the protein leading to a stable, well-equilibrated complex.
From Figure 5(d), a comparison between unbound_NS1 and D4Complex clearly depicts the decrease in R g of complex in correspondence to the unbound protein.Both the R g curves for unbound and D4Complex exhibit steady and stable regions over the major of the 100 ns simulation.However, the slight fluctuations of unbound protein in the initial 30 ns could be an indication of structural instability.Consequently, D4Complex maintains a constant, plateau-like region with only minor fluctuations during the entire simulation.This demonstrates that D4Complex (formed between D4 ligand and NS1) attains a properly folded, compact, stable state and maintains it throughout the simulation.D4Complex displays a decrease of 0.1446 nm in R g in the bound state.This indicates the increase in compactness, stability, and decrease in flexibility of the protein upon ligand  -benzoyl-7 0 binding owing to its higher binding affinity between the entities of the complex.

RMSF analysis
RMSF is indicative of the residual fluctuations in the generated protein conformations as a result of macromolecular motions over the entire period of simulation and hence the conformational flexibility/stability.RMSF depicts local changes all along the protein chain (Bornot et al., 2011).For the inhibitor obtained against the second cavity, a comparison of RMSF values of the interacting residues of both the unbound and bound forms of NS1 is done by superimposing the RMSF graphs of unbound_NS1 and bound NS1  This fall in the RMSF values signifies that atomic fluctuation decreases on the interaction of the target protein with the ligand specifically due to contact at the related residue(s).Herein, the RMSF result analysis can be viewed in favor of the target ligand stability giving rise to a conformation with low flexibility.Consequently, such a drop in the atomic fluctuations of the interacting residues signifies that the ligand formed a secure interaction giving proper stability to the protein structure with reduced flexibility.Furthermore, the standard deviation (SD) of the RMSF values calculated for the unbound NS1 protein is 0.13484391nm which is higher than all of the SD values of the four chosen inhibitors, i.e. 0.067684, 0.077655, 0.068684 and 0.077381nm for D1, D2, D3 and D4, respectively.This data suggests that there is lesser fluctuation in the bound protein-ligand complexes than in the unbound protein.Conclusively, an increase in the stability or compactness in the bound structure is observed which is an advantage for considering these four ligands as a potential drug against dengue.

SASA analysis
The amount to which an amino acid interacts with the solvent is related to the surface area exposed to these environments (Ali et al., 2014).The better the protein-ligand interaction, the lower the SASA value and the more stable the complex will be (Du et al., 2016).The average SASA of unbound protein is 166.713nm 2 , whereas for D1Complex, D2Complex, D3Complex, and D4Complex the average SASA amounted to 147.324, 146.5092, 147.894, and 147.352 nm 2 , respectively over the 100 ns simulation.The SASA graphs were superimposed from the 100 ns simulation runs of both the bound and unbound_NS1.
The curve obtained for the unbound state exhibits prominent oscillations, but the curve obtained for the D1Complex has minimal fluctuations (Figure 7(a)), indicating the optimal and constant interaction of the D1 ligand with the NS1 protein.Over the 100 ns simulation period, the unbound_NS1 generated an average SASA value of 166.713 nm 2 , whereas the D1complex generated an average SASA value of 147.324 nm 2 .D1 interaction with the NS1 protein resulted in a 19.389 nm 2 decrease.The NS1 protein's 199th Glycine residue was discovered to create a hydrogen bond with the D1 ligand, resulting in an efficient drop in the SASA value.The NS1 protein residues Ile21, Thr22, Thr27, Trp28, Tyr32, Trp201, Lys272, and Leu273 were shown to be engaged in hydrophobic interactions with the D1 ligand and exhibited a significant decrease in SASA value when compared to the unbound state (Supplementary Figure 1(a)).This infers that the D1 bound NS1 protein interacted more effectively than the unbound_ NS1 protein, suggesting a more stable and compact complex devoid of flexibility.
Both the SASA curve in D2Complex and unbound state remain majorly non-fluctuating throughout the simulation (Figure 7(b)).The average SASA of D2Complex as determined was 20.2036 nm 2 less than NS1 in its unbound form.The highest obtained SASA value for unbound NS1 is 179.408nm 2 whereas the maximum SASA value for D2Complex is 171.659nm 2 .All this data suggests that solvent accessibility of NS1 protein declined on complexation.Interacting residues Trp28, Thr29, Arg62, Leu66, Lys69, and Leu169 show a steady decline in SASA values due to receptor domain-ligand binding (Supplementary Figure 1(b)).The decline in the SASA values of the interacting residues corresponds to their decrease in RMSF values due to a decrease in the atomic fluctuations caused by restricted molecular motions.
The average SASA for the D3complex is 147.894nm 2 over the 100 ns simulation.Both the SASA curve in bound and unbound state displays mostly non-fluctuating parts throughout the simulation (Figure 7(c)).A gradual decrease is observed in the SASA value when the D3 ligand binds to the protein.This high decrease in SASA value is due to the compactness of the bound complex which arose due to the active and steady interaction with the residues of the D3Complex.The interacting protein residues Gly9, His26, Trp28, Arg62, Asn65, Asp92, Ile93, Lys94, Gly95, Ile96, Met97, Ile218, Val220, His269, Leu270, Gly271, and Lys272 showed a significant decrease in the SASA values (Supplementary Figure 1(c)).Thus, it can be concluded that the two entities of the complex developed a very stable and strong interaction between them which lead to the formation of a more compact protein structure.
Both the SASA curves for unbound and D4Complex (Figure 7(d)) exhibit steady and stable, non-fluctuating regions over most of the 100 ns simulation.D4Complex displays a decrease of 19.361 nm 2 in SASA in a bound state in comparison to the unbound molecule.This high decrease in SASA value in the bound state denotes a high number of residual interactions in the complex.Residues like Leu13, Lys14, Thr22, Asp23, and Lys189 which were involved in maximum hydrophobic and hydrogen bond interactions displayed a decrease in their SASA values after complexation (Supplementary Figure 1(d)).This signifies the formation of stable, compact, steady complex characterized by strong interactions between the entities.

H-bond analysis
Hydrogen bonds are one of the key non-covalent interactions that contribute to the stability of protein-ligand complexes.H-bonds have also been potentially identified as critical to protein folding simulation (Gao et al., 1988).During the protein-ligand binding process, besides the existing intramolecular H-bonds of these two molecular entities, new H-bonds are created at the protein-ligand interaction interface (intermolecular interactions).The formation of such hydrogen bonds renders drug-like molecules with the potential to execute specialized interactions as well as an affinity for the target protein.The amount of hydrogen bonds formed (both intramolecular and intermolecular), during the 100 ns simulation run are utilized to assess the stability of the complex, by superimposing both the bound and unbound curve.The protein maintains a steady curve in both its bound and unbound states (Figure 8).The average number of hydrogen bonds for the unbound protein was determined to be 214 whereas the average value for the D1Complex, D2Complex, D3Complex, and D4Complex amounted to 221, 221, 217 and 219, respectively.This increment in the average number of hydrogen bonds upon ligand association indicates the development of hydrogen bonds in the complex, resulting in high binding affinity and strong interactions of the ligand with the protein.

Evaluation of conformational changes during the simulation
Analyzing the conformational changes of the secondary structure during the simulation is critical for the protein complex stability.For the unbound NS1 protein (Supplementary Figure 2), all the secondary structural elements fluctuate with almost no constant trend.The sudden variances in the conformational states indicate that the unbound protein fails to achieve a well-equilibrated stable structure over the entire simulation time.The poor arrangement of the secondary elements in correspondence to the fluctuating RMSD and R g value further strengthens the conformational instability of the NS1 protein during the simulation.From Supplementary Figure 3, it was observed that in D1Complex isolated b bridge, extended bridge and turns composition displayed fluctuations throughout the simulation.Extended bridge (21.2 !19.8%) and turns composition (35.5 !32.1%) showed a noticeable decrease during the second 25 ns of the simulation period.The isolated b bridge exhibited fluctuating patterns with rapid increase and decrease.A steady trend was observed in the case of 3 10 helix until a sudden rise (1.7 !3.4%) was recorded at 75 ns.a-helix followed a steady trend till 50 ns, but a slight inclining trend was noticed from 75 ns, whereas, random coil showed an increase in the second 25 ns of the simulation but exhibited a declining trend from the 75 ns.Such a diminution in the percentage of the random coil characterizes a reduction in the flexibility and henceforth, the dynamic motion of the protein complex.The trend demonstrated by the two forms of helix structures and also the RMSD fluctuations, implies that a transition of 3 10 helix to a-helix (Armen et al., 2003) has occurred, thereby increasing the interaction of the ligand with the protein.
As depicted by the STRIDE analysis, the extent of increase of b-pleated sheets in D2Complex (Supplementary Figure 4) varies as in 20.3% for 25 ns, 23.5% for 50 ns, 22.6% for 75 ns, and 22.1% for 100 ns.The percentage contribution of a-helix is almost constant all along the MD simulation ranging from 6.9 to 7.2%.This constant trend of a-helix is influenced by inter and intramolecular hydrogen bonds developed during the formation of the secondary structure.This a-helix and b-pleated sheets procure a well-built structure of D2Complex.Lower variation in the structure accounts for structural rigidity as an outcome of the higher level of receptor-ligand interaction.Lack of variation implies accession of an equilibrated conformational state and attainment of a plateau-like region in the RMSD curve.A fair amount of random coil contribution is observed in the complex ranging from 26.6 to 32.4% (30.7% for 25 ns, 26.6% for 50 ns, 30.1% for 75 ns, and 32.4% for 100 ns).The presence of 3 10 -helix (g) in the bound complex conveys helix !coil transition (Armen et al., 2003).It is probable that conversion of the 3 10 -helix (an intermediate structure between a-helix and random coil) to a-helix (g !h) rather than to a coiled structure (g !c) which subsequently favors structure stabilization (Armen et al., 2003).Turns and Isolated b bridges have also been traced in the concerned structure.The stability of the protein is conferred by the presence of strong inter and intermolecular hydrogen bonds via helix and extended strand evolution.The flexibility or frequent fluctuations occur as an outcome of the presence of loops and coils in MD simulation trajectories.Throughout the simulation, numerous patterns are being observed in the secondary structural conformations of the D3Complex.From Supplementary Figure 5, it is observed that extended strands, isolated b bridge, and turns are continuously fluctuating throughout the simulation.Initially from 25 to 50 ns, there was an increase in isolated b bridge percentage, with an instantaneous decline (4.6 !2.9) during the last 50 ns.The composition of the extended bridge and random coil got decreased during the initial 25 ns but remains significantly constant for the rest of the simulation.3 10 helix composition fluctuates throughout the simulation with a declining trend with a major decline during 50 and 100 ns (4.6 !1.7 and 4 !0.9%, respectively) of the simulation.For the initial 50 ns a-helix shows a declining trend but later on remains steady and non-fluctuating throughout the simulation.All together a proportionate and steady composition of a helix, isolated beta bridge, turns, extended bridge, and random coils contribute towards the formation of a very steady, well-equilibrated, and stable interaction between the ligand and the protein receptor which ultimately gives rise to a stable and compact protein-ligand complex.
Observation drawn from Supplementary Figure 6, the secondary structure conformations of D4Complex show various trends throughout the simulation.Extended strands, isolated b bridge, and turns composition fluctuates throughout the simulation.A declining trend is observed in isolated b bridge percentage till 75 ns, with a sudden increase (2.6 !4%) during the last 25 ns.Extended strand and turns composition decrease during the first 25 ns but later shows an inclining trend during the rest of the simulation.a-helix and 3 10 helix composition remain majorly steady throughout the entire simulation.A declining trend is observed in the percentage of random coils.The steady composition of a-helix, an increasing proportion of extended strands and turns along with decreasing random coil represents decreasing flexibility, increase in stability, steadiness, compactness, and organization of the complex structure with stable interactions occurring in the protein structure.

Dynamical cross-correlation matrix (DCCM) analysis
The relationship between the motion correlation and the structural change has been investigated to fully grasp the dynamic communication induced by the hydrogen and hydrophobic binding in the protein-inhibitor complexes (Zhang & Lazim, 2017).The correlative motions of the four different systems were analyzed by the dynamic cross-correlation matrix (DCCM).The cyan color indicates a strong positive correlation; however, the pink color indicates a negative correlation.Residues exhibiting a concerted movement in the same direction are said to be correlated, i.e. the parallel movements represent correlated motion.Conversely, the regions that show fluctuation or are anti-parallel represent the residues that are anti-correlated.The correlations occur as a result of specific contacts or long-range interactions.Interestingly, all the protein-inhibitor complexes exhibited different movements than their unbound contender (Figure 9).The variance in the movements of the complexes is proposed to be induced by an allosteric transition.The DCCM of target regions of the protein (both unbound and bound form) are analyzed to understand the effect of inhibitors on the residual movement of the target region.For D1 (Figure 9(b)), the target region (C313, C316, and C329) undergoes a drastic change from its unbound form.In the unbound (Figure 9(a)) form, the target region is mainly interspersed with a deeper range of colors indicating strong correlations.It exhibits negative correlations (anti-correlated) with residues 160-210 and positive correlations with residues 210-280.In protein-ligand form, the same target region displayed majorly low anti-correlation with the residues 160-280.For D2 (Figure 9(c)) the target region (N130) displays both high correlated and anti-correlated movements with a residue range of 1-250 in the unbound form.In the bound form, the residues in the region exhibit high correlation with residues 1-120 range whereas, exhibits anti-correlated movements with residues 155-250.For D3 (Figure 9(d)), the target region (C55) displays anti-correlated movements with residues 40-50 and 250-300 in the unbound form.It also exhibits high positive correlations with residues 20-40.In the ligandbound form, the target regions show majorly low positive correlations with residues 20-50.It also displays insignificant or no residual correlative movements with 250-300 after ligand binding.For D4 (Figure 9(e)), the target region (N207) displays highly correlated movements with residues 1-230 in the unbound form.In complex form, the target region parallels the residual motion of the unbound state residues 170-230 but exhibits major anti-correlated movements with unbound residues 1-160.Further, in comparison to unbound residual movements, the dynamic matrix of the ligand-bound state of the protein was mainly dominated by white color which denotes no correlation between the residues.This depicts that residual movements are restricted in regions of the bound protein.This observation corroborates the fact that stronger ligand binding will induce compactness in the protein, thus overall movement of the protein will decrease (Limongelli et al., 2012).Conclusively, the loss of unbound correlations in the ligand form confers the stability of the protein-inhibitor complexes.

Principal component analysis (PCA)
Principal Component Analysis (PCA) is a multivariate statistical technique applied to systematically decrease the number of dimensions required to define protein dynamics through a decomposition process that screens observed motions from the largest to smallest spatial scales (Sittel et al., 2014).PCA helps develop a perception of the dominant motions in a protein structure based on the eigenvector values collected from the covariance matrix of the simulation followed by filtering the trajectories with their respective principal components or eigen vectors (Maisuradze et al., 2009).We examined the first 20 eigenvalues and the first two principle components, PC1 and PC2, to analyze about the systems' dynamic movements and fluctuations.The eigenvalue plot against the proportion of variance of the four complexes is plotted in Supplementary Figure 7.A comparative analysis of the 2D-PCA was performed for PC1 and PC2, to review the diffusive properties of the NS1 protein in its unbound and D1, D2, D3, and D4 during various stages of its folding pathway.
Unbound NS1 exhibited dispersed motion with a higher degree of atomic variation along both PC1 and PC2, whereas D1Complex (Figure 10) exhibited confined motion along those two principle components.The majority of conformational movements were restricted to a limited part of the broad conformational space for D1Complex, whereas the unbound state demonstrated a greater degree of travelling throughout the conformational space.The overall analysis of dynamic motions would suggest that binding of the ligand .In all cases, correlations are represented in the color spectrum, from cyan to pink.Parallel motion of the residues indicates correlated movement, however, antiparallel motion of the residues exhibits anticorrelated movement.The positive correlation between residues, that is, residues move in the same direction, is represented as cyan color.Alternatively, residues having opposite motion are denoted as the anti-correlative motion marked as pink color.
accounts for the change in protein atomic motion with lowering of structural flexibility of the protein and acquisition of lesser phase space of the bound protein (David & Jacobs, 2014).The ligand's association enhanced the efficiency of the protein to attain its native state at a more rapid pace than the unbound protein would take to reach its folded state, therefore leading to the accomplishment of a compact, stable D1Complex.
Visual representation illustrates the compactness of D2Complex and more structural rigidity in the bound state than in the unbound form (Figure 11).It can be observed from the graph that in all the quadrants, the motion of the atoms and molecules of the complex is more restricted in the conformational space when compared to the unbound protein NS1.In other words, a clustered appearance of the D2Complex is seen wherever present.The reduction in intermolecular motions causing the structural rigidity of the complex reinforces the stabilizing nature of NS1 bound to D2Complex rather than its unbound form.Over and above, PCA involves the extraction of the most integral variables (in the form of components) from a large set of variables that are present in a data set.The objective is to incorporate as much information as possible.Corroborative to this fact, it can be said that the availability of a larger dataset would have made the visualization of a more stabilizing nature of D2Complex more equitable.
The highly haphazard and scattered motion of the unbound NS1 was noticed along the XY plane of two principle components accounting for its higher degree of dynamic motion.However, the plot obtained for the D3Complex exhibited highly constricted movement of all its conformations, throughout the vast conformational space (Figure 12).The central clustering and the majority of structural movement of the D3Complex were confined to an extremely limited area, in contrast to the unbound complex, indicating a decrease in the overall motion of the protein upon ligand binding (David & Jacobs, 2014).This signifies that ligand binding induces the formation of a stable, well-  equilibrated complex which is more steady and compact displaying lesser fluctuations.
To represent the smaller atomic motions of a protein, PCA includes both larger atomic motions and smaller subspace motions.PCA is not restricted to harmonic motions because it does not assume any harmonicity.Indeed, the conformational changes that can be examined as a result of PCA can diverge significantly from the harmonic assumption since PCA is independent of the model invoked during the simulation to construct the trajectory (Liu et al., 2018).Though comparing the unbound and bound D4Complex macromolecule from Figure 13, it seems bound protein conformational motion occupies more space in the 2D matrix, but understanding the PCA of a protein is much more complex than that.The overall major large-scale motion of the protein complex is more restricted with fewer fluctuations and deviations than its unbound contender in the second and third quadrants of the PCA matrix signifying a strong interaction between the D4 and unbound protein.Before analyzing the upsurge in the motion of the complex in the first and fourth quadrant subspace, we have to keep in mind that individual PCA mode directions are subject to errors related to finite sampling of conformations (Liu et al., 2018).The empirical C-matrix should provide a reliable approximation of the population C-matrix (infinite samples).In reality, the existence of outliers in a dataset can have a significant impact on PCA (Liu et al., 2018).The primary issue is that the outliers can distort the first few mode directions.Additionally, as long as biologically important motions are described using a superposition of a small set of dominating modes (second and third quadrant) rather than concentrating on one subspace motion, this form of intrinsic mistake does not present a significant issue.Our discussion was focused on the first two principal components; thus, it can be hypothesized that the upsurge of protein motion in the small subspace might be its failure to attain full conformational stability during the first two PCs but an increase in the mode number could increase the core part of this subspace to become stable against sampling noise.

Total energy analysis
The translation and rotation (or transformation) of one protein chain (here the ligands) with respect to the other chain (NS1 as the receptor protein) drive the energy function of the total system.The total energy or internal energy of a protein represents the summation of total kinetic energy and potential energy associated with the motion of the atoms, and vibrational and electrical energy of the molecules in the complex.A macromolecule must possess a low internal energy to be considered as a stable structure (Deller et al., 2016).For a protein molecule to attain thermodynamic stability, it must fold into its native structure, manifesting a state of minimum energy among all possible conformations.The native structure of a protein molecule, is thereby, considered to be in its utmost biologically active form, due to its structural stability.Construction of total energy landscape and surfaces using the simulation trajectories is one way of determining consistent complex conformations.The energy surface and landscape for both complexes were obtained by mapping the total energy to the corresponding radii of gyration and C-alpha atom RMSD.The total energy surface value is interpreted from the graph, in which, the x-axis represents the energy in kJ/mol and y-axis represents the radius of gyration values in nm.For the unbound NS1 protein (Figure 14), the lowest energy (À 712,747.875 kJ/mol) was extracted at R g value 2.17827 nm.
However, for D1Complex (Figure 15), the lowest energy (À 804,978.1875kJ/mol) extracted at R g value 2.09747 nm, showed a noteworthy decrease as corresponded to the unbound state.Furthermore, conformations of various energy levels for the unbound state were seen to be more scattered in the energy surface compared to that of the bound state, which is representative of the fact that D1 induced considerable steadiness to the protein.Additionally, a less scattered and extremely compact representation of the energy landscape for the D1Complex was seen, suggesting that most of the conformation attained the lowest energy minima making D1 a potential drug.
For D2Complex (Figure 16), the lowest energy conformation (where energy ¼ À 805,469.1549KJ/mol) is obtained at R g ¼ 2.09509 nm.Evidently, the complex conformations are products of less fluctuating residues due to less variation in their energy values as obtained from the scatter diagram, and hence D2Complex follows a stabilizing trend in view of structural orientation.In contrast to unbound NS1, a less outspread structure implies that most of the conformations are liable for attaining the lowest energy minimum within a given timescale to reach a state of equilibrium.
The lowest energy of the protein in its D3Complex (À 805,553.3125kJ/mol) was extracted at R g value 2.1591 nm.The graph for the protein's unbound state depicted scattered trends for all the conformations on the energy surface, while, for the D3Complex, compactness at a higher degree was noticed over the energy surface, establishing a profound interaction of the inhibitor D3 which thereby induced great stability to the protein (Figure 17).Also, the energy landscape for the D3Complex displayed a potential compact structure for all the conformations as compared to the unbound state.
D4Complex displayed a minimum energy conformation with energy value À 805,173 kJ/mol at 2.077 nm radius of gyration.The complexes display a great number of  conformations occupying a lower energy state than its unbound contender, owing to its more stable and compact structure.Comparing the energy surface and landscape of both the unbound and bound forms in Figure 18 displays less scattering of D4Complex with a higher number of conformations achieving a low energy state at a lower R g range than its unbound state.
Overall, the results are corroborated in the formation of a narrower basin in the energy landscape of all the proteininhibitor complexes where red to blue color denotes a high energy state to a lower energy state of the conformations.The conformational energy, upon the association of the ligand, is reduced to a great extent, which is characterized by an alteration in the color of the energy landscape in the unbound and bound states.For the unbound state, the major conformations fail to attain the lowest energy as depicted by the distribution of rainbow spectrum.Whereas, for complexes, the landscape is dominated by colors representing lower energy structures.This denotes how maximum conformations achieved stability at a lower energy state in its complex state.Thus, binding of the ligands to the receptor protein led to an increase in overall stability, steadiness of protein, and a decrease in the free energy of the system showing that ligand binding is a favored process in comparison to its unbound form.Conclusively, D1, D2, D3, and D4 as ligands can be termed to be energetically favorable as a dengue inhibiting drug.
A comparison of the minima conformations of the unbound and bound NS1 reveals a considerable amount of changes in the conformational state of the proteins.For D1Complex minima, an upsurge of turns is noticed along with a decrease in destabilizing random coil composition in comparison to unbound minima.The overall composition of turns increased by 4.8% and coil decreased by 0.6%.These display an increase in compactness of the protein with a subsequent decrease in flexibility of the protein.A significant change is observable in the percentage of turns in D2Complex (30.7-35.5%).An increase in the percentage of turns (b-turns) signifies a stabilizing structure due to the presence of hydrogen bonds between the first and the fourth residue of the b-turn (Marcelino & Gierasch, 2008).This stability is accompanied by the decline in the percentage contribution of the random coil from 34.1% in the unbound form of NS1 to 33.5% in its bound state.Lessening of random coil denotes decreasing structural disorderliness (Khrustalev et al., 2013).An increase in 3.1% isolated beta bridge composition in D3Complex with respect to its unbound minima indicates a better protein organization.Moreover, the percentage of extended strands increases by 2.9% with a subsequent decrease of 2% random coil in the D4complex than its unbound contender, demonstrating a significant enhancement in the structure's stability since b-Sheets are substantially more stable than turns and coils (Cheng et al., 2013).In both D3Complex and D4Complex, a higher proportion of 3 10 -a-helix composition is obtained in correspondence to their unbound conformation.3 10 -Helices have been proposed to be intermediates in the folding/unfolding of a-helices.In the helix-coil theory, it had been suggested that a/3 10 /coil conformations remain at equilibrium and 3 10 often acts as the transition between helices and coil  ( Rohl & Doig, 1996).Molecular dynamics studies of helices in peptides have proposed that a helices and 3 10 co-occur along the same peptide and that interconversion between the two conformations occurs rapidly in the nanosecond timescale (Rohl & Doig, 1996).The presence of a higher proportion of 3 10 and a helices together often induces higher stability in the protein (Peterson et al., 1999).Their cumulative effect in the D3 and D4 complexes induces a stronger hydrogen network in the complex, bringing higher compactness, steadiness, and stability to the complexes.

Study of interaction energy of the complexes
The Lennard-Jones potential simulates the van der Waals interaction; however, the coulombic potential simulates the electrostatic interaction between protein and ligand atoms depending on their separation (Asthagiri et al., 1999).For a protein molecule to be in its fully folded state, it should project the lowest amount of interaction energy among all of its conformational states, owing to the attainment of its absolute stable configuration.For 100 ns simulation, the interaction energy graphs for both unbound and bound proteins have been constructed in Figure 19, respectively.The graph for the unbound NS1 displayed more variations, but for the bound system, the energy curve exhibits the trend of remaining in a stabilizing condition with limited atomic motion along the entire protein chain, with only minor deviations suggesting an efficient interaction of the ligand with the protein.The average total value of LJ-SR and C-SR for unbound is À 955,104.4688kJ/mol (Supplementary Figure 10) and for D1Complex, D2Complex, D3Complex, and D4Complex retrieved from the simulation, amounted to À 1,019,039.715,À 1,019,349.482,À 1,019,960.785,and À 1,019,296.614kJ/mol, respectively.A significant reduction in the average sum of the LJ-SR and C-SR interaction energy values denotes that the ligand brought an appreciable amount of stability to the protein, subsequently aiding the protein to gain its native state at a much more rapid degree.

Binding affinity analysis of the complexes
The average was calculated for the 20 minima conformations of the complexes (Figure 20) and tabulated in Supplementary Table 8.The average binding affinity obtained for the complexes is À 8.785, À 8.62, and À 9.385 kcal/mol for D1Complex, D2Complex, D3Complex, and D4Complex, respectively.The data calculated further co-related with the docking score (Table 1) evaluated formerly.Thus, the ligands display high binding affinity towards the NS1 targets.Therefore, the ligands can be proposed as potent inhibitors with high binding with NS1 active site saturation and good binding orientation.

Interaction pattern of the complexes
Hydrophobic interactions are thought to be one of the most critical elements in the stability of specific folded configurations in natural proteins (Meyer et al., 2006).Hydrogen bonds and hydrophobic interactions are the major classes of molecular interactions that optimally secure the energetically favorable leads (ligands) in their positions when in contact with the receptor surface (Patil et al., 2010).Hydrophobic contacts have exhibited their ability to heighten target receptor and inhibitor interaction in the complex.
For D1Complex (Figure 21) 16 hydrophobic interactions and 3 hydrogen bonds formed between the ligand and protein.Protein residues,specifically,Ile21,Thr22,Asp23,Val25,Thr27,Trp28,Tyr32,Asp197,Met198,Tyr200,Trp201,Ile217,Ile218,Lys272,Leu273,and Glu274 were discovered to be engaged in hydrophobic interactions with the ligand.Consequently, an augmented number of hydrophobic interactions is suggestive of the fact that the ligand exhibits a higher binding affinity towards the protein.Residues Asp24 and Gly199 were involved in hydrogen bond formation.Asp24 was involved in 2 hydrogen bonds with the ligand of bond length 2.87 Å and 3.21 Å.
For D2Complex (Figure 22) 12 hydrophobic interactions formed between the entities of the complex.The interacting residues Trp28, Glu30, Tyr32, Phe34, Arg62, Leu63, Leu66, Met67, Lys69, Gln70, and Leu169 formed hydrophobic contact with NS1.Naturally, these residues directly determine the nature of binding or binding affinity of the receptor with the active ligand.
A total number of 17 interacting residues were involved in the interaction with the D3complex (Figure 23).From them, two of the residues, i.e.Val220 and Met97 formed hydrogen bonds with the NS1 protein.The bond lengths of hydrogen bonds formed by the residues are 3.13 and 2.96 Å, respectively.The residues that formed hydrophobic bonds with the target protein are Gly9, His26, Trp28, Arg62, Asn65,  Asp92, Ile93, Lys94, Gly95, Ile96, Ile218, His269, Leu270,  Gly271, and Lys272.The interaction resulted in the gain of a more stable and steady state of the complex.

ADMETox evaluation
A good ADME profile of lead compounds has been seen to be associated with an increase in the efficacy and therapeutic action of the drug along with acceptably low toxicity levels, which ultimately results in their success in subsequent in-vitro and clinical trials (Guan et al., 2019).All the ADME data related to the four predicted inhibitors are present in Supplementary article.The data for physiochemical properties, lipophilicity models, water solubility, pharmacokinetics, drug-likeness, and medicinal chemistry for the four inhibitors are tabulated in Supplementary Tables 2-7, respectively.

Bioavailability, lead-likeliness, synthetic accessibility, and PAINS and Brenk filters
Drug-likeness is a quantitative parameter that measures a compound's oral bioavailability (Bickerton et al., 2012).The bioavailability score is calculated based on Lipinski's rule of five.Using semi-quantitative data, compounds are divided  into four likelihood score groups according to rule-based scoring, namely 11, 17, 55, and 85% (Martin, 2005).Our inhibitors have shown a bioavailability score of 55% which lies within the acceptable probability score (Supplementary Tables 2 and 6).It also can be interpreted that our inhibitors have passed the Lipinski rule of five.Synthetic accessibility evaluation is a process to assess the ease of synthesis of compounds (Baba et al., 2018).The lower the value of the synthetic accessibility score, the easier it is to synthesize the drugs (Ertl & Schuffenhauer, 2009).D1, D2, D3, and D4 exhibit a score of 2. 76, 4.78, 5.16, and 4.22, respectively.Conclusively, D1 poses as the easiest to synthesize among all the four inhibitors designed and discovered.Ease of drug synthesis ascends in the order D3-D2-D4 complex.
The Brenk method and PAINS (Pan Assay Interference Compounds) method are used to identify potentially problematic molecular fragments that might produce biological  activity results that are falsely positive (Alam & Khan, 2019;Dahlin et al., 2015).No alert was found for D1 signifying the presence of no problematic chemical fragment in it, D2 has alerts of phenol and ester groups, D3 has an alert of a nitro group and D4 has an alert for beta keto anhydride and nitro group.A compound's lead-likeness is predicted using factors like MW (250 g mol-1 MW 350 g mol-1), coefficients of octanol/water partition (XLOGP 3.5), and the total number of rotatable bonds.Only D1 fell under the set criterion for leadlikeness, whereas D2, D3, and D4 failed to reach the optimum value (Supplementary Table 7).

Absorption, distribution, metabolism, and excretion properties assessment of the top four ligands
The effective solubility of a substance in a non-aqueous media, or lipophilicity, is associated with pharmacological properties, including adsorption, distribution, metabolism, and excretion (Shekhawat & Pokharkar, 2017).The rate of solubility of the compounds in aqueous and non-aqueous conditions is a key factor affecting their absorption in their body (Shekhawat & Pokharkar, 2017).
The efficient solubility of a compound in non-aqueous medium is termed as its lipophilicity capability and is  correlated to various models of drug properties, including toxicity, distribution, adsorption, and metabolism.Five available predictive models, including iLOGP (implicit log P o/w ), WLOGP (Wildman and Crippen log P o/w ), XLOGP3 (enhanced atomic/hybrid log P o/w 3), SILICOS-IT and MLOGP (quantitative-structure log P o/w ) were used to assess the lipophilicity of the ligands (Cheng et al., 2007).The consensus log P o/w which is the mean predicted lipophilicity values were obtained using these methods.Lower the consensus log P o/w value, the more soluble the compound (Ogata et al., 2018).The lipophilicity outcomes for the inhibitors D1, D2, D3, and D4 are 2. 55, 3.41, 3.02, and 2.70, respectively (Supplementary Table 3).All the values are in the lower range, indicating the four inhibitors have good solubility in a non-aqueous medium.Some drugs need to be highly water soluble to distribute enough of their active ingredients, which is further related to absorption by the body.There were 3 models of SwissADME in-house solubility predictor to predict the water solubility of the compounds namely ESOL (Estimated Solubility), Ali, and SILICOS-IT (Sluga et al., 2020).A qualitative estimation of solubility according to log S scale is: < À 10 (poorly soluble), < À 6 (moderately soluble), < À 4 (soluble), < À 2 (very soluble), and <0 (highly soluble) (Savjani et al., 2012).The water solubility of the four predicted inhibitors was analyzed by opting for the average of the above three parameters.D1 and D2 displayed moderate solubility in aqueous medium whereas D3 and D4 depicted poor solubility based on the predicted results (Supplementary Table 4).
As the drug is absorbed by the body, it encounters a variety of membrane barriers, including the target cell, hepatocyte membrane, gastrointestinal epithelial cells, blood capillary wall, and restrictive organ barriers (such as the blood-brain barrier).If the value of log K p (logarithmic skin permeation coefficient) is more negative, a molecule is said to be less skin permeant (Potts & Guy, 1992).The log K p values for D1, D2, D3, and D4, respectively are -6.09,-7.04, -6.31, and -5.99 cm/s.D4 exhibited maximumpermeability and D2 posed as the least skin permeable according to the predicted values.Moreover, additional parameters, such as gastrointestinal (GI) absorption or human intestinal absorption (HIA) were assessed to evaluate the bioavailability of the predicted inhibitors through a BOILED (Brain or IntestinaL EstimateD permeation) egg model (Daina & Zoete, 2016) (Supplementary Figure 9).All four inhibitors are found to have high GI absorption and HIA values indicating that are well absorbed in the gut system (Supplementary Table 5).No compound was found to have the tendency to cross the Blood Brain Barrier (BBB).This indicates that the compounds suggested here are relatively big in their overall chemical structure (Supplementary Figure 8) and cannot cross the blood-brain barrier.Additionally, a drug's less blood-brain permeability indicates, when metabolized, it will not enter or cause damage by entering the blood circulation to the brain cells (Pardridge, 2012).The remaining substances were expected to neither penetrate the brain nor be absorbed.The drugs that are administered into the patient's system, undergo absorption entering the systemic circulation, undergo metabolism, and finally are excreted out as nontoxic compounds.Additionally, it is critical to comprehend whether a particular compound functions as a substrate or an inhibitor of the P-glycoprotein (P-gp).This P-glycoprotein is a member of the class of transporters called ATP-binding cassette transporters required to assess active efflux through biological membranes (Goebel et al., 2021).It is also crucial to know about the process by which molecules interact with cytochrome P450(CYP) enzymes, which are involved in the removal of drugs via metabolic modification (Ogu & Maxa, 2000).According to a research study (Lin & Yamazaki, 2003), CYP and P-gp can process small molecules (drugs) in a collaborative manner.P-glycoprotein functions as an efflux pump, pumping the drug substrates from inside to outside the cells of the gut lumen as well as removing them from the kidney and liver.On the other hand, CYP enzymes are adapted as the body's main line of defense against xenobiotics (synthetic foreign molecules), and as part of this process, they also bioactivate drugs and toxins to produce more reactive intermediates (Lin & Yamazaki, 2003).Induction or inhibition of CYP enzymes is a major mechanism that underlies drug-drug interactions.Drug bioavailability is significantly affected by the natural defense system of an individual's body.Blocking these isoenzymes could cause side effects by lowering the drug's or its metabolites' solubility and rendering accumulation in the system.To comprehend the inhibitor's mechanism of action better, the top hit compounds were assessed for toxicity, deposition, and checked if they could both interact with P-gp and CYPs as a substrate or an inhibitor.All four inhibitors are found to be substrates of P-gp.For CYP1A2, D1 and D2 both act as inhibitors; for CYP2CI9, D2 and D3 act as inhibitors; for CYP2C9, D2, D3 and D4 act as inhibitor; for CYP2DC6, D1 and D2 act as inhibitor and lastly for CYP3A4, D2, D3 and D4 act as potential inhibitor (Supplementary Table 5).Overall, the oral bioavailability profile of the four inhibitors falls entirely into the bioavailability radar plot's pink area (Figure 25).This indicates that the inhibitors display a good pharmacodynamic and pharmacokinetic profile.

Inhibition constant of the four compounds
K i (inhibition coefficient) is used to denote the concentration of the inhibitor or the proposed drug that would lessen the activity of a target or reduce the rate of the reaction by half of the maximal rate (Bachmann & Lewis, 2005).It highlights the potency of an inhibitor in its action.K i is beneficial for specifying the likelihood that a specific drug will inhibit a specific enzyme and cause a clinically significant drug interaction with a substrate.The smaller the K i value of a drug against a definite receptor, the lesser the quantity of the drug required to restrict the specific activity of a particular target (Bachmann & Lewis, 2005).It further depicts a stronger binding affinity of the drug with a substrate for the enzyme to limit its functionality.This is because the lower K i means that the drug can occupy 50% of those receptors even when the drug is present in a lower concentration.The predicted K i values for the four proposed inhibitors (D1, D2, D3, and D4) are 1.34724 lM, 0.9609 lM, 0.578812 lM, and 0.1497 lM, respectively (Table 1).All the values obtained henceforth are in the lower range thus it showcases their high efficacy towards drug action in the biological system.Thus, these inhibitors can be administered in lesser quantities to the patient.As these inhibitors are all P-gp substrates, a lesser quantity of the inhibitors in the system makes them easier to be expelled.This prevents accumulation of the drug in the body and lesser chances of accumulation led toxic side-effects.

In silico toxicity evaluation of the top hit compounds
ADME properties' analysis is a primary and crucial step during the in-silico development of novel drugs.A drug should undergo additional toxicity tests to know if it is safe to be consumed by a patient.Here, an in-silico toxicity assessment of the top hit compounds was performed using Pro-Tox.LD50 of a drug is called the median lethal dose.This value signifies the amount of drug which when administered as a whole to the test subjects' body kills 50% of the test group (Raj et al., 2013).There are six different toxicity classes based on the LD50 values.On a scale from 1 to 6, 1 is the most toxic at the extreme left end and 6 is the least toxic towards right.The overall predicted toxicity class for D1 and D4 is 4 with obtained LD50 values of 650 and 400 mg/kg, respectively, for D2 it is 5 with a value of 10000 mg/kg and for D3 the class assigned is 6 with 3100 mg/kg as the obtained LD50 value, which indicates all the inhibitors have less to non-toxic nature.The Pro-Tox online server also predicts four toxicological endpoints, such as cytotoxicity, mutagenicity, carcinogenicity, and immunotoxicity.Hepatotoxicity is anticipated to determine whether the drug will result in liver dysfunction after consumption (David & Hamilton, 2010).If this attribute is found active in a drug, it might alter or hinder the normal liver functions of an individual.Hepatotoxicity is found to be active for D1 and D4.Immunotoxicity checks for substances that are known to alter the function of the immune system by inhibiting B cell proliferation (Phadnis-Moghe & Kaminski, 2017).D2 has minor immunotoxic outcomes.The active carcinogenicity might result in triggering cancerous proliferation of tissues and cells resulting in abnormal and uncontrolled growth (Williams & Weisburger, 1985), while the active mutagenicity indicates that it has the potent effect of causing mutagenesis and changing the information of genetic material like DNA (Barbour et al., 2006).Carcinogenicity and mutagenicity are found to be active for D3 and D4.Cytotoxicity which measures the toxic effect of drugs in the living cells is found to be negligible for all four inhibitors.Furthermore, the variety of enzymes present in a biological system could either allow the drug to be therapeutically metabolized or cause the production of toxic metabolites.This depicts that all the inhibitors have good pharmacodynamic and pharmacokinetic properties.The pink area represents the optimal range for each property (lipophilicity, size, polarity, TPSA, solubility: log S not higher than 6, and flexibility: no more than 9 rotatable bonds). in this case, the compounds are predicted to be orally bioavailable as they are within the radar, not being too flexible and too polar as predicted in the diagram.
domain-containing protein 5 (ATAD5).From the data in Table 2, it can be concluded that none of the targets are majorly affected negatively by the inhibitors predicted in the study.Only D1 mildly affects MMP.

Conclusion
The NS1 protein, encoded by the DENV2 strain, was detected to play a crucial role in the fatality of the disease and was involved in worsening the condition of the patients in the majority of the cases.Unveiling the enigmatic role of NS1 protein through numerous prior findings brings us to an understanding that the protein polymerizes through its various folds, enhancing its stability and secretion rate, corroborating a major participation in the immune evasion.Designing an inhibitor against Dengue CPE has been a major alarm because a potential therapeutic strategy is yet to be developed.A residual study of the NS1 protein was conducted which revealed that the cysteine residue at the 55th position was reckoned for its participation in intramolecular disulphide linkages, and the three terminal cysteine residues (at 313th, 316th, and 329th position) were in authority for the homodimer formation.However, between the two glycosylation sites, N130 was accountable for providing stability to the hexameric form of NS1, whereas, N207 facilitates secretion, provides extracellular protein stability, and assists cell surface attachment through interaction with C4BP.Thus, the objective of this investigation is to develop four new viable inhibitors-D1, D2, D3, and D4 for, restricting the sites responsible for the polymerization and cell surface attachment of the NS1.Following the pharmacophore-based virtual screening of the chemical library, a total of 4 ligand libraries (one for each cavity) were generated possessing a very strong affinity towards the target molecule.Molecular docking was carried out and the top ligand, exhibiting a strong interaction and stable ADMET profile, was chosen for the molecular dynamic' simulation run, to understand the extent up to which the ligand-induced stability to the protein.Molecular dynamics simulation was performed for each ligand towards their respective binding sites, which resulted in considerable favorable findings.Both RMSD and R g reveal that ligand binding has a noteworthy influence on the protein structure, resulting in lower deviations and minimal fluctuations throughout the 100 ns simulation period.The evaluation of the RMSF and of the overall SASA demonstrated a decrease in flexibility in addition to reduced atomic fluctuations, accompanied by an increase in the compactness of the protein, thereby leading to the attainment of correctly folded and stable conformations.The four complexes depicted the establishment of new correlated movements and a profound reduction in the overall unbound correlations, suggesting the enhancement in the stability of the complexes.A significant increment in the number of the non-covalent interactions-H-bonds and hydrophobic interactions, was suggestive of the fact that the association of the ligand with the protein established a higher degree of binding affinity between the two entities.All the ligand-bound complexes exhibited a less scattered and more compact representation of both energy surface and energy landscapes suggesting the ligand association minimized the overall energy of the protein.All four inhibitors were revealed to have a good ADMET profile and their lower K i values demonstrated their higher efficacy towards drug action.As a result, all four ligands have the potential to operate as a denguesuppressing medication and may be biologically examined in the laboratory for making it commercially obtainable to mankind.

Figure 1 .
Figure 1.The molecular mechanism approach is applied for drug designing in the research.(1) Infection spread through the bite of viral dengue mosquito vectors.(2) The genome organization of DENV2 consists of three structural proteins (capsid [C], membrane [M], and envelope [E]) and seven non-structural proteins (NS1, NS2 a/B, NS3, NS4 a/B, NS5).(3)NS1 protein with the target residues highlighted in purple (C55), yellow (N130), pink (N207) and green color (C313, C316, C329).(4), (5), (7) monomeric (4), dimeric (5), and hexameric forms (7).NS1 is first expressed in infected cells as a monomer.It translocates co-transcriptionally to the ER lumen.(6) Further oligomerization leads to the formation of NS1 hexameric structure formed with glycosylated combinations with lipids.(4a, 4b) D3, blocking the C55 residue position that is involved in the formation of disulphide linkages, thus in turn stopping the CPE.(5a) D1 will block C313, C316, and C329 that is majorly involved in homodimer formation.The dimers further modify cellular lipid dynamics during CPE and/or lead to the recruitment of cholesterol and triglycerides to the viral-replication complex.(7a) Through interactions with the plasma complement regulator C4BP, NS1 binds to the cell surface through N207.D4 would be targeting and blocking the active glycosylation site at N207. (7b) N130 stimulates the oligomeric assembly of secreted NS1 along with the stabilization of the secreted hexamer.D3 would be targeting and blocking the N130 site to disrupt the hexametric stability of the NS1 protein vital to viral stability, contagiousness, and antigenicity.

Figure 2 .
Figure 2. The overall methodology applied in the study.

Figure 3 .
Figure 3.The pharmacophore features of (a) Cav7, (b) Cav10, (c) Cav14, and (d) Cav15 are extracted and depicted in This figure.The colored mesh-like spheres exhibit different features of the pharmacophore.The red sphere represents a negative ion center, while the white spheres represent the presence of a hydrogen donor, the yellow spheres denote hydrogen acceptors, and the green sphere symbolizes the hydrophobic group.

Figure 4 .
Figure 4. RMSD curves of the (a) D1Complex, (b) D2Complex, (c) D3Complex, and (d) D4Complex, superimposed against the RMSD curve of the unbound NS1, extracted from the 100 ns simulation run.The orange line in all four graphs symbolizes the fluctuation of the unbound protein, however, the curves superimposed against the orange curve represented by yellow, green, brown, and blue depicts the fluctuations for D1Complex, D2Complex, D3Complex, and D4Complex, respectively.

Figure 5 .
Figure 5. R g of the (a) D1Complex, (b) D2Complex, (c) D3Complex, and (d) D4Complex, superimposed against the R g of the unbound NS1 protein.For each graph the fluctuation exhibited by the unbound protein is displayed by the orange line as against the orange line in all four graphs symbolizes the fluctuation of the unbound protein, however, the curves superimposed against the orange curve represented by yellow, green, brown, and blue depicts the fluctuations for D1Complex, D2Complex, D3Complex, and D4Complex, respectively.

Figure 6 .
Figure 6.RMSF graphs generated for the unbound protein and the (a) D1Complex, (b) D2Complex, (c) D3Complex, and (d) D4Complex, superimposed, have been depicted.The orange line in all four graphs symbolizes the fluctuation of the unbound protein, however, the curves superimposed against the orange curve, represented by yellow, green, brown, and blue depicts the fluctuations for D1Complex, D2Complex, D3Complex, and D4Complex, respectively.

Figure 7 .
Figure 7. Solvent accessible surface area (SASA) analysis of the unbound protein and the ligand-bound complexes are displayed.The curve for unbound_NS1, represented by orange line, is superimposed over the curves of (a) D1Complex, (b) D2Complex, (c) D3Complex, and (d) D4Complex, represented by yellow, green, brown, and blue, respectively.

Figure 8 .
Figure 8. Graphs for the number of hydrogen bonds formed during the 100 ns simulation has been depicted.The graphs for (a) D1Complex, (b) D2Complex, (c) D3Complex, and (d) D4Complex have been superimposed on the curve for the unbound_NS1.The orange line for all four graphs represents the NS1 protein in its unbound state, whereas, the yellow, green, brown, and blue lines represent the D1Complex, D2Complex, D3Complex, and D4Complex, respectively.

Figure 9 .
Figure 9.The residual correlative motion in protein structures was visualized by the dynamic cross-correlation Map analysis, for (a) unbound_NS1, (b) D1Complex, (c) D2Complex, (d) D3Complex, and (e) D4Complex.In all cases, correlations are represented in the color spectrum, from cyan to pink.Parallel motion of the residues indicates correlated movement, however, antiparallel motion of the residues exhibits anticorrelated movement.The positive correlation between residues, that is, residues move in the same direction, is represented as cyan color.Alternatively, residues having opposite motion are denoted as the anti-correlative motion marked as pink color.

Figure 10 .
Figure 10.2D-PCA of unbound protein and the D1Complex, superimposed, over the 100 ns simulation.The orange dots depict the unbound state of the protein, whereas, the yellow dots represent the protein in its D1 bound state.

Figure 11 .
Figure 11.Superimposed depiction of the 2D-PCA of the unbound protein and the D2Complex over the 100 ns simulation.The orange dots indicate the unbound protein, whereas, the green dots represent the D2Complex.

Figure 12 .
Figure 12.Over a 100 ns simulation, 2D-PCA of unbound protein and the D3Complex were superimposed.The orange dots represent the protein in its unbound form, whereas the brown dots indicate the protein in its D3 bound state.

Figure 13 .
Figure 13.Superimposed depiction of the 2D-PCA of the unbound protein and the D4Complex over the 100 ns simulation.The unbound protein is exhibited by the orange dots, whereas the D4Complex is represented by the blue dots.

Figure 14 .
Figure 14.(a) Total energy landscape and (b) total energy surface of the unbound_NS1 protein.The minima depict the lowest energy unbound conformation during the simulation.

Figure 15 .
Figure 15.(a) Total energy landscape and (b) total energy surface of the D1Complex.The minima depict the lowest energy D1Complex conformation during the simulation.

Figure 16 .
Figure 16.(a) Total energy landscape and (b) total energy surface of the D2Complex.The minima depict the lowest energy D2Complex conformation during the simulation.Figure 17.(a) Total energy landscape and (b) total energy surface of the D3Complex.The minima depict the lowest energy D3Complex conformation during the simulation.

Figure 17 .
Figure 16.(a) Total energy landscape and (b) total energy surface of the D2Complex.The minima depict the lowest energy D2Complex conformation during the simulation.Figure 17.(a) Total energy landscape and (b) total energy surface of the D3Complex.The minima depict the lowest energy D3Complex conformation during the simulation.

Figure 18 .
Figure 18.(a) Total energy landscape and (b) total energy surface of the D4Complex.The minima depict the lowest energy D4Complex conformation during the simulation.

Figure 19 .
Figure 19.Total short range interaction energy analysis of (a) D1Complex, (b) D2Complex, (c) D3Complex, and (d) D4Complex.These graphs represent the total sum of the Lennard Jones energy and the Coulombic interactions in kJ/mol, generated throughout the 100 ns simulation.

Figure 20 .
Figure 20.Binding affinity of the first 20 minima conformations of the four complexes to estimate the average binding affinity throughout the simulation.

Figure 21 .
Figure 21.(a) Illustration of the interaction pattern analysis of the D1Complex.The red arcs represent the protein interacting residues and the red dashes emerging from them indicate hydrophobic interactions with the ligand.Ligand is exhibited in purple, and the green dashes represent hydrogen bonds with the ligand.Carbon atoms are shown in black, whereas nitrogen, oxygen, and fluorine atoms are depicted in blue, red, and light green colors, respectively.(b) The binding pose of the D1 ligand.The pink region denotes H-bond donor surface and green region denotes H-bond acceptor surface.

Figure 22 .
Figure 22.(a) an Illustration of the D2Complex interaction pattern analysis.The red arcs show protein-interacting residues, and the red dashes that emerge from them suggest hydrophobic interactions with the ligand.The ligand is shown in purple.Carbon atoms are portrayed in black, whereas nitrogen and oxygen atoms are depicted in blue and red, respectively.A fluorine atom, evidently present in the ligand, is denoted by light green color.(b) The binding pose of the D2 ligand.Pink region denotes H-bond donor surface and green region denotes H-bond acceptor surface.

Figure 23 .
Figure 23.(a) Illustration of the D3Complex interaction pattern analysis.The red arcs show interacting residues of the protein, and the red dashes emerging from them suggest hydrophobic interactions with the ligand.The ligand is shown in purple, while the green dashes depict hydrogen bonding with it.Carbon atoms are portrayed in black, whereas nitrogen and oxygen atoms are depicted in blue and red, respectively.Sulphur atoms are represented in yellow color.(b) The binding pose of the D3 ligand.Pink region denotes H-bond donor surface and green region denotes H-bond acceptor surface.

Figure 24 .
Figure 24.(a) Representation of the D4Complex interaction pattern analysis.The red arcs show interacting protein residues, and the red dashes that emerge from them suggest hydrophobic interactions with the ligand.The ligand is shown in purple; however, the green dashes denote hydrogen bond formed between the two entities.Carbon atoms are illustrated in black, whereas nitrogen and oxygen atoms are depicted in blue and red, respectively.(b) The binding pose of the D4 ligand.Pink region denotes H-bond donor surface and green region denotes H-bond acceptor surface.
Figure 25.Schematic diagram of bioavailability radar: bioavailability radar of the top four inhibitors based on their suitable physicochemical indices lie under the pink colored region.This depicts that all the inhibitors have good pharmacodynamic and pharmacokinetic properties.The pink area represents the optimal range for each property (lipophilicity, size, polarity, TPSA, solubility: log S not higher than 6, and flexibility: no more than 9 rotatable bonds). in this case, the compounds are predicted to be orally bioavailable as they are within the radar, not being too flexible and too polar as predicted in the diagram.

Table 1 .
Docking based binding energy, average binding affinity throughout the simulation and inhibition constant of the top potential ligands of the target protein 'NS1' of dengue virus.