Exploring Quercetin Hydrate’s Potential as an Antiviral Treatment for Oropouche Virus

: The Oropouche virus is an orthobunyavirus responsible for causing Oropouche fever, a disease that primarily affects thousands of people in South and Central America. Currently, no speciﬁc antiviral treatments or vaccines are available against this virus, highlighting the urgent need for safe, affordable, and effective therapies. Natural products serve as an important source of bioactive compounds, and there is growing interest in identifying natural bioactive molecules that could be used for treating viral diseases. Quercetin hydrate is a compound classiﬁed as a ﬂavonoid, which has garnered scientiﬁc attention due to its potential health beneﬁts and its presence in various plant-based foods. In this study, we aim to evaluate the in vitro antiviral activity of quercetin hydrate against the Oropouche virus (OROV). Furthermore, we intend to explore its mode of action through in silico approaches. The cytotoxicity and antiviral activity of the compound were assessed using Vero cells. In addition, in silico studies were also performed through molecular docking, molecular dynamics simulations, Molecular Mechanics Poisson–Boltzmann surface area (MM/PBSA), and quantum-mechanical analysis in order to evaluate the interaction with the Gc protein of OROV. The assay revealed that the compound was highly active against the virus, inhibiting OROV with an EC50 value of 53.5 ± 26.5 µ M under post-infection treatment conditions. The present study demonstrates that the compound is a promising antiviral agent; however, the mechanisms of action proposed in this study need to be experimentally veriﬁed by future assays.


