In-silico evidences on filarial cystatin as a putative ligand of human TLR4

Abstract Cystatin is a small molecular weight immunomodulatory protein of filarial parasite that plays a pivotal role in downregulating the host immune response to prolong the survival of the parasite inside the host body. Hitherto, this protein is familiar as an inhibitor of human proteases. However, growing evidences on the role of cystatin in regulating inflammatory homeostasis prompted us to investigate the molecular reasons behind the explicit anti-inflammatory trait of this protein. We have explored molecular docking and molecular dynamics simulation approaches to explore the interaction of cystatin of Wuchereria bancrofti (causative parasite of human filariasis) with human Toll-like receptors (TLRs). TLRs are the most crucial component of frontline host defence against pathogenic infections including filarial infection. Our in-silico data clearly revealed that cystatin strongly interacts with the extracellular domain of TLR4 (binding energy=-93.5 ± 10 kJ/mol) and this biophysical interaction is mediated by hydrogen bonding and hydrophobic interaction. Molecular dynamics simulation analysis revealed excellent stability of the cystatin-TLR4 complex. Taken together, our data indicated that cystatin appears to be a ligand of TLR4 and we hypothesize that cystatin-TLR4 interaction most likely to play a key role in activating the alternative activation pathways to establish an anti-inflammatory milieu. Thus, the study provokes the development of chemotherapeutics and/or vaccines for targeting the cystatin-TLR4 interaction to disrupt the pathological attributes of human lymphatic filariasis. Our findings are expected to provide a novel dimension to the existing knowledge on filarial immunopathogenesis and it will encourage the scientific communities for experimental validation of the present investigation. Communicated by Ramaswamy H. Sarma


Introduction
Lymphatic filariasis (LF) is a debilitating mosquito-borne parasitic disease of humans and it is caused due to the infection of pathogenic filarial nematodes Wuchereia bancrofti, Brugia malayi and B. timori. W. bancrofti is considered is the principal causative parasite as it causes 90% of the total cases of human filariasis worldwide (World Health Organization (WHO), 2020). Currently the epidemiological map constitutes over 120 million infected people in 72 countries across the globe while 893 million people are at risk of infection (World Health Organization (WHO), 2020). India alone contributes 40% of the total diseases burden and the disease is spread in 17 states and 6 union territories (Riches et al., 2020). Disruption of inflammatory homeostasis and dysfunction of the lymphatic system including lymph stasis, lymphedema, lymphangitis, and elephantiasis are major clinical manifestation of LF (Mukherjee, Huda, et al., 2019). LF patients are also prone to secondary microbial infections caused by bacteria, fungus, and protozoa resulting in the impairment of important organ functions (Mukherjee et al., 2014). Albendazole (a broad-spectrum anthelmintic), diethylcarbamazine (a parasitic metabolic inhibitor), and ivermectin (a Clion channel blocker) are used for treating filarial infection (Mukherjee et al., 2018). However, these drugs possess several constraints like limited efficacy against adultstage parasite, emergence of resistance, low absorption, and toxic side effects (Mukherjee et al., 2018). Global control program strategy adopted by World Health Organisation (WHO) has not reached success due to the lack of complete knowledge in understanding the pathogenesis of the filarial nematode, especially W. bancrofti (Mukherjee et al., 2018). In general, the infective stage of the filarial parasite (L3) enters the human body via mosquito bites and subsequently matures to the adult-stage (L5) through several developmental changes. Adults parasites produce larvae (microfilariae or L1) and are transmitted to mosquito to complete the cycle. Adult filarial parasite usually stays for several years (5-10 years) and the major determinant of their survival inside the human body is the explicit efficiency of the parasite to modulate the immune response of human (Mukherjee, Joardar, et al., 2019;. In this context, filarid exploits a number of somatic (cystatin, chitohexose) and excretory-secretory proteins (ES-62)/glycoproteins/glycans that trigger suppression of anti-parasitic immunity through various ways (Babu & Nutman, 2014;Bisht et al., 2019;Hoerauf et al., 2005).
Filarial Cystatin, an immunomodulatory glycoprotein having a molecular weight of $15 kDa protein has been characterized as a strong reversible competitive inhibitor of human cysteine proteases viz. cathepsins B, L, and S (Gregory & Maizels, 2008). A total of three cystatin homologs namely Bm-CPI-1, Bm-CPI-2, and Bm-CPI-3 have been reported in B. malayi (Bisht et al., 2019) while a single gene coding for cystatin has been documented in W. Bancroft (Small et al., 2019). Cystatins of parasitic nematodes are well-described pathogenicity factors that contribute to the downregulation of T-cell proliferation of their hosts and induce anti-inflammatory cytokine responses (Khatri et al., 2020). Considering intense immunomodulatory property, it was hypothesized that the functional evolution of cystatin could have played a major role in the parasitic adaptation of the filarial parasites (Schierack et al., 2003). Unlike other immunomodulatory molecules, filarial cystatins trigger multiple mechanisms to subdue the detrimental proinflammatory responses induced by the host immune system. Filarial cystatins have been reported to suppress the expansion of human peripheral blood mononuclear cells (PBMC) and murine spleen cells. Interestingly, filarial cystatin can significantly upregulate the expression and secretion of interleukin (IL)-10 from human PBMC (Schierack et al., 2003) and inducible nitric oxide by IFN-gamma-stimulated mouse macrophages (Hartmann et al., 2002). The upregulation of IL-10 and induction of nitric oxide in macrophages by cystatin clearly indicated towards the involvement of alternative activation of macrophages which governs the suppression of host immaune response and creates a parasite-friendly environment (Gordon & Martinez, 2010). Intriguingly, alternative activation of macrophages is known to be induced as well as regulated by the Toll-like receptors (TLRs), mostly by cell surface TLRs viz TLR1, TLR2, TLR4, TLR5, and TLR6 (Gordon & Martinez, 2010;Mukherjee et al., 2016;Panda et al., 2012). Therefore, the possible role of cystatin in the regulation of TLR activity is an open question to answer. The present study is a maiden report demonstrating W. bancrofti as a putative ligand of human TLR4 through in-silico approaches.

Retrieval of sequence
The amino acid sequence of the cystatin of W. bancrofti (Accession Id: EJW82673.1) was retrieved in FASTA format using the tool, available in NCBI (National Center for Biotechnology Information) (https://www.ncbi.nlm.nih.gov/) for further in-silico analysis. The gene sequence was also deduced from the amino acid sequence using CDS (Protein coding sequence) search given in NCBI.

Construction of phylogenetic tree
Cystatin sequences from different filarial parasites including W. bancrofti and other organisms were subjected to multiple sequence alignment (MSA) employing MEGA-X application and an unrooted phylogenetic tree was constructed using maximum parsimony method.

Physiochemical characterization
The physico-biochemical properties of cystatin were explored through different online servers. ProtParam server provided in Expasy database (https://web.expasy.org/protparam/) was accessed to analyze different physiochemical properties like molecular weight (MW), theoretical isoelectric point (pI), number of positive and negative charged residue, extinction coefficient (EC), estimated half-life, instability index (II), aliphatic index (AI) and grand average of hydropathicity (GRAVY) while solubility was predicted from ProteinSol web server (https://protein-sol.manchester.ac.uk/).

Evaluation of antigenicity and allergenicity
Antigenicity was predicted by VaxiJen v2.0 server and later was validated using AntigenPro server. Where AllerTop v2.0 and AllergenFP v1.0 both servers were used to evaluate the allergenic property of the cystatin.

Molecular docking and analysis of protein-protein interaction
In order to examine the candidature of cystatin as a possible ligand of human TLRs, a molecular docking study was performed. A molecular docking experiment was conducted using the predicted structure of W. bancrofti cystatin and crystal structure of human TLR4-MD2 complex (PDB ID: 3FXI_AC), TLR1-TLR2 complex (PDB ID: 2Z7X_AB), and TLR5 (PDB ID: 3J0A) were obtained from RCSB protein data bank (https://www.rcsb.org/). Cystatin was separately docked with each target TLR using HEX 8.0.0 software and ClusPro web server (https://cluspro.bu.edu/login.php?redir=/queue.php) (Singh et al., 2019). PyMOL was used to visualize the docked complex and interactions were analyzed using Discovery Studio 2020 Client software.

Molecular dynamics (MD) simulations
MD simulations were carried out for the human TLR4 with cystatin native and 5 mutants (Asn52 to Ser52, Arg87 to His87, Lys88 to Glu88, Lys97 to Glu97, and Asp101 to Gly101) for a period of 50 ns using GROMACS v5.1 molecular dynamics package Pronk et al., 2013). Hþþ server (http://biophysics.cs.vt.edu/index.php) was accessed first to generate protonation states of the cystatin residues at pH 7 and the resulted cystatin was then employed for MD simulation (Kumari et al., 2020).
The unit cell defined as a cubical box, with a minimal distance of 10 Å from the protein surface to the edges of the box, was solvated using the Simple Point Charge (SPC) water model GROMOS96 53a6 force field (Oostenbrink et al., 2004). Counter-ions were added to make every system electrically neutral at a salt concentration of 0.15 mol/ L. Before the MD run, each system was subjected to energy minimization by employing the steepest descent integrator for 5000 steps with force convergence of <1000 kcal mol À1 nm À1 .
Thereafter, the protein-ligand complex was equilibrated for 5 ns using the canonical (NVT) and isothermal-isobaric (NPT) ensembles. During equilibration, the system was coupled with the Berendsen temperature and the Parrinello-Rahman pressure controllers, respectively, to maintain a temperature of 300 K and pressure 1 bar. The Particle Mesh Ewald (PME) algorithm (Essmann et al., 1995) was employed to deal with the long-range Coulomb interaction with a Fourier grid spacing of 0.12 nm. The short-range van der Waals interaction was given by the Lennard-Jones potential with a cut-off distance of 1 nm. All bond lengths were constrained by the linear constraint solver (LINCS) method (Hess et al., 1997).
Subsequently, 50 ns was considered for production under the micro-canonical ensemble by relaxing the couplings with the thermostats. A time step of 2 fs was used and the coordinates were saved at every 10 ps during the production run. For the structural analyses of every system, the resultant MD trajectories saved at the interval of 100 ps were analyzed using the built-in modules of GROMACS and visual molecular dynamics (VMD 1.9.1) (Schuler et al., 2001). The 2 D plots depicting the intrinsic dynamical stabilities captured by the root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), solvent accessible surface area (SASA), and radius of gyration (R g ) of the complexes were generated by Grace 5.1.23 program.

Principal component analysis (PCA)
Essential dynamics (ED) or principal component analysis (PCA) is an efficient statistical method applied to reduce the number of dimensions needed to describe protein molecular dynamics through a decomposition process (Humphrey et al., 1996). The GROMACS module tool called g_covar was used to yield the eigenvalues and eigenvectors by calculating and diagonalizing the covariance matrix, whereas the g_anaeig tool was used to analyze and plot the eigenvectors Swain et al., 2018).

Calculation of binding and per-residue interaction energy
Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/ PBSA) is widely used for free energy calculation from last 5 ns (45-50 ns) MD trajectory (Aalten et al., 1995;Malik et al., 2019;Wang et al., 2019). The binding free energy (DG bind ) in a solvent medium was calculated as follows: Where G complex is the total free energy of the substrate-protein complex, G protein and G ligand are the total energies of protein and substrate alone in a solvent, respectively. The free energies for each individual G complex , G protein and G ligand were estimated by: Where p can be protein, ligand, or complex. E MM is the average molecular mechanics potential energy in the vacuum and G solv is the solvation free energy.
The molecular mechanics potential energy was calculated in the vacuum as follows: Where E bonded or (E int ) is the total bonded interaction, which includes bond angle, dihedral, and improper interactions; E non-bonded is the total non-bonded interaction consisting of both van der Waals (E vdw ) and electrostatic (E elec ) interactions. E bonded is always taken as zero.
The solvation free energy (G solv ) was estimated as the sum of electrostatic solvation free energy (G polar ) and nonpolar solvation free energy (G non-polar ) as given below: Where G polar , the polar solvation energy, was determined using the Poisson-Boltzmann (PB) linear equation and the nonpolar contribution, G non-polar , was estimated from the solvent accessible surface area (SASA) as per the following equation: Where c (coefficient related to surface tension of the solvent) ¼ 0.02267 kJ/mol/Å 2 and b ¼ 3.849 kJ/mol.
The binding free energies of the complex were calculated based on 5000 snapshots taken at an equal interval of time from 50 ns MD simulations. The per-residue energy contribution was also computed to understand the contribution of individual amino acids to the total binding energy.

Assessment of conformational changes in TLR4 after binding of cystatin
Conformational changes in TLR4 after binding of cystatin were examined by superimposing the structure of unbound and bound TLRs obtained from molecular docking and MDS following Liou et al., 2014. In addition, a comparative study on TLR4-cystatin complexes before and after MD simulation was also performed to analyze the conformational variations of important residues required for protein-protein binding Dhankhar, Dalal, Mahto, et al., 2020).

Determination of the relative importance of different amino acids
In addition to the MD simulation study, the relative contribution of the key amino acids present in the binding surface and their role in mediating the binding of cystatin to TLR4 was investigated through the induction of mutation in-silico. Different mutant cystatin structures were separately examined for their binding to TLR4 and binding efficiency was investigated.

In-silico analyses of phylogeny and physico-biochemical properties of W. bancrofti cystatin
Cystatin of W. bancrofti contains 127 amino acids out of which frequency of isoleucine (8.7%) and lysine (9.4%) are copious (Figure 1(a and b)). Cystatin possesses an isoelectric point (pI) of more than 7 and a lower instability index that collectively indicated towards the stability of the protein ( Table S1). The solubility of the cystatin was then achieved from ProteinSol webserver that was about 0.591 (Table S1). The predicted localization studies have suggested W. bancrofti cystatin as a membrane-bound extracellular (secreted) protein containing membrane specific signal peptide and non-cytoplasmic domain. In addition, this protein also contains a single 23 amino acids long (MFFPIVWFSVLLIISKSFAREIR) transmembrane domain. Filarial cystatin was found to contain three protein kinase C phosphorylation sites, a casein kinase II phosphorylation site, and prenyl group binding site which indicated phosphorylation and prenylation as the principal post-translational modifications of this protein (Table S2). The predicted antigenic score of the cystatin was found to be in the range of 0.7627-0.808068 while analyses using AllerTop v2.0 and AllergenFP v1.0 indicated cystatin as a non-allergen.
We have studied the microevolution of cystatin protein across the different filarial nematodes and humans ( Figure  1(c)). An unrooted phylogenetic tree generated from the phylogenetic reconstruction resembles B. malayi cystatin as the closest homolog of W. bancrofti cystatin and these two homologs share a common monophyletic group supported by a high bootstrap value (BV) value of 96. Ascaris lumbricoides was inferred as the most distantly related nematode inferred through the maximum parsimony tree (Figure 1(c)). Considering the position of humans in the aforesaid tree, alignment of the amino acid sequences of W. bancrofti cystatin was compared with human and 24.4% similarity was observed (Figure 1(d)). The proportion of alpha-helix, strand, beta-turn, and coil in filarial cystatin is 42.52%, 20.47% 2.36%, and 34.65% respectively while in cases of human it is 27.48%, 15.27% 3.05%, and 54.2% ( Figure S1, Figure S2A and S2B, Table S3). We have also deduced the gene sequence of W. bancrofti cystatin from the amino acid sequence and again compared with the human cystatin gene (Figure 1(e and f)). Sequence alignment study revealed 45% gene homology between human and filarial cystatin (Figure 1(f)).
Homology modeling, post-modeling analyses, and in-silico characterization Homology modeling generated the predicted 3 D structure of W. bancrofti cystatin and the structure is presented both in space-filling and cartoon models (Figure 2(a)). Stereochemical quality and energy parameters of this modelled structure were respectively evaluated using ProSA as well QMEAN4 web servers (Figure 2(b) and Figure S2C and S2D). A Z-score value of À4.54 in ProSA and QMEAN4 score of À1.29 collectively indicated the energetic and stereochemical quality of modelled structure of cystatin (Figure 2(b) and Figure S2C and S2D). Moreover, a comparison study on Z-score reflected the relatedness of the structure quality between W. bancrofti modelled cystatin with the experimentally developed structure of human cystatin ( Figure S4A). In addition, the stereochemical quality of filarial cystatin structure was also confirmed in the Ramachandran plot (Figure 2(c)). Ramachandran plot depicted that 96.19% of the residues are in the favourable region (shown as red and yellow) which resembles good model quality ( Figure  2(c)). Moreover, an ERRAT value of 87.058 and further evaluation through Verify 3D (determined using SAVES server) validated the conformational stability (Figure 2(d) and S4B). In addition, structural superposition of the modelled cystatin of W. bancrofti with human cystatin has 0.121 Å RMSD for 415 atoms. In order to finalize the model structure, MD simulation was employed and the RMSD plot reflects a very satisfying report ( Figure S5). After confirming the model quality, the developed structure was used for further experiments.

Molecular docking study and analysis of the binding of cystatin to human TLRs
The present study seeks to determine the role of cystatin in filarial immunomodulation. A series of previous studies in this direction put forwarded the postulation that cystatin may trigger the function of human cell surface TLRs to establish a parasite-friendly environment inside the host (Bisht et al., 2019;Gregory & Maizels, 2008;Khatri et al., 2020;Schierack et al., 2003). We have investigated this postulation by determining the interaction of filarial cystatin with the native structure of each of the cell surface TLRs through molecular docking ( Figure S3). Interactions amongst filarial cystatin and human cell surface TLRs viz. TLR1-TLR2 dimer, TLR4-MD2 complex, TLR5, and TLR2-TLR6 complex are presented in Table 1. The order of strength of binding between cystatin and TLRs is cystatin-TLR4 > cystain-TLR5 > cystatin-TLR6 > cystatin-TLR1-TLR2. The comparative interaction studies revealed the interaction between filarial cystatin and TLR4-MD2 complex as the strongest interaction with docking/interaction energy (E TOTAL ) of À634.9 (Table 1) and the predicted results were validated by the weighted scores determined from ClusPro web server. Therefore, further studies on the protein-protein interactions were emphasized on the cystatin-TLR4 complex.

Analysis of the binding of cystatin to human TLR4
Cystatin-TLR4 interaction was explored for protein-protein interaction using PyMOL for visualizing the docked complex in spherical appearance which evidences the binding of cystatin specifically at the extracellular domain of TLR4 ( Figure  3(a and b)). Such binding of cystatin at the extracellular domain of native TLR4-MD2 complex clearly indicated cystatin as a putative ligand of TLR4 and this inference is supported by other studies conducted on demonstrating TLR4 ligand from filarial parasite (Mukherjee et al., 2017;Panda et al., 2012). The major interacting forces involved in the interaction between the two proteins were determined and visualized using Discovery Studio 2020 Client software (Figure 3(c)). We have found that the binding interface between cystatin and TLR4 is consisting of hydrogen bonding and hydrophobic interactions (Table 2). Lys97 and Arg87 residues were predominantly present at the binding surface of cystatin and involved in forming hydrophobic interactions with the amino acids namely Tyr403, His426, and Ile50 present at the extracellular domain of TLR4 (Figure 3(c) and Table 2). On the other hand, Asn52, Arg87, Lys88, Lys97, and Asp101 of cystatin were found to form hydrogen bonds with Gln73, Asp405, Asp428, Glu474, and Arg355 respectively ( Table 2).
We have examined the relative importance of the aforesaid amino acids in cystatin, especially at the level of ligandreceptor interaction. We have introduced five random mutations in the amino acids (Asn52 to Ser52, Arg87 to His87, Lys88 to Glu88, Lys97 to Glu97, and Asp101 to Gly101) present in the interacting surface of cystatin and mediating the binding of cystatin to TLR4. Each mutant cystatin was docked with TLR4 separately and interaction energies were compared with the wild type (Figure 3(d)). We have found that the change of amino acid no. Lys88 to Glu88 in cystatin resulted in a significant weakening of the binding of cystatin to TLR4 (Figure 3(d) and Table 3).

Biophysical basis of cystatin-TLR4 interaction
Taking clues from the molecular docking, cystatin has been found to be a strong interacting partner of human TLR4. In order to study the biophysical basis of the interaction, normal mode analysis (NMA) and molecular dynamic (MD) simulation were studied for the cystatin-TLR4 complex.      NMA study NMA analysis was performed to investigate the biophysical stability and alterations in the cystatin-TLR4 complex (Figure 4). The domain motion in the protein complex is represented as a curved arrow (Figure 4(a)). High deformability (shown as hinges) throughout the chain of cystatin-TLR4 complex was observed and this is due to high helical content (Figure 4(b)). B-factor values representing the relative amplitude of atomic displacements also revealed higher deformability in the protein complex which indeed indicates a greater flexibility of the complex (Figure 4(c)). Eigenvalue indicates the impact of each deformation movement in the total protein motion (L opez- Blanco et al., 2014). A greater eigenvalue represents localized displacement while a low eigenvalue corresponds to collective conformational changes in the protein complex (L opez- Blanco et al., 2014). Herein, a low eigenvalue of 3.086933e À05 that indicated lower energy requirement in deforming the structure which in turn reveals a greater stability of the complex and an overall conformational change in the complex (Figure 4(d)). This inference is supported by a previous report (Abdelli et al., 2020). In NMA-based determination of molecular flexibility, the amplitude of the fluctuation of each mode varies inversely to the eigenvalue and it is quantitatively expressed as variance (Abdelli et al., 2020). Percent variance for each normal mode is represented by red colours and cumulative variance by green (Figure 4(e)). 80% of the variance is achieved by the first 10 modes out of the total 20 modes (Figure 4(e)). The covariance matrix is presented through the graphical representation where red represents correlated motion between a pair of residues, white colour for uncorrelated motion, and blue for anticorrelated motion (Figure 4(f)). Elastic network map graph represents the pairs of atoms connected by springs and each grey dot represents one spring (Figure 4(g)). In this study, an abundance of grey dots in the elastic network map reveals the motion stiffness in the protein complex (Figure 4(g)).

Molecular dynamics trajectory analysis of cystatin-TLR4 complex
MDS provides great advantages to analyze the internal motions, conformational changes, stability, etc. of protein-ligand complexes (Sen Gupta, Biswal, et al., 2020;Singh et al., 2020). Using the MDS trajectory, the root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (R g ), solvent accessible surface area (SASA), principal components (PC), and binding free-energy of the cystatin-TLR4 complex were determined to explore the structural stability, binding modes, and binding strength ( Figure 5).

Conformational stability of the cystatin-TLR4 complex
RMSD is a useful parameter to investigate the stability of a protein in a complex (Sen Gupta, Biswal, et al., 2020;Singh et al., 2020). A larger RMSD value is indicative of a lower stability of a protein complex. The RMSD plot for cystatin-TLR4 complex against the simulation time displays a minute deviation and less oscillation with an average value of 0.52 nm that indicate the attainment of a stable conformation, strong binding, and large stability of the complex throughout the whole dynamics ( Figure 5(a)). We have also compared between the RMSD profiles of cystatin-TLR4 complex and free TLR4 (Figure 5(b)). Our data indicate the average RMSD value of the former as 0.52 nm, which is smaller than that of the later (0.56 nm), suggesting better binding as well as the stability of the cystatin-TLR4 complex ( Figure 5(b)). After the initial rise of RMSD till 6 ns, both cystatin-TLR4 complex and free TLR4 attain the convergence and remain stable with very little fluctuation throughout the 50 ns MD simulation ( Figure 5(b)). RMSD results of mutant cystatin (Arg87 to His87, and Asp101 to Gly101)-TLR4 complexes showed higher RMSD as compare to native cystain-TLR4 complex, suggesting these (Arg87, and Asp101) amino acids are very essential in forming cystatin-TLR4 complex ( Figure S7A).

Residue flexibility analysis
RMSF is a dynamic parameter that measures the average main chain flexibility at each residue position. Analysis of the flexibility of the cystatin-TLR4 complex in terms of fluctuations of each amino acid as a function of residue number is presented by an RMSF plots ( Figure 5(c)). A large RMSF value indicates the flexibility of the protein structure, loose bonding, or the presence of loops. By contrast, a small RMSF value indicates the stability and the presence of the secondary structures like sheets and helices. While analyzing the TLR4cystatin complex, similar to the trend of RMSD, overall a smaller fluctuation as well as a mean RMSF value of 0.33 nm demonstrate greater stability of the complex (Figure 5(c)). Some terminal residues of TLR4 in the cystatin complex lead to a very negligible fluctuation of RMSF ( Figure 5(c)). RMSF plot disclosed 4 mutant cystatins (Asn52 to Ser52, Arg87 to His87, Lys88 to Glu88, and Asp101 to Gly101) -TLR4 complexes among the 5 possess high RMSF value than native cystatin-TLR4 complex ( Figure S7B).

Compactness analysis
The radius of gyration (R g ) describes the level of compaction of protein. It is defined as the mass-weighted root-meansquare distance for a collection of atoms from their common center of mass. Hence, the trajectory analysis for R g provides conformational variation, compactness, and tertiary structural volume of the protein-ligand complex. Lower average R g value of 3.14 nm of the TLR4-cystatin complex akin to RMSD and RMSF indicates very compact and stable complex formation between cystatin with TLR4, necessarily arising due to a very strong interaction between them ( Figure 5(e)). While comparing between the cystatin-TLR4 complex and free TLR4, the average value of the former is lower than that of the latter, respectively, having 3.14 nm and 3.19 nm, indicative of the greater stability of the complex than the free TLR4 ( Figure 5(f)). R g plot of mutant cystatin-TLR4 complex indicates Arg87, Lys88, Lys97, and Asp101 are highly responsible for overall compactness of cystatin ( Figure S7C).

Solvent accessible surface area
Solvent Accessible Surface Area (SASA) accounts for the surface area of a protein-ligand complex that directly interacts with the solvent molecules. The increase in SASA denotes structural expansion and more exposure to the solvent. Herein, the average value of SASA for the cystatin-TLR4 complex is found to be 238 nm 2 with no major fluctuation throughout the 50 ns MD simulation implying stability of the complex (Figure 5(g)). SASA plots of mutant cystatin-TLR4 complexes displayed that all the mutant complexes possess higher SASA value than native one ( Figure S7D).

Principal component analysis (PCA)
For principal component analysis (PCA), essential dynamics (ED) was employed to extract the directions of principal motions by a set of eigenvectors (EV) called principal or essential modes. PCA is a linear transformation for obtaining the most important spatial data using a covariance matrix of the Cartesian coordinates that describe the accessible degrees of freedom (DOF) of protein in the MD trajectory. This covariance matrix is diagonalized to extract a set of eigenvectors and eigenvalues that reflect a correlated, concert motion of the molecule (Khan et al., 2016). Each MD trajectory was projected onto the phase space to yield a spectrum of EVs, which depicted the vectorial representation of each single component in motion. Each EV holds an eigenvalue that describes the energetic contribution of each component to the motion. The PCA results are represented in Figure 6(a-d).
It was observed that 90-95% of the backbone motion was covered by the first 30 EVs where an exponentially decaying curve of eigenvalues was obtained against the EVs ( Figure  6(a-d)). The differential scattering of atoms in contour-like plots specifies the occurrence of conformational changes in the complex in agreement with the other MD analyses.
Taken the first two PCs into consideration, simulation results reveal a higher subspace dimension for the Free TLR4 (Figure 6(b)). On the contrary, the TLR4-ligand complex (Figure 6(a)) covers the least subspace and shows the least variation. In the 2 D projection plots of trajectories, the free TLR4 shows higher trace values of the covariance matrix than the complex. The plot of eigenvalues vs. eigenvector indices of the complex has lower values and a faster decay (Figure 6(c)) than that of the free protein molecule (Figure 6(d)). The overall analysis manifests again that the cystatin-TLR4 complex has greater stability as reflected from its least conformational changes due to decreasing collective motions. Altogether, the analysis of the RMSD, RMSF, Rg, SASA, and PCA indicates better stability of the cystatin-TLR4 complex compared to the protein itself (free TLR4).

Analysis of binding free energy and residue interaction energy
The binding free energy (DG bind ) was computed with MM/PBSA (Kumari et al., 2014), which is considered as a very key approach to study crucial molecular recognition processes (Kumari et al., 2014;Sen Gupta, Biswal, et al., 2020;Singh et al., 2020). Currently, this method has been widely applied in various studies including stability checks, protein-ligand interactions, proteinprotein interactions, mutational studies, and drug design and rescoring (Genheden & Ryde, 2015;Sen Gupta, Islam, et al., 2020;Singh et al., 2020;Wang et al., 2019).
The contributions of different interactions to the overall binding free energy are presented in Table 4 for the cystatin-TLR4 complex. Decomposition into separate energy terms indicates that the polar solvation energy opposes the binding of cystatin to TLR4 significantly and thereby increasing the total binding energy (Table 4). Among the various interactions, the electrostatic energy (DE elec ) is the most and the apolar solvation energy is the least favourable contributions towards the negative binding free energy of the complex (Table 4). In addition, the contributions of individual amino acids to the binding free energy (DG bind ) were also computed using the MM/PBSA method and are presented in Table 5. For clarity, only actively contributing residues with either positive or negative binding energy are presented. A per-residue interactional energy profile reveals that Asp405, Arg355, Asp428, His426, Phe377, Asp379, Tyr403, Tyr451, Ile450, and Ser381 from TLR4 actively participate in the interaction to give rise to a very strong binding and large stability of the complex. Interestingly, these residues were also interacting during the docking analysis, validating the molecular docking result. The four active residues contributing lesser than À5 kJ/mol to the cystatin-TLR4 complex, such as Asp405, Arg355, Asp428, His426, are the major contributors.

Cystatin induced conformational changes in TLR4
In addition to the above studies, we have compared the native structure of unbound TLR4 with the structures of cystatin-TLR4 complex taken from different timepoints of MDS (Figure 7). Our study clearly revealed the changes in TLR4 after cystatin binding at the extracellular domain of TLR4 (Figure 7). Surprisingly, the superimposed image of TLR4-cystatin complexes before and after MD simulation reveals key residues of TLR4 and cystatin required for the protein-protein binding maintaining a high stability throughout the MD simulation ( Figure S6).

Discussion
Cystatin is a small molecular weight immunomodulatory protein of filarial parasite that plays a pivotal role in downregulating the host immune response to prolong the survival of the parasite inside the host body (Bisht et al., 2019;Gregory & Maizels, 2008). Hitherto, this protein is familiar as an inhibitor of human cysteine proteases viz. cathepsin L, B, and S (Hartmann et al., 2002). Inhibition of the host proteases is pivotal in establishing a 'worm friendly' environment required for their survival inside the human (Babu & Nutman, 2014;Gregory & Maizels, 2008;Hoerauf et al., 2005). However, growing evidences on the role of cystatin in regulating inflammatory homeostasis prompted us to investigate the molecular reasons behind the explicit anti-inflammatory trait of this protein. In fact, a series of studies have already evidenced the strong anti-inflammatory activity of filarial cystatin and its isoforms (Gregory & Maizels, 2008;Khatri et al., 2020;Schierack et al., 2003). Owing to its explicit anti-inflammatory nature, filarial cystatin has been examined as a possible therapeutic molecule for treating inflammatory diseases like bacterium-induced sepsis, experimentally induced colitis, rheumatoid arthritis, etc. (Bisht et al., 2019;Gregory & Maizels, 2008;Khatri et al., 2020). However, the mechanistic insight of the induction of anti-inflammatory response has not been addressed to date. But previous researchers have reported that cystatin may have a role in regulating the function of the pattern recognition receptors (PRRs) of the host (Bird et al., 2009;Klotz et al., 2011).
TLRs are considered as the most influential PRRs in the human immune system (Mukherjee et al., 2016;Mukherjee, Huda, et al., 2019). Human TLRs have been grouped into two major classes viz. cell surface (TLR1, TLR2, TLR4, TLR5 and TLR6) and intracellular TLRs (TLR3, TLR7, and TLR8) (Mukherjee et al., 2016;Mukherjee, Huda, et al., 2019). Each TLR recognizes distinct pathogen-associated molecular pattern from the invading pathogens and shape the immune response against the pathogens (Mukherjee et al., 2016;Mukherjee, Huda, et al., 2019). In filarial infection, cell surface TLRs majorly TLR2, 4, and 6 play a crucial role in identifying the distinct filarial antigens and directing the immune network to eliminate the infection (Mukherjee et al., 2016;Mukherjee, Huda, et al., 2019). For example, TLR4 located on the innate immune cells like macrophages and dendritic cells recognizes filarial surface antigen namely bestrophin resulting in the induction of proinflammatory responses through classical macrophage induction, activation of dendritic cells and T-cells (Th1 and Th17) (Mukherjee et al., 2017;Mukherjee, Karnam, et al., 2019). However, filarial parasite also causes anti-inflammatory responses through activating alternative TLR4 signaling pathways (Goodridge et al., 2005;Maizels et al., 2018). Therefore, two major events in filarial immunomodulation are the activity of filarial cystatin and alternative activation of human TLR4. The question of whether these two proteins interact with each other or not is an open question to answer.
The present study is a maiden report that investigates the possible molecular interactions between filarial cystatin and human TLR4 which could be one of the reasons behind the induction of anti-inflammatory responses by the filarial parasite. Amino acids sequences of the filarial parasite W. bancrofti and the crystal structures of the human cell surface TLRs are available in the public domains (NCBI and PDB). So, we have employed in-silico approaches for studying the interaction between cystatin and TLRs. Firstly, phylogeny and biophysical characteristics of W. bancrofti cystatin were studied and filarial cystatin was found to be different both phylogenetically and structurally than that of the human (Figure 1 and 2). Next, we have exploited molecular docking to explore the interaction between W. bancrofti cystatin and human TLRs. Our in-silico data clearly revealed that filarial cystatin interacts with most of the cells surface TLRs of humans but the most significant binding occurs with TLR4 ( Figure 3 and Figure S3). Our in-silico data clearly revealed that cystatin strongly (binding energy=-93.5 ± 10 kJ/mol) interacts with the extracellular domain of TLR4 and the binding position is adjacent to MD2. The biophysical basis of the interaction i.e. binding of cystatin to TLR4 i.e. mediated by hydrogen bonding and hydrophobic interaction ( Figure 3 and Figure S3, Table 3). Furthermore, we have also identified Lys88 of filarial cystatin and Asp405 of human TLR4 as the predominant contributors of binding (Figure 3, Figure S3 and Figure S7, Table 3 and Table 5). All these evidences collectively demonstrate the strong interaction between cystatin and TLR4 which primarily indicates that cystatin acts as a putative ligand of TLR4. Taking a clue from the molecular docking data, we performed NMA and MD simulation studies to confirm the stability of the protein complex formed by cystatin and TLR4. NMA initially indicated towards greater stability and flexibility in the protein-protein interaction between cystatin and TLR4 ( Figure  4). MD simulation supported the postulations of NMA and revealed excellent stability of the cystatin-TLR4 complex (Figures 5 and 6). Both NMA and MD simulation indicated towards the change in the conformation of TLR4 after the binding of cystatin (Figure 4-7). Induction of conformational changes in TLR4 structure is considered as another important trait of a TLR4 ligand and conformational change in TLR4 plays a crucial role in activating the downstream signaling cascade (Latty et al., 2018). Herein, our experimental data clearly revealed that binding of cystatin does induce conformational changes in TLR4 (Figure 7) which could be correlated experimentally with its explicit anti-inflammatory activities in near future. Intriguingly, the mode of conformational changes in TLR4 by cystatin was found to be similar to the different TLR4 ligands reported from filarial parasites like ES-62 (Rajasekaran et al., 2017), chitohexose (Panda et al., 2012), and bestrophin (Mukherjee, Karnam, et al., 2019).
Hitherto, cystatin was known to be a key regulator in filarial immunomodulation and the immunomodulation has been reported to be executed via the induction of anti-inflammatory responses inside the host (Goodridge et al., 2005;Gregory & Maizels, 2008;Khatri et al., 2020;Maizels et al., 2018;Schierack et al., 2003). However, the molecular route behind the induction of such an anti-inflammatory response was unclear. Therefore, the present study links both the aforesaid studies by demonstrating cystatin as a TLR4 ligand which binds to human TLR4, induces conformational changes, and most likely activates the anti-inflammatory/ alternative pathway leading to immunomodulation. The possible mode of action of cystatin has been presented in Figure 8.

Conclusion
Taken together, our data indicated that cystatin appears to be a ligand of TLR4 and we hypothesize that cystatin-TLR4 interaction most likely to play a key role in activating the alternative activation pathways to establish an anti-inflammatory milieu. Intriguingly, the characteristics of the filarial cystatin were found to be dissimilar to that of the human which appears to be a key advantage in targeting this protein for therapeutic intervention of bancroftian filariasis. Thus, the study provokes the development of chemotherapeutics for targeting the cystatin-TLR4 interaction to disrupt the pathological attributes of human lymphatic filariasis. Similarly, the design of a vaccine against the filarial cystatin can also provide useful results in combating the infection. Our findings are expected to provide a novel dimension to the existing knowledge on filarial immunopathogenesis and it will encourage the scientific communities for experimental validation of the inferences of the present investigation.