The novel anti-cancer feature of Brazzein through activating of hTLR5 by integration of biological evaluation: molecular docking and molecular dynamics simulation

Many of plant proteins exhibit the properties similar to the antitumor proteins although the anticancer activity of Brazzein on modulating the autophagy signaling pathway has not been determined so far. The present study aimed to develop a simplified system to enable the rational design of the activating extracellular domain of human Toll-like receptor 5 (hTLR5). To identify the anticancer effect of Brazzein, HADDOCK program and molecular dynamics (MD) simulation were applied to examine the binding of the wild type (WT) and p.A19K mutant of Brazzein to the TLR5. The expression of MAP1S and TNF-α genes was estimated based on real-time PCR. The results clearly confirmed that the WT of Brazzein activated hTLR5 in the MCF-7 cell line since the genes were more and significantly less expressed in the cells treated with the WT and p.A19K mutant than the control, respectively. The snapshots of MD simulation exhibit the consistent close interactions of hTLR5 with the two helices of Brazzein on its lateral side. The results of per residue-free energy decomposition analysis substantiate those of intermolecular contact analysis perfectly one. We propose that the WT of Brazzein can act as an antitumor drug candidate.


Protein-protein interaction (PPI). The interactions between either the WT or mutant with hTLR5
were assessed through applying the molecular docking using a HADDOCK algorithm 22 . Table 2 lists the free energy of docking values for the WT-and mutant-hTLR5 complexes. Based on the results, the minimum Z-score was observed in the WT, reflecting that the p.A19K substitution can destabilize the brazzein-hTLR5 complex ( Table 2). Thus, MD simulation was employed to confirm this hypothesis.
The 3D structure of the favorable complex is displayed in Fig. 1. As shown, the concave β-sheet domains of TLR5 and conserved α2-helixe sites of Brazzein form the binding interface, leading to higher interface distortion and more buried residues, which strengthens the binding of the two proteins. Further, Ala substitution to Lys reverses the structure of Brazzein in the binding pocket of hTLR5 (Fig. 1A,B), which is likely to unstabilize the interaction between p.A19K mutant and hTLR5. The subsequent analyses could confirm this hypothesis.
Effects of mutation on the binding pocket size. CASTp server measures all surface and internal pockets, and their exact volume and area, as well as the size of openings 23 . Before the mutation, the area of the binding pocket 1 was equal to 178.725, which elevated to 2351.578 following the mutation and covered almost the entire molecule (Table 3).

MD simulation analysis.
The Rg values for the protein were plotted against the 100 ns of MD simulation to analyze the compactness, shape, and folding of the overall structure of WT and p.A19K mutant in complex with hTRL5 ( Fig. 2A). In fact, the Rg exhibits the compactness and overall dimension of protein or macromolecular www.nature.com/scientificreports/   Fig. 2A, the Rg plot fluctuated during 30-80 ns simulation time. It seems that the structure of the complexes may be altered during simulation. Furthermore, the WT-hTLR5 complex had a mean Rg of around 3.63 nm, while a slightly lower Rg of about 3.55 nm was obtained in another, which represents the greater compactness of the system during the 100-ns simulation. The similar Rg for both of the structures suggested the tight packing of the WT-hTLR5 complexes, making the structure relatively stable ( Fig. 2A). The Rg curve compression of Brazzein with TLR5 is similar to the RMSD and RMSF parameters, indicating that Brazzein tries to reach internal configuration in hTLR5. Two MD simulations were performed for WT and p.A19K mutant at 37 °C to compare proteins in terms of structure and dynamic behavior. Additionally, the RMSD values for the backbone atoms of the proteins were monitored relative to the starting structures over the simulation time. As depicted in Fig. 2B, the WT structure reflects a rapid rise to 1.14 nm following the beginning of the simulation run. The substitution of A 19 with K leads to lower fluctuations than the WT (≈ 50 Å) although the greatest flexibility belongs to the WT during the 50-70 ns of running, which is in line with the results of Rg plot after 80 ns. Further, the fluctuations are very small and both trajectories have almost the same values. Therefore, the structures are in their equilibrium states, which are highly correlated with certain interaction energies such as van der Waals force and non-polar desolvation 25 .
The RMSF analysis can be used as a tool to describe local flexibility differences among residues throughout the MD trajectory 26 . It was expected that the catalytic domain of hTLR5, including the residues 23-664 in complex with Brazzein, exerts various structural effects on the different regions of the protein (Fig. 2C,D). The residues had well-known interactions with hTLR5. The peak fluctuation values were 0.5 and 1.8 nm for the WT and p.A19K mutant, respectively. In both complexes, the N-terminal (NT) end was quite flexible than the C-terminal (CT) one (Fig. 2C), while WT-hTLR5 complex possessed a slightly more stable CT end region in comparison with another. The CT end region is a fragment of K15-N20 in the Brazzein structure, which was important for binding to hTLR5 and underwent hotspot mutation. The result is in good agreement with that of Ghanavatian et al. 5,27 . which revealed a region of protein containing position 19 as one of the critical points for interacting with the receptor. Seemingly, any change in this position may alter the stability of the corresponding complex. The B-factors or thermal fluctuations of the complexes modeled from MD simulation were computed to corroborate the results of RMSF analysis. The results of RMSF analysis are perfectly correlated with the B-factor one where the central loops illustrate a high degree of deformity along the NT and CT regions. Furthermore, the PCA was carried out on the MD trajectory during 100 ns to understand the motion changes in both complexes in more details (Fig. 2E,F). www.nature.com/scientificreports/ To comprehend the motion directions captured by the two top-ranked principal components (PCs) quantitatively, porcupine plots were generated by using the extreme projections on first (PC1) and second PC (PC2) (Fig. 3). In the cases, the length and direction of arrow denote the strength and direction of motion, respectively. The porcupine plots suggest the high-degree inward motion in the CT and NT regions of WT and p.A19K mutant in PC1. In the case of PC2, inward motion is only observed in WT, while the CT region of p.A19K mutant has outward motion. The results of PCA suggested the terminal end of motion modes, as well as central loops as the major flexible regions in p.A19K-hTLR5 complex, which are completely consistent with those of RMSF and B-factor analysis.
Essential dynamics. In the present study, the PCA was employed to describe fluctuation, flexibility, and atomic positional fluctuations for quantifying the main correlated motions of the WT and mutant, as well as identifying major dynamic protein regions (Fig. 4). The covariance matrix and 2D projection of simulations were mathematically assessed for the brazzein-hTLR5 complex by using two PCA eigenvectors during the 100ns simulation (Fig. 4). In this plot, a positive correlation between two residues reflects a determined motion along the same direction (blue), while a negative one refers to the opposing directions of motion (red) 26 . Based on the results in porcupine curves, the most correlated motions were found in the WT, the majority of which was centered on the interaction between Brazzein and TLR5 ( Fig. 4A,B). As shown in Fig. 4, the PC1 and PC2 between interacting residues indicate a twisting motion. The range of motion varies − 0.534 to 1.35 in the WT, while it is larger in the mutant (− 0.545 to 3.24). In general, the blue and red spots are more focused on the WT although they are broadly scattered and less intense in the mutant. In addition, more uncorrelated motions are detected in the mutant one. The results of correlation matrix analysis of the p.A19K mutant represent more scissoring motions along the PC1 axis compared to the WT (Fig. 4B). At the qualitative level, the PCA scatter plots obtained by using the first two PCs (PC1 and PC2) for the both models demonstrate the differences between the two types of eigenvectors ( Fig. 4C-E). The calculated eigenvectors are displayed as bars in Fig. 4C, which clearly exhibits that the first four eigenvectors account for more than 80% of the collective motions of Cα protein atoms. Further, the first two PCs possess the highest range of motion. The values related to PC1 and PC2 are respectively equal to 0.03174 and 0.01554 nm 2 in the WT, as well as 0.03628 and 0.01441 nm 2 in the p.A19K mutant. Furthermore, the comparison of the values reflects a greater range of motion in the simulation of p.A19K mutant. The scatter plot of PC1 against PC2 in phase illustrates the occupation of less conformational space by WT-hTLR5 complex than the p.A19K-hTLR5 one. Finally, p.A19K-hTLR5 complex occupies wider conformational space than the experimental complex due to the high degree of flexibility (Fig. 4F).   However, two residues like Ala 283 and Thr 310 were involved in the hydrophobic interaction with hTLR5 following the replacement of residue Ala 19 with Lys. Further, major differences were observed with respect to the number of hydrophobic contacts so that the WT-hTLR5 complex had more electrostatic and pi-alkyl contacts compared with another ( Fig. 5A-D). Given the high thermal stability and low chemical stability of Brazzein reported in the previous studies, the presence of disulfide bonds profoundly influenced the structural integrity of the protein 5 . In summary, the antithetical properties of p.A19K mutants on the protein activity revealed that the presence of charged residue at this position may exert a profound effect on the interaction between protein and corresponding receptor. Furthermore, the number of the positively-charged amino acids in the WT can enhance the interaction strength toward increasing the sweetness of protein. Seemingly, residues Ser, Gln, Asn, Tyr, and Glu in the p.A19K mutants were directly implicated in forming hydrogen bonds with other residues and polar water molecules. However, the interactions distorted β-strands, and decreased the stability of β-sheets or overall conformational stability of protein. The above-mentioned issues and gene expression results introduced local flexibility and positive charge at the protein surface as the more important factors affecting the protein interaction with the G-protein-coupled human sweet receptor. Therefore, this mutation destabilized the CT of hTLR5 binding to Brazzein. The results are in line with those of the previous studies 5 . Figure 5E demonstrates the WT structure (green) superimposed on either the p.A19K mutant (yellow) in complex with hTLR5 (WT in blue and mutant in red colors) to compare fluctuation level in the LRR and CT domains during the simulation. As depicted, the p.A19K mutant structurally changes over 30-80 ns, and exhibits more flexible structure compared to the WT, which are confirmed by using the results of RMSD curve.

Interaction of Brazzein with hTLR5. To gain deep insights into the interactions mediated by
The results of detailed structural analysis indicated the presence of more H-bonds in the WT than the p.A19K, which may alter the conformation of several overlying residues at the LRR and CT domains.
Regarding both the complexes, the interatomic distance between the crucial interacting pairs was measured from the MD trajectories (Fig. 5F,G). The mean distance between the interacting pairs of atoms in Brazzein and hTLR5 are outlined in Fig. 5H. The distance profile represented the continuous participation of the important amino acid residues in the molecular recognition process of Brazzein by hTLR5. The mean distance of the residues Phe 200 and Glu 254 elevated in the presence of hotspot mutation. Generally, the MD simulation of both hTLR5 complexes revealed a stronger interaction between the WT receptor and hTLR5.  www.nature.com/scientificreports/ of the PCA scatter which reflected the probability distribution of the first two eigenvectors and fluctuations of Brazzein mutant. The comparison between the conformational entropy of the complexes based on the PC1 and PC2 produced from the PCA indicated the highest entropy in the p.A19K mutant (Fig. 6B). As demonstrated in Fig. 6C, the minimum energy clusters are highlighted with black dotted circles. The results of the analysis exhibited a centralized energy distribution, suggesting the overall conformational stability of this system. Thus, the mutation of the residue, p.A19K, influenced the overall conformational stability of the brazzein-hTLR5 complex. The variations in the cumulative motion correlation profile for the p.A19K-hTLR5 complex revealed the relative differences in both landscapes.
Effect of Brazzein on the MCF-7 cell viability. n this study, the cells were incubated with the various concentrations of the WT and p.A19K mutant for 48 and 72 h. Based on the results, the cell viability increased in WT (85, 72, 58, and 47%) and p.A19K mutant (88, 78, 61, and 54%) after incubating for 48 h (Fig. 7A). Therefore, 80 μg/ml level of Brazzein was better in the WT than the p.A19K mutant. A comparison between the WT protein and p.A19K mutant treated for 72 h is provided in Fig. 7B. The WT at 20, 40, 60, and 80 μg/ml amounts reduced cell viability by 68, 62, 48, and 34%, respectively. However, 72, 66, 52, and 42% decline was observed following the incubation with the same concentrations of p.A19K mutant, respectively. Given that both the WT and p.A19K mutant had IC50 of 60 μg/ml level in the 72-h treatment, this concentration was used for further evaluation.
Expression of TNF-α and MAP1S genes. The TLR5 signals through the MYD88 pathway, which drives the production of inflammatory cytokines such as tumor necrosis factor-alpha (TNF-α) and interleukin-12 (IL-12) 29 . The antitumor mechanisms of hTLR5 signaling in breast cancer cells suggest MAP1S as an important autophagic adapter, which is associated with tumor suppression via autophagy regulation. MAP1S level promotes in response to flagellin stimulation in the MCF-7 cells 16 .
Compared to the control group, TNF-α gene was significantly less and more expressed in the p.A19K mutant (P < 0.05) and WT, respectively (Fig. 7C,D). In terms of MAP1S gene expression, the results reflected a significantly decrease in the p.A19K mutant (P < 0.01) and an enhancement in the WT (P < 0.05) compared to the control (Fig. 7E,F).

Discussion
This study provided a molecular rationale for understanding the conformational changes resulted from binding Brazzein to hTLR5 in the breast cancer cell line better. It presented the binding mode of Brazzein, a sweet protein, to PRRs for the first time. Despite the possible direct effect of the mutant on binding properties, the study considered that it is more likely to affect the conformational change of hTLR5 required to activate second messenger signaling within the cell. Similarly, the point mutation in Brazzein failed to act as presumed. In other words, the mutation did not cause autophagy by reducing gene expression. However, the WT protein led to the autophagy of breast cancer cells by elevating gene expression and activating the TLR5 pathway, which may be a good candidate as an anticancer peptide. Both the p.A19K mutant and WT similarly influenced breast cancer proliferation, while the p.A19K mutant did not enhance the hTLR5 signaling genes. Further, the strong correlation between the computational estimations and gene expression obtained by RT-PCR through cell-based responses upon exposure to Brazzein and hTLR5 strengthens the idea of using molecular docking and dynamics for designing and predicting the stability of their interaction. Furthermore, the WT can be introduced as the hTLR5 agonist, as well as can better respond to hTLR5 than the p.A19K, which is in line with the modeling results. It is proposed that future studies focus on the experimental assessment of this result through using the final Brazzein-receptor model to design the new mutants of the protein.

Materials and methods
Materials. Recombinant E. coli strain Shuffle ® T7 was utilized for gene expression in the wild type (WT) and p.A19K mutant of Brazzein 5 . In addition, MCF-7 human breast cancer cell lines were obtained from the Pasteur Institute of Iran. DMEM-high-glucose, fetal bovine serum (FBS), trypsin-EDTA, antibiotic (Pen-Strep), and phosphate buffered saline (PBS) were purchased from Inoclon (Iran). Further, MTT was provided by Merck (Germany), and RNAse extraction and cDNA synthesis kits were acquired from GeneAll (Korean) and Pishgam (IRAN), respectively. The other chemicals applied in this study were commercially available with the highest analytical grade.
Homological modeling. The molecular docking analysis was employed to examine the interaction between Brazzein and hTLR5 molecules. The 3D structure of the WT (PDB code: 2LY5) as a ligand 30 and human TLR5 (PDB code: 3J0A) as a receptor 31 was extracted from the protein data bank (http:// rcsb. org/). Furthermore, a mutation of alanine (Ala) to lysine (Lys) at the 19th position was generated in silico in the WT crystal structure using Modeller 9.20 software 32 . The general stereochemical quality of the final modeled protein structure was evaluated based on the parameters of Verify 3D, Z-DOPE, and RMSD obtained using SAVES (https:// saves. mbi. ucla. edu/) and ModEval servers, (https:// modba se. compb io. ucsf. edu/ evalu ation/) and SpdbViewer software, respectively 33,34 . The structures were also analyzed by the Ramachandran Plot Server (https:// zlab. umass med. edu/ bu/ rama/) to be sure that the selected structure has the minimum questionable positions, as compared with the other structures 35  www.nature.com/scientificreports/ The protein structural calculations were guided by defining active residues based on the previous report 28 , followed by selecting the structures with the least HADDOCK score and best binding site. Then, the program for automatically plotting protein-protein interactions (Ligplot-Dimplot, v 2.2.5) was applied to determine intermolecular contacts such as hydrogen bonds and non-bonded ones 36 . Finally, the 3D models of docking and molecular dynamics (MD) results were analyzed by using Chimera 1.13.1 software 36 .

MD simulation.
The computational studies on the biological macromolecules, provides a plethora of various data types as the output that can be more analyzed and processed 37,38 .
In this study, MD simulation was performed to evaluate the motions and fluctuations of two systems under the effect of the internal and external forces generated by hotspot mutant residue in Brazzein. Additionally, 100-ns simulations of the p.A19K mutant, detecting the conformational changes in the functional binding domain of Brazzein in complex with hTLR5, were independently created by using the GROMACS MD package (v 2021.2) 39 in the CentOS Linux system, and implemented with the GROMOS96 force field (v 43a1). This MD protocol has been previously described in the literature 40 . Briefly, simulated proteins were solvated in a 1.0 nm cubic box in the presence of Na + and Cl − ions through using SPC/E water model. The number of solvent molecules and total charge added to each system were 66,030 and 10 ions for wild type and 75,034 solvent and 10 ions for mutant. All of the covalent bonds to hydrogen atoms were constrained based on the SHAKE algorithm. Further, electrostatic interactions were calculated according to the particle-mesh Ewald (PME) algorithm 41 , and pre-equilibration was carried out with approximately 50,000 steps of steepest-descent minimization. After minimizing energy, the system was equilibrated by using the 100-ps simulation of a NVT (canonical) followed by a NPT simulation (P = 1 bar) and gradual heating from 310 K. The Berendsen thermostat 42 algorithm was utilized to maintain constant temperature and pressure 43 . Furthermore, 50,000 steps of minimization and 100 ns of unrestrained MD were performed. Then, short-range electrostatic and van der Waals interactions were computed with a distance cut-off of 1.0 nm 18,44,45 . After equilibration, coordinates were saved every 2 ps during the sampling process. The conformations generated from the simulations were used in further analyses. Finally, the coordinates and periodic corrections were examined by installing trjconv tools in the GROMACS software package. We have not capped the N-and C-terminals during simulation, as terminal residues are seen in the structure and are present in the PDB coordinates.
Trajectory analysis. The comparative equilibration was determined from a flattening of the radius of gyration (Rg), as well as root mean square deviation (RMSD), residue-based root mean square fluctuation (RMSF), secondary structure modifications, between-atom distances, free-energy landscape (FEL), entropy, principal component analysis (PCA) 46 , and eigenvectors (Evs) over time after a 100-ns interval for WT and mutant model. Data values are expressed as % of control values. Significant change (*P < 0.05, **0.05 < P < 0.01 and ***P < 0.001) in comparison with control group). Real-time PCR analysis for TNF-α and MAP1S genes in MCF-7 cells. The values put onto each graph represent the relative fold change calculated by calibrating the ΔΔCt data TNF-α (C,D) and MAP1S (E,F) gene expressions, respectively. Error bars indicate mean ± standard deviation (SD). The significant level was define as *P ≤ 0.05 and **0.05 < P < 0.01. cDNA synthesis. The total RNA from each sample was used to form cDNA with oligo (dT) primers according to the manufacturer's protocol (cDNA synthesis kit, Pishgam, IRAN).

Scientific
Oligonucleotide primers. This study utilized TNF-α, MAP1S, and β-actin (housekeeping gene). All of the primers were designed with Gene Runner, the sequences of which are presented in Table 4.

Molecular analysis: quantitative real-time polymerase chain reaction assay. The TNF-α and
MAP1S genes were evaluated using the housekeeping gene of β-actin. Further, the primers were previously checked through employing conventional real-time polymerase chain reaction (RT-PCR) and agarose gel electrophoresis (1.5%). The quantitative RT-PCR (qRT-PCR) was performed on a final volume of 20 μl containing 4 μl of Eva Green qPCR Mix, 1 μl of cDNA, and 1 μl of each of forward and reverse primers ( Table 4). The cycling protocol involved an initial denaturation step at 95 °C for 15 s, followed by 40 cycles of denaturation at 95 °C for 15 s, annealing at 52 °C for 15 s, and extension at 72 °C for 15 s. Furthermore, the CT values obtained based on the CT relative quantification method were normalized for endogenous 48 . The fold-change in gene expression was computed using the melt curve approach and normalized to β-actin. Ultimately, the relative gene expression levels were calculated with reference to the control 48 .
Statistical analysis. The data were analyzed through implementing one-way ANOVA, followed by Tukey's post-hoc test for comparing the mean of each group and an unpaired t-test between groups by using the Prism software (v 9.00) for Windows (GraphPad Software, RRID:SCR_002798, La Jolla, CA) for curve generation and statistical analysis. The data were expressed as mean ± the standard error of the mean (SEM) and p values less than 0.05 were considered as statistically significant difference. Each experiment was repeated in triplicate.

Data availability
The datasets of the 3D structure of the WT as a ligand and human TLR5 used during the current study are available in the PDB code: 2LY5 and PDB code: 3J0A, respectively (http:// rcsb. org/). The general stereochemical quality of the final modeled protein structure was analyzed based on the parameters of Verify 3D, Z-DOPE, and RMSD obtained using SAVES (https:// saves. mbi. ucla. edu/) and ModEval servers, (https:// modba se. compb io. ucsf. edu/ evalu ation/) and SpdbViewer software, respectively.