Introduction
Viral infections pose a significant public health concern due to their impact on morbidity rates. Among these infections, Oropouche fever, caused by the zoonotic arbovirus Oropouche virus (OROV) [1,2], is particularly noteworthy. OROV belongs to the Peribunyaviridae family and is predominantly transmitted to humans through the bites of infected midges (Culicoides paraensis) prevalent in South and Central America [3]. Oropouche fever is a significant public health concern, primarily affecting populations in South America and the pan-Amazon region, particularly in low-income countries [1,2]. Currently, there are no specific antiviral drugs or vaccines available for the treatment of this infection. Hence, there is an urgent need to discover effective, safe, and affordable therapeutic options.
The genome of OROV is composed of three distinct segments: large (L), medium (M), and small (S). The L segment carries the genetic information for essential proteins involved in viral replication, including the RNA-dependent RNA polymerase (RdRp) and endonuclease proteins. The M segment acts as a precursor for glycoproteins and encodes two major structural glycoproteins, namely glycoprotein N (Gn) and glycoprotein C (Gc), in addition to two non-structural (NS) proteins, NS1m and NS2m. Finally, the S segment contains the nucleocapsid protein (N), which plays a crucial role in encapsulating the viral genomic RNA [4]. The Gc and Gn proteins are essential for facilitating the key processes of attachment and fusion, which are critical for the virus to enter the host cell. Consequently, these proteins become important targets for drug discovery, offering promising possibilities for therapeutic intervention [5].
Protein structures play a crucial role in elucidating the molecular mechanisms underlying biological processes at the atomic level. In this regard, bioinformatics tools have emerged as powerful resources for conducting in silico studies. One such tool, Al-phaFold2 [6,7], has gained considerable attention as a promising artificial intelligence-based approach for predicting three-dimensional protein structures. Additionally, molecular docking of compounds has been recognized as an important step in early drug discovery, enabling the assessment of key interactions and the binding mode of a ligand to its target [8][9][10][11][12]. Complementary to this, molecular dynamics simulations offer valuable insights into the stability of protein-ligand complexes in solution [13,14]. Furthermore, quantummechanical methods (QM) have demonstrated high accuracy in predicting binding energies and elucidating the critical interactions between ligands and receptors [15][16][17]. These methods not only identify key binding protein residues but also shed light on important regions within the ligand.
Throughout history, plants have been recognized for their therapeutic properties. Flavonoids, a diverse and abundant class of plant secondary metabolites [18], are widely distributed in various plant species and are responsible for the vibrant pigmentation observed in fruits, vegetables, flowers, and leaves [18,19]. Among the flavonoid compounds, quercetin hydrate belongs to the subclass of flavonols [20], a distinct group within the broader class of flavonoids, and has been found to possess a wide range of biological activities [21,22]. Quercetin, a flavonoid compound, has garnered attention for its notable antioxidant properties [23], which contribute to the neutralization of harmful free radicals and the protection of cells against oxidative damage [23][24][25]. Furthermore, quercetin exhibits anti-inflammatory effects [26], supports cardiovascular health [27,28], and displays promising potential in cancer prevention and management [29,30]. Notably, quercetin has been the subject of research investigating its antiviral properties, including its potential inhibitory effects against a range of viruses [12,[31][32][33][34][35][36][37][38][39].
In this study, our objective is to assess the antiviral activity of quercetin hydrate against the Oropouche virus (OROV). In addition, we conducted a comprehensive theoretical, computational analysis of the compound and employed in silico molecular dynamics simulations to explore the potential interaction between quercetin hydrate and the Gc viral protein. Our investigation aims to shed light on the mechanism underlying the antiviral activity of quercetin hydrate against OROV, providing valuable insights for future therapeutic interventions.

Cells, Virus, and Chemical
The Vero cell line (CCL-81 ATCC) was cultured in Minimal Essential Medium (MEM) (Gibco, Waltham, MA, USA) supplemented with 10% (v/v) heat-inactivated fetal bovine serum (FBS) (Gibco, Waltham, MA, USA), 100 U·mL −1 of penicillin, 0.1 mg·mL −1 of streptomycin, and 0.5 µg·mL −1 of amphotericin B (Gibco, Waltham, MA, USA). The cell cultures were maintained at 37 • C under a humidified atmosphere containing 5% CO 2 . In parallel, the C6/36 cell line was cultivated in Leibovitz-15 medium (L-15) supplemented with 10% FBS at 28 • C. The Oropouche virus (strain BeAn 19991) stocks were propagated in C6/36 cells and subsequently titrated in Vero cells to assess plaque formation unit, following the protocols detailed in Section 2.4 of this study. To conduct the experiments, quercetin hydrate (2-(3,4-dihydroxyphenyl)-3,5,7-trihydroxy-4H-chromen-4-one hydrate) from Sigma-Aldrich (Saint Louis, MO, USA) was prepared as a 200 mM concentration stock solution in dimethyl sulfoxide (DMSO) and stored at −20 • C until its intended use. These methods and materials provide the necessary foundation for the subsequent investigations in this research, and due acknowledgment has been given to the respective suppliers and previous research works cited in Section 2.4.

Evaluation of Cell Cytotoxicity
The present methodology was executed following established procedures, as previously reported [12]. In brief, a total of 5 × 10 4 Vero cells were cultured in 96-well plates and treated with varying concentrations of quercetin hydrate, ranging from 500 µM to 31.25 µM, for a duration of 48 hours. Afterward, 1 mg/mL of 3-(4,5-dimethyl-2-thiazolyl)-2,5-diphenyltetrazolium bromide (MTT) from Sigma-Aldrich (Saint Louis, MO, USA) was introduced to the cells, which were then incubated for 1 hour. Following the incubation period, formazan crystals were dissolved in dimethyl sulfoxide (DMSO), and the absorbance was measured at 550 nm using a Spectramax Plus microplate reader (Molecular Devices, Sunnyvale, CA, USA). The obtained results were expressed as the percentage of viable cells in the quercetin-hydrate-treated group relative to the untreated control cells. To ensure robustness, three independent experiments were conducted in triplicate. The 50% cytotoxic concentration (CC 50 ), indicating the concentration of the compound that reduced cell viability by 50%, was calculated from the dose-response curve using four-parameter curve-fitting in GraphPad Prism software (version 8.00).

Dose Response Assay
To elucidate the dose-response curve of quercetin hydrate on OROV, 24-well plates were utilized, with each well containing 1 × 10 5 cells. The cells were subsequently exposed to OROV (MOI = 0.1), and various concentrations of quercetin hydrate, ranging from 250 µM to 31.25 µM, were added one hour after infection. Controls comprised cells treated with 0.5% dimethyl sulfoxide (DMSO) and untreated cells. Consequently, postinfection treatment conditions were employed in Vero cells, and 48 hours post-infection (hpi), the supernatant virus plaque assay (as described in Section 2.4) was conducted. Three independent experiments were performed in triplicate, pool-titled. A four-parameter curve fitting analysis using GraphPad Prism software (version 8.00) was employed to calculate the half-maximal effective concentration (EC 50 )-the concentration of the compound inhibiting 50% of infection-along with its corresponding confidence interval (CI). Moreover, to assess the compound's selectivity, the selectivity index (SI) was determined as the ratio of CC 50 (cytotoxic concentration causing a 50% reduction in cell viability) to EC 50 .

Virus Plaque Assay
Vero cells grown in a 24-well culture plate were infected with 0.1 mL of 10-fold dilutions of supernatants. Following 1 h of incubation at 37 • C, 0.5 mL of culture medium supplemented with 2% FBS and 1.5% carboxymethylcellulose sodium salt (Sigma-Aldrich, Saint-Quentin-Fallavier, France) was added, and the incubation was extended for three days at 37 • C. After removing the media, the cells were fixed with 10% formaldehyde and stained with 2% crystal violet diluted in 20% ethanol. Plaques were counted and expressed as plaque-forming units per milliliter (PFU·mL −1 ) [12].

Protein Modeling
The OROV Gc protein was modeled using the AlphaFold2 Collab platform (https:// colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb# scrollTo=AzIKiDiCaHAn, accessed on 25 May 2023), a state-of-the-art deep learning-based method for protein structure prediction. The input sequence contained 385 amino acids retrieved from the Oropouche virus (strain BeAn 19991), and default settings were used. For our purpose, the transmembrane segments were removed. The residues ranged from 942 to 1326, based on the M protein sequence. The AlphaFold2 Collab utilizes a two-step process involving multiple neural network architectures to generate accurate protein structure models. In the first step, the network predicts residue-level distances and orientations. In the second step, the network refines the initial model by incorporating additional structural information. The AlphaFold2 Collab provides an efficient and reliable approach for protein structure prediction, leveraging the power of deep learning algorithms and large-scale training data. By utilizing this platform, we obtained high-quality protein structure models for further analysis and interpretation.

Quercetin Hydrate Exhaustive Docking
The 3D structure of quercetin hydrate was obtained from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/compound/16212154, accessed on 26 May 2023) in SDF format. The structure was converted to pdbqt format using the Open Babel 3.1.0 software to make it compatible with the AutoDock Vina 1.1.2 software [40,41]. For docking the molecule, a grid box was established based on previous studies [42] involving the Gc protein to define the binding site. Additionally, the modeled protein was converted to pdbqt format using the AutoDock Tools platform. An exhaustive docking analysis comprising 1000 runs was conducted, resulting in 10 ligand binding conformations per run. The best Vina Score obtained was selected for molecular dynamics (MD) simulation.

Molecular Dynamics Simulation
Two different systems were subjected to molecular dynamics (MD) simulations: one consisting of the AlphaFold2 output structure without any ligand and the complex obtained from molecular docking. For the first system, MD simulations were performed using GROMACS 2020.3 software [43] with the Amber ff99SB-ILDN force field. The protonation of histidines was verified using the PBD2PQR server (https://server.poissonboltzmann.org/pdb2pqr, accessed on 26 May 2023) at pH 7.4. The systems were placed in a cubic box with the TIP3P water model extended 12 Å away from the solute atoms and neutralized with Na + ions. Initial clashes in the starting structure were resolved through two rounds of energy minimization. The first minimization step involved a maximum of 500 steps or termination when the maximum force on any atom dropped below 50 kJ/mol/nm, employing the steepest descent algorithm with protein restraint to focus on solvent relaxation. The second minimization step was carried out without protein restraint in flexible water, using the same steepest descent algorithm. The number of maximum steps was increased to 10,000 steps or termination when the force on any atom fell below 250 kJ/mol/nm.
The system's pressure and temperature were initialized to 1 atm and 310 K, respectively, using two consecutive 100 ps steps known as the NVT ensemble (for temperature setting) and the NPT ensemble (for pressure setting). The modified Berendsen algorithm [44] was employed to control temperature, while the Parrinello-Rahman algorithm [45] was used to maintain pressure to achieve this. During both steps, hydrogen bonds were constrained using the LINCS algorithm [46], and positional restraints were applied to the protein to facilitate the equilibration of the solvent around the solute.
The long-range electrostatic interactions were computed using the Particle Mesh Ewald summation method, and a non-bonded interaction cut-off of 1 nm was applied. Integration of the equations of motion was carried out using the leap-frog algorithm [47] with a time step of 0.2 fs. Prior to the MD simulation, a short 1 ns NPT ensemble was conducted with no restrictions on the protein position. Subsequently, the production run for each identical system lasted 200 ns, also with no restrictions on the protein conformation. A total of 2000 protein frames were generated from each MD run.
Before conducting the MD simulation for the second system comprising the protein and ligand, the ligand parameters were generated. Initially, hydrogen atoms were added at pH 7.4 using the Avogadro program [48]. The parameters for the ligand were determined by submitting the ligand to the ACPYPE server (https://www.bio2byte.be/acpype/, accessed on 26 May 2023) [49]. The ACPYPE server utilized the BCC charge method and GAFF2 atom type to generate the topology and parameter files necessary for MD simulation using the GROMACS 2020.3 software. The MD steps of the first system were then conducted for the complex. Hence, two MD simulations (protein and complex) were obtained.
For the analysis of the MD simulation, we utilized native GROMACS packages (gmx rms, gmx rmsf, and gmx hbond). Additionally, we calculated the contact time between the residue and the ligand using an in-house Python script, where a contact was defined as a distance of up to 4 Å between the alpha carbon (Cα) of the residue and the ligand. The contact analysis was performed for each of the 2000 frames obtained in the trajectory, totaling 200 ns.

MM/PBSA Analysis
Following the generation of the MD trajectory, the analysis focused on the last 500 frames (50 ns) of the complexes. The Molecular Mechanics Poisson-Boltzmann surface area (MM/PBSA) free energy calculation was employed to evaluate the energetics, excluding all solvent molecules and Na + ions. The MM/PBSA was conducted by the gmx_MMPBSA program [50]. This software enables the direct calculation of MM/PBSA from the output of GROMACS MD simulations.

QM Binding Analysis Calculations
QM is a highly accurate method for studying intermolecular interactions; however, its application to large protein complexes can be computationally demanding. To address this challenge, Zhang and Zhang (2003) [51] developed the Molecular Fractionation with Conjugated Caps (MFCC) approach. This method involves splitting the system into individual amino acids by breaking the peptide bonds. In order to maintain the proper valence of the N-and C-terminals, the amino acids preceding and succeeding the target amino acid are included in the calculations. These neighboring amino acids are referred to as "caps". By performing calculations individually between each amino acid and the binder rather than for the entire system, the computational cost is significantly reduced. The following calculation seen in Equation (1) is performed to calculate the interaction energy (IE) between the binder (LIG) and the amino acid (R i ) to remove the influence of the caps: The ligand in this study is denoted as LIG, while the main residue is represented as R i , where i refers to the ith amino acid in the protein chain. To consider the influence of neighboring residues covalently bound to the amine and carboxyl group of R i , two caps, C i−1 and C i+1 , were introduced. In Equation (1), the first term (E(LIG + C i−1 R i C i+1 )) quantifies the interaction energy between the ligand and the main residue with caps. The second and third terms were subtracted to eliminate the contribution of the caps interactions. The second term (E(C i−1 R i C i+1 )) calculates the energy of the residue bound to the caps alone, while the third term (E(LIG + C i−1 C i+1 )) accounts for the interaction energy between the ligand and the caps. The fourth term (E(C i−1 C i+1 )) is included because the caps were removed twice (in the second and third terms) to isolate the effect of the interaction energy between LIG and R i from the caps' energy.
The interaction energy between the receptor and the binder was calculated using the Gaussian 16 package [52], employing density functional theory (DFT) formalism [53,54] and the B97D generalized gradient approximation (GGA) functional [55]. This GGA-type density functional has been demonstrated to be an efficient and accurate QM method for large systems, particularly those involving significant dispersion forces [56]. To expand the Kohn-Sham orbitals for all electrons, the 6-311 + G(d,p) basis set with triple split valence (valence triple-zeta) and an additional diffuse function (+) was utilized. The conductor-like polarizable continuum model (CPCM) [57] was employed to consider solvent effects in the QM calculations. The CPCM takes into account the dielectric constant (ε) of the solvent, and in this study, ε value of 40 was used. Subsequently, each term of Equation (1) underwent (DFT) calculations to obtain the corresponding results. To avoid overlooking important interactions, an analysis of the convergence of the total binding energy as a function of radius r was performed for the ligand binding pocket, limiting the number of amino acid residues to be studied. Thus, imaginary spheres were elaborated (considering r = n/2; n = 1, 2, 3, 4,. . .), limiting increasing distances of the ligands to obtain the sum of the individual energies of the residues present in each of these spheres. In this sense, convergence is achieved when the energy variation in the successive radius is less than 10%.

Cytotoxicity and Antiviral Activity
Cytotoxicity of the quercetin hydrate was determined in Vero cells using an MTT assay and performed in parallel with the replication inhibition activity measures by PFU to establish the antiviral effects. Cell viability corresponds to 93.6 ± 12.3% at 500 µM, which consequently makes it unable to determine a CC 50 value ( Figure 1A). No cytotoxicity is observed in cells treated with 0.5% DMSO (the final concentration of solvent used to dissolve the compound in a cell culture medium). The antiviral effect of the quercetin hydrate was investigated in vitro against OROV, which may threaten human health. Vero cell monolayers were treated with 500-31.25 µM quercetin hydrate or the equivalent volume of DMSO 1 h after OROV (MOI = 0.1) penetration into cells to assess the activity of the compound in a post-entry stage of the virus replicative cycle. Then, virus yields were measured by viral titration (PFU/mL). The presence of DMSO does not affect the production of progeny infectious virus particles. The EC 50 of quercetin hydrate is 53.2 µM ± 26.5, and the SI > 9.3. As demonstrated in Figure 1B, these results indicate a significant dosedependent decrease in the production of infectious OROV particles in the presence of increasing quercetin hydrate concentrations.
Biophysica 2023, 3, FOR PEER REVIEW 6 The interaction energy between the receptor and the binder was calculated using the Gaussian 16 package [52], employing density functional theory (DFT) formalism [53,54] and the B97D generalized gradient approximation (GGA) functional [55]. This GGA-type density functional has been demonstrated to be an efficient and accurate QM method for large systems, particularly those involving significant dispersion forces [56]. To expand the Kohn-Sham orbitals for all electrons, the 6-311+G(d,p) basis set with triple split valence (valence triple-zeta) and an additional diffuse function (+) was utilized. The conductor-like polarizable continuum model (CPCM) [57] was employed to consider solvent effects in the QM calculations. The CPCM takes into account the dielectric constant (ε) of the solvent, and in this study, ε value of 40 was used. Subsequently, each term of Equation (1) underwent (DFT) calculations to obtain the corresponding results. To avoid overlooking important interactions, an analysis of the convergence of the total binding energy as a function of radius r was performed for the ligand binding pocket, limiting the number of amino acid residues to be studied. Thus, imaginary spheres were elaborated (considering r = n/2; n = 1, 2, 3, 4, ...), limiting increasing distances of the ligands to obtain the sum of the individual energies of the residues present in each of these spheres. In this sense, convergence is achieved when the energy variation in the successive radius is less than 10%.

Cytotoxicity and Antiviral Activity
Cytotoxicity of the quercetin hydrate was determined in Vero cells using an MTT assay and performed in parallel with the replication inhibition activity measures by PFU to establish the antiviral effects. Cell viability corresponds to 93.6 ± 12.3% at 500 µM, which consequently makes it unable to determine a CC50 value ( Figure 1A). No cytotoxicity is observed in cells treated with 0.5% DMSO (the final concentration of solvent used to dissolve the compound in a cell culture medium). The antiviral effect of the quercetin hydrate was investigated in vitro against OROV, which may threaten human health. Vero cell monolayers were treated with 500-31.25 µM quercetin hydrate or the equivalent volume of DMSO 1 h after OROV (MOI = 0.1) penetration into cells to assess the activity of the compound in a post-entry stage of the virus replicative cycle. Then, virus yields were measured by viral titration (PFU/mL). The presence of DMSO does not affect the production of progeny infectious virus particles. The EC50 of quercetin hydrate is 53.2 µM ± 26.5, and the SI > 9.3. As demonstrated in Figure 1B, these results indicate a significant dosedependent decrease in the production of infectious OROV particles in the presence of increasing quercetin hydrate concentrations.

OROV Gc Protein Modeling
The protein structure predicted using the AlphaFold2 collab platform yielded a plDDT score of 86 and a pTM value of 0.669. The plDDT score represents the confidence level of the predicted local structure for each residue in the protein, ranging from 0 to 100. A higher plDDT score indicates a higher degree of confidence and more accurate predictions for the local structure. Similarly, the Predicted Template Modeling Score (pTM) quantifies the structural similarity between two folded protein structures, with a range of 0 to 1. A pTM value greater than 0.5 typically indicates a substantial level of similarity, enabling meaningful inferences to be made. Hence, the results obtained here describe a high-quality modeled protein. The structure of the protein in PDB format can be seen in Supplementary File S1.
In Figure 2A, the sequence identity of the query is depicted. The majority of sequences displayed a sequence identity greater than 0.5 for all amino acid positions. Figure 2B illustrates the plDDT score for each residue in the five models (ranked as rank_1 to rank_5). A significant portion of the residues exhibits scores higher than 80. Notably, the top-ranked model (rank_1) displays a predominance of residues with high plDDT values, as depicted in Figure 2C in ribbon structure representation. This structure was submitted to MD simulation, and the results will be discussed later.
was evaluated using the MTT assay. All data are presented with error bars, and three independent experiments were performed. (B) Quercetin hydrate causes dose-dependent inhibition of OROV yields. Vero cells were treated with only medium with 2% FBS (black bar), DMSO solvent control (white dotted bar) or 250-31.25 µM (pink bars) of quercetin hydrate in indicated conditions prior to infection with OROV at MOI = 0.1. The supernatant was titrated at 48 h post-infection. Mean values with standard deviation from three independent experiments are shown over each bar (*, p < 0.05).

OROV Gc Protein Modeling
The protein structure predicted using the AlphaFold2 collab platform yielded a plDDT score of 86 and a pTM value of 0.669. The plDDT score represents the confidence level of the predicted local structure for each residue in the protein, ranging from 0 to 100. A higher plDDT score indicates a higher degree of confidence and more accurate predictions for the local structure. Similarly, the Predicted Template Modeling Score (pTM) quantifies the structural similarity between two folded protein structures, with a range of 0 to 1. A pTM value greater than 0.5 typically indicates a substantial level of similarity, enabling meaningful inferences to be made. Hence, the results obtained here describe a high-quality modeled protein. The structure of the protein in PDB format can be seen in Supplementary File S1.
In Figure 2A, the sequence identity of the query is depicted. The majority of sequences displayed a sequence identity greater than 0.5 for all amino acid positions. Figure  2B illustrates the plDDT score for each residue in the five models (ranked as rank_1 to rank_5). A significant portion of the residues exhibits scores higher than 80. Notably, the top-ranked model (rank_1) displays a predominance of residues with high plDDT values, as depicted in Figure 2C in ribbon structure representation. This structure was submitted to MD simulation, and the results will be discussed later.

Quercetin Hydrate Docking in OROV Gc Protein
The docking of the quercetin hydrate compound was conducted within the highlighted region delimited in Figure 3. The selected binding pocket corresponded to the one observed in the Gc protein of the Rift Valley virus [42]. The most favorable quercetin pose, as illustrated in Figure 3, exhibits an AutoDock Vina score of −6.268 kcal/mol. This complex was subsequently subjected to molecular dynamics (MD) simulation to assess its stability in a solution environment. A summary of the results obtained from the exhaustive docking is presented in Table 1.
The docking of the quercetin hydrate compound was conducted within the highlighted region delimited in Figure 3. The selected binding pocket corresponded to the one observed in the Gc protein of the Rift Valley virus [42]. The most favorable quercetin pose, as illustrated in Figure 3, exhibits an AutoDock Vina score of −6.268 kcal/mol. This complex was subsequently subjected to molecular dynamics (MD) simulation to assess its stability in a solution environment. A summary of the results obtained from the exhaustive docking is presented in Table 1.

OROV Gc Protein MD Simulations and MM/PBSA Analysis
The Root-Mean-Square Deviation (RMSD) analysis of the OROV Gc protein revealed that the protein is more stable when complexed with quercetin hydrate. This is evident in Figure 4A, where, after 125 ns of simulation, the RMSD of OROV Gc-ligand is below 0.55 nm, while the protein alone exhibits a high fluctuation ranging from 0.25 nm to 0.85 nm. These findings indicate that the binding of the ligand helps to stabilize the protein dynamics in solution. Additionally, when analyzing the per-residue fluctuation ( Figure 4B), it is observed that most of the residues exhibit similar or lower fluctuations in the complex compared to the protein without the ligand, which further supports the results obtained from the RMSD analysis.
During the MM/PBSA analysis ( Figure 4C), it was observed that the majority of frames from the last 50 ns of simulation exhibited negative values of binding energy, indicating the stability of the complex in solution ( Figure 4B). The average binding energy, calculated to be −7.62 kcal/mol, suggests that the complex is energetically favorable. Furthermore, the lowest observed value of −18.04 kcal/mol (frame 1844) highlights the exceptional stability and strong binding affinity of the complex. Based on these findings, the most stable structure was selected for further analysis, specifically for conducting QM analysis using DFT. The selected structure in PDB format and the ligand in MOL2 format can be seen in Supplementary File S2 and Supplementary File S3, respectively.

OROV Gc Protein MD Simulations and MM/PBSA Analysis
The Root-Mean-Square Deviation (RMSD) analysis of the OROV Gc protein revealed that the protein is more stable when complexed with quercetin hydrate. This is evident in Figure 4A, where, after 125 ns of simulation, the RMSD of OROV Gc-ligand is below 0.55 nm, while the protein alone exhibits a high fluctuation ranging from 0.25 nm to 0.85 nm. These findings indicate that the binding of the ligand helps to stabilize the protein dynamics in solution. Additionally, when analyzing the per-residue fluctuation ( Figure 4B), it is observed that most of the residues exhibit similar or lower fluctuations in the complex compared to the protein without the ligand, which further supports the results obtained from the RMSD analysis.
During the MM/PBSA analysis ( Figure 4C), it was observed that the majority of frames from the last 50 ns of simulation exhibited negative values of binding energy, indicating the stability of the complex in solution ( Figure 4B). The average binding energy, calculated to be −7.62 kcal/mol, suggests that the complex is energetically favorable. Furthermore, the lowest observed value of −18.04 kcal/mol (frame 1844) highlights the exceptional stability and strong binding affinity of the complex. Based on these findings, the most stable structure was selected for further analysis, specifically for conducting QM analysis using DFT. The selected structure in PDB format and the ligand in MOL2 format can be seen in Supplementary File S2 and S3, respectively.
To evaluate the stability of the ligand, we analyzed the RMSD of the ligand itself. A gradual increase in the RMSD of the ligand during the trajectory, but the value consistently remained below 1.0 nm ( Figure 5A). These findings suggest that the ligand remained bound within the binding pocket throughout the entire MD simulation. To provide a more comprehensive visualization, we have included Supplementary Video S1, which shows the MD trajectory with the protein and ligand. The per-residue analysis selected the residues by distance (6 Å) from the ligand based on the first frame of the analysis (frame 1500-150 ns): Lys1028, Ile1029, Pro1030, Ala1031, Lys1032, Glu1033, Gly1034, Trp1035, Leu1036, Thr1037, Gln1066, Thr1146, Thr1160, Ile1161, Gly1162, Thr1163, Gly1164, and Phe1186. All energy contributions from the last 50 ns of these residues were calculated. The average result can be seen in Figure 4D. The residue that presents lower energy is Ile1161 (−1.56 kcal/mol). The other residues present binding energy higher than −1.0 kcal/mol in the following order of increasing energy: Gly1162, Leu1036, Thr1146, Phe1186, Trp1035, Thr1163, Gln1066, Glu1033, Ala1031, Lys1028, Thr1037, Ile1029, Pro1030, Gly1164, and Thr1160. Lys1032 (0.01 kcal/mol) and Gly1034 (0.3 kcal/mol) present positive values, indicating repulsive interaction.
To evaluate the stability of the ligand, we analyzed the RMSD of the ligand itself. A gradual increase in the RMSD of the ligand during the trajectory, but the value consistently remained below 1.0 nm ( Figure 5A). These findings suggest that the ligand remained bound within the binding pocket throughout the entire MD simulation. To provide a more comprehensive visualization, we have included Supplementary Video S1, which shows the MD trajectory with the protein and ligand.
Moreover, we analyzed the number of hydrogen bonds (Hbonds) formed between the ligand and protein ( Figure 5B). Throughout the trajectory, the Hbond count ranged from 0 to 6. However, it was observed that the ligand and protein predominantly formed two or three hydrogen bonds most of the time, as illustrated in the density histogram ( Figure 5B). Moreover, we analyzed the number of hydrogen bonds (Hbonds) formed between the ligand and protein ( Figure 5B). Throughout the trajectory, the Hbond count ranged from 0 to 6. However, it was observed that the ligand and protein predominantly formed two or three hydrogen bonds most of the time, as illustrated in the density histogram ( Figure 5B).
To increase the comprehension of protein-ligand interactions, the time of contact between the amino acids and the ligand was analyzed ( Figure 5C). A contact was defined as a distance less than 4 Å between the Cα of the residue and the ligand to simplify the analysis and reduce computational cost. For better readability in the graph, only those with more than 50 contacts are shown in the image. It is worth noting that residues Ala1031, Lys1032, Glu1033, Gly1034, Trp1035, Leu1036, Gln1066, Thr1160, Ile1161, Gly1162, and Phe1186 remained in contact with the ligand for more than 150 ns, with the ligand staying close to the amino acids Lys1032, Glu1033, Gly1034, Trp1035, and Ile1161 throughout almost the entire simulation (200 ns). This indicates that the ligand remained stable and coordinated by these residues with average binding affinity, as shown in Figure 4D.

Quantum Interactions Analysis of OROV Gc and Quercetin Hydrate
The interaction between the OROV Gc protein and quercetin hydrate complex was analyzed using DFT and MFCC to determine the interaction energy of each amino acid residue. The cumulative energy was calculated by summing the energy contributions from residues presented in each radius (r). Convergence criteria are observed to be met at approximately r = 4.5 Å (Figure 6A), where the cumulative energy is −27.23 kcal/mol. To ensure that all relevant residues were included in the analysis, the analysis was extended to r = 9.5 Å. As a result, a total of 39 residues are observed within a distance of 9.5 Å from the ligand, and the cumulative energy is −27.79 kcal/mol. To increase the comprehension of protein-ligand interactions, the time of contact between the amino acids and the ligand was analyzed ( Figure 5C). A contact was defined as a distance less than 4 Å between the Cα of the residue and the ligand to simplify the analysis and reduce computational cost. For better readability in the graph, only those with more than 50 contacts are shown in the image. It is worth noting that residues Ala1031, residue. The cumulative energy was calculated by summing the energy contributions from residues presented in each radius (r). Convergence criteria are observed to be met at approximately r = 4.5 Å (Figure 6A), where the cumulative energy is −27.23 kcal/mol. To ensure that all relevant residues were included in the analysis, the analysis was extended to r = 9.5 Å. As a result, a total of 39 residues are observed within a distance of 9.5 Å from the ligand, and the cumulative energy is −27.79 kcal/mol.  Figure 6B). The energy difference observed between MM/PBSA and QM arises from two primary factors. Firstly, the values presented in Figure 4D represent an average of the energy binding contributions from the last 50 frames of the MD trajectory in MM/PBSA, whereas QM ( Figure 6B) relies on a single frame. Furthermore, the calculation methods employed in MM/PBSA and QM significantly differ. MM/PBSA is based on classical mechanics, whereas QM utilizes quantum mechanics, which is known for providing more accurate calculations.  Figure 6B). The energy difference observed between MM/PBSA and QM arises from two primary factors. Firstly, the values presented in Figure 4D represent an average of the energy binding contributions from the last 50 frames of the MD trajectory in MM/PBSA, whereas QM ( Figure 6B) relies on a single frame. Furthermore, the calculation methods employed in MM/PBSA and QM significantly differ. MM/PBSA is based on classical mechanics, whereas QM utilizes quantum mechanics, which is known for providing more accurate calculations.

Discussion
Natural products have emerged as a source of active biomolecules for the treatment of various diseases, making them of great importance in drug discovery and development, with their offering a vast array of bioactive compounds that can interact with biological targets and pathways [58]. In recent years, there has been a growing interest in natural products as they continue to provide valuable leads for drug discovery. Numerous clinically relevant drugs have originated from natural sources, highlighting the significance of natural product-based research [58,59]. Consequently, natural products have demonstrated considerable promise in addressing challenging therapeutic areas such as cancer, infectious diseases, and neurological disorders [60,61]. The exploration of natural products as a source of active biomolecules holds immense potential in drug discovery and represents a promising avenue for the development of novel and effective treatments against a wide range of diseases [58,59]. Quercetin hydrate, a flavonoid compound abundantly found in various fruits, vegetables, and plants, holds significant importance as a natural product and potential therapeutic agent for the treatment of diseases.
Extensive research in the past years has shed light on its diverse biological activities and health benefits. Quercetin hydrate exhibits potent antioxidant properties, capable of scavenging free radicals and protecting against oxidative stress-induced damage [23,26,62,63]. Furthermore, it demonstrates anti-inflammatory effects by modulating key inflammatory signaling pathways, which are implicated in the pathogenesis of chronic inflammatory diseases [26,62,64]. Additionally, quercetin hydrate has been shown to possess anticancer properties through its ability to inhibit cell proliferation, induce apoptosis, and impede angiogenesis [29,30,65]. Moreover, its potential antiviral activity against a variety of viruses, including respiratory viruses, has gained attention [12,31,33,34,36,37,39]. Quercetin hydrate's broad spectrum of biological activities makes it a promising candidate for the prevention and treatment of diseases, and for this reason, we investigated its potential antiviral activity against OROV.
Results from our study suggest that quercetin hydrate exhibits potent in vitro antiviral activity against intracellular virus replication (SI > 9.3) without relevant cytotoxicity. These observations can suggest that one of the possible mechanisms for quercetin hydrate intracellular activity against OROV could be attributed to its ability to bind and/or to inactivate important structural and/or non-structural protein(s) of OROV. Thus, we proposed the interaction of the compound with the Gc protein of OROV as a possible mechanism of action. The Gc protein is a crucial component of the OROV virion envelope and plays a significant role in viral entry and fusion, budding, and infectivity [66]. Thus, the Gc protein can be an interesting target for the rational development of drugs, and the present study was a pioneer in proposing such a target for antiviral studies. However, it is essential to acknowledge that our current study's scope is limited to confirming the binding of quercetin hydrate solely to Gc or excluding its interaction with other viral proteins or subunits. Further investigations of a more intricate nature would be required to address these assumptions comprehensively, and we encourage the pursuit of such research endeavors in the future. By conducting these additional studies, a deeper understanding of the multifaceted role of quercetin hydrate in the context of viral infection could be elucidated, shedding light on its potential implications and applications in the field of virology. While OROV is an important cause of mosquito-borne viral disease in South and Central America, research on antiviral interventions against this specific virus may be limited compared to more extensively studied viruses.
The in silico analysis revealed that the ligand binds to the same region as previously observed in studies on Gc of the Rift Valley virus. The MD and MM/PBSA analyses demonstrated the stability of the complex in solution, with favorable binding energy throughout the last 50 ns of the MD trajectory. QM analysis using the DFT functional was performed to evaluate the binding energy of the best complex observed in the MD trajectory to provide a more accurate assessment. The results revealed a binding energy of −27.79 kcal/mol, indicating that Gc is a promising target for inhibiting OROV membrane fusion. Additionally, the analysis highlighted important residues for ligand stability, offering insights into the inhibition site of OROV Gc.
Computational studies are widely recognized as an important initial step in drug design and discovery. However, it is important to acknowledge that these studies have limitations. Therefore, further research is essential to fully understand the mechanisms of action of the ligand, optimize its bioavailability, and evaluate its efficacy and safety in clinical settings.

Conclusions
In summary, our study demonstrated that quercetin hydrate exhibited potent anti-OROV activity in vitro. In addition, using methods of virtual modeling, we have acquired a model of Gc protein that can be used for further development of new inhibitors of OROV reproduction by means of molecular docking. The interaction model of OROV-Gc with quercetin hydrate was observed by molecular docking, molecular dynamic simulation, and mechanical-and quantum-interaction analysis. Results from this study warrant future in vivo antiviral, toxicity, and pharmacokinetic studies as part of the developmental process for the development of quercetin hydrate as a potential anti-OROV therapeutic.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/biophysica3030032/s1, File S1: OROV Gc protein structure modelled by AlphaFold; File S2: OROV Gc protein in complex with quercetin hydrate obtained in the MD simulation for QM analysis; File S3: Quercetin hydrate structure after protonation at 7.4 pH. Video S1: MD trajectory movie of OROV Gc in complex with quercetin hydrate.

Data Availability Statement:
The data for OROV Gc protein modelling are openly available in NCBI/GenBank at https://www.ncbi.nlm.nih.gov/protein/, reference number AJE24679. The modelled protein structure and the frame selected for QM analysis is available in Supplemenaty Materials. The data that support the findings of MD/MM-PBSA and QM/MFCC are available from the corresponding author, Umberto Laino Fulco, upon reasonable request due to large size files.