key: cord-0022292-p3td176k
authors: Chatterjee, Shilpa; Kim, Choon-Mee; Kim, Dong-Min
title: Potential efficacy of existing drug molecules against severe fever with thrombocytopenia syndrome virus: an in silico study
date: 2021-10-21
journal: Sci Rep
DOI: 10.1038/s41598-021-00294-7
sha: 3808a4ea85d6fb9d646c24a432b15d493d59f4cb
doc_id: 22292
cord_uid: p3td176k

Severe fever with thrombocytopenia syndrome (SFTS) is a zoonotic disease caused by the SFTS virus (SFTSV). SFTS can be considered a life-threatening notifiable infectious disease. The unavailability of specific therapeutics encourages the investigation of potential efficacy of existing drugs against this infection. Drug repurposing was done by performing virtual screening of already established drug molecules followed by 100 ns molecular dynamics simulations and molecular mechanics Poisson–Boltzmann surface area–based binding-energy calculation by targeting the SFTS L protein. On the basis of binding energy and protein–ligand interactions, top 10 promising hits were identified, showing stable binding with SFTS L protein. Further 100 ns atomistic MD simulation refined the hits from top 10 to top 4 with docking-based binding energy lesser than −8.0 kcal/mol toward the SFTS L protein and engaged in π–π interactions with pivotal amino acid residues. Various parameters and binding affinity of top 4 ligands towards L protein was computed. Ligand zaltoprofen exhibited best binding energy −220.095 kJ/mol. The present work is the first in silico study to assess bromfenac, cinchophen, elliptinium, and zaltoprofen; four promising hits against SFTS. Nonetheless, further proper biological evaluation is necessary to determine their efficacy against SFTS.

www.nature.com/scientificreports/ a capped, mostly nonpolyadenylated viral mRNA that contributes to viral replication and transcription. The mechanism behind SFTSV genome replication is believed to be a de novo process, but the transcription of SFTSV genes is solely dependent on the phenomenon of cap-snatching, where cellular RNAs are involved 8 . Targeting the L segment is believed to suppress SFTSV replication by inhibiting viral transcription and replication 5 . Hence, the multifunctional and multidomain characteristics of the L segment of bunyaviruses make it an ideal drug target in SFTSV. Vogel et al. (2020) reported the structure of the cap-binding domain of the SFTSV L protein complexed with a triphosphate inhibitor: 7-methyl-guanosine-5′-triphosphate (MGP) 6 .

It is a well-established fact that drug discovery is expensive and time consuming. Therefore, considering the urgency of finding effective drug candidates against an emerging disease such as SFTS, we decided to perform a drug-repurposing study by screening a database of FDA-approved small-molecule. Drug repurposing is an approach to drug development whereby we can reuse triphosphate like drugs already proven to be effective against other diseases 9 . This approach is becoming an efficient, cost-effective, and universal strategy with numerous advantages. In the present study, structure-based virtual screening of the DrugBank database was performed against the SFTSV L protein. Molecular dynamics (MD) simulation was performed to investigate the stability of ligand-receptor complexes. Binding free energies in terms of molecular mechanics Poisson-Boltzmann surface area (MMPBSA) were calculated too. The present work outlines the findings about four promising hits-bromfenac, cinchophen, elliptinium, and zaltoprofen-featuring possible binding to the SFTSV L protein.

Protein data preparation. This step is an important part of the in silico drug design approach. The SFTSV L protein's structure was retrieved from the Protein Data Bank (PDB; https:// www. rcsb. org/ search) (PDB ID: 6XYA). The protein crystal structure contains bound co-crystalline water molecules that were removed by using the PyMol software 10 . This co-crystalline protein structure contains a bound sodium ion (Na + ), which was also removed. The protein structure was bound with co-crystalline ligand MGP which was extracted, and its data were saved separately. Then, the co-crystalline ligand, water, and ion-free protein structure were saved in a pdb file and imported into software called AutoDock tools 11 . After that, polar hydrogen and Kollman charges were added. The uniform distribution of charges was taken into account, and the protein coordinates were saved in pdbqt file format.

Small-molecule database preparation. The DrugBank (https:// go. drugb ank. com/) database (which is a chemical library) was used for the present work. Two-dimensional (2D) chemical structures of molecules were retrieved from DrugBank in SDF file format. The Open Babel software 12 was employed to convert all 2D chemical structures into 3D structures, followed by energy minimization and structure optimization. An in-house bash script was used to run the Open Babel software. For structure optimization, an MMFF94 force field built into Open Babel was utilized. Using the steepest descent algorithm 13 , each ligand structure was minimized for 10,000 steps. Then, each minimized structure was saved in the pdbqt file format.

Molecular-docking-guided virtual screening. For this purpose, the AutoDock Vina software 14 was used. An in-house bash script was executed to implement the molecular-docking-based virtual screening process. A confined co-crystalline ligand-binding site was regarded as a receptor grid. By means of receptor grid X, Y, and Z coordinates of 2.809, 0.512, and 16.668, respectively, and 20 Å grid box dimensions, an AutoDock Vina configuration file was prepared with an exhaustiveness setting of 8. After successful execution of the dockingbased virtual screening, 10 conformers were retained for each ligand. The PyMol software was employed for visual inspection of the docking results. The Maestro-v12.3 visualization ware (Schrödinger Release 2020-1; Maestro, Schrödinger, LLC, New York, NY, 2021) was used for rendering images.

Docking validation. This validation is a crucial step in molecular-docking-based virtual screening 15 . The 3D structure of the co-crystalline ligand MGP was retrieved from the DrugBank database (https:// go. drugb ank. com/) (DrugBank ID: DB02716) in the SDF format. By the process described in above subsection, this molecule was prepared and saved in the pdbqt format. Next, the prepared molecule DB02716 was docked with the target protein by means of X, Y, and Z grid coordinates of 2.809, 0.512, and 16.668, respectively, and a grid box dimension setting of 20 Å. Ten docked conformations were recorded. Each docked conformation was superimposed upon its native pose (co-crystalline MGP pose), and root mean square deviation (RMSD) was calculated. For superimposition, the RMSD calculation pair_fit plugin script for the PyMol software (https:// pymol wiki. org/ index. php/ Pair_ fit) was used.

MD simulations. These simulations were conducted by means of GPU-accelerated Gromacs 2018.1 16 software. The Charmm36 17 force field was applied to prepare the protein topology. The SwissParam online serverbased software 18 was used to generate the ligand parameterization topology. Each system was solvated via the TIP3P 19 water model resulting in a 10 × 10 × 10 Å cubic box. An adequate amount (0.15 M) of Na + and Cl-ions was added to neutralize each solvated protein-ligand system. The steepest descent algorithm was executed to minimize each system with a maximum of 100,000 steps, and the force was set to the value lesser than 10.0 kJ/ mole. Two-stage equilibration steps were analyzed. In the 1st step, i.e., the NVT ensemble step, volume, temperature, and the number of particles were kept constant and maintained for 2 ns. The 2nd step is the NPT ensemble step with constant pressure along with temperature, and the numbers of particles were equilibrated for 10 ns. For each equilibration step, a 100 ns positional restraint was applied to Cα atoms. To maintain a solvent equilibrium, free movements were allowed for the solvent molecules. The linear constraint solver algorithm 20 was executed to constrain covalent bonds of the system. The particle mesh Ewald 21 method was applied to long-range www.nature.com/scientificreports/ electrostatic interactions with a cutoff of 1.2 nm and a Fourier spacing of 1.2 nm. The V-rescale weak coupling technique was used to regulate the temperature (300.00 K) of the system. The Parrinello − Rahman method 22 was utilized to regulate 1 atm pressure, density, and total energy of the system. Each equilibrated system with acceptable geometry, solPyMOL > pair_fit 6xya_entry_00001_conf_02, ccl-min, and Executive RMS: RMS = 1.301 (33 to 33 atoms) vent orientation was subjected to a 100 ns production run without application of any restraint, followed by a 2 femtosecond (fs) step. The structural coordinates were recorded at 2 picosecond (ps) intervals. After successful completion of an MDS, water and ions were stripped out, followed by Periodic boundary correction (PBC) to refine the trajectories. From the refined trajectories, various parameters, such as RMSD 23 , root mean square fluctuation (RMSF) 24 , radius of gyration (Rg) 25 , solvent-accessible surface area (SASA) 26 , and the number of hydrogen bonds (H-bonds) between the ligand and protein were calculated. The VMD software 27 was utilized to visualize the trajectory and render the images, and the Grace software (https:// plasma-gate. weizm ann. ac. il/ Grace) to plot the data.

Binding-Energy Calculation. The ligand-protein binding interaction was quantitatively estimated by the widely accepted MMPBSA 28 approach. The MMPBSA method for the GROMACS (g_mmpbsa) script program 29 was used to perform the MMPBSA-based binding-energy calculation. To explain the basic working principle of this program, a few known equations are given below. It is a well-known fact that for stability of a protein-ligand complex, the energy of the system must be less than the energy of its individual components. Suppose the term ΔGbind represents the binding free energy of a protein with a ligand. Then, ΔGbind can be expressed as

where Gcom denotes the free energy of the protein-ligand complex. Terms GP and GL represent the free energy of the unbound protein and ligand in a solvent, respectively. Now, the individual components GP and GL can be expressed in Eqs. (2a) and (2b), respectively.

where T and S denote the temperature and entropy, respectively. The term EMMP represents molecular mechanics potential energy in vacuum, and Gsolv is the free energy of solvation. EMMP can be calculated from the molecular mechanics force field parameters via the formula where EB represents bonded interactions, angle dihedrals, and other parameters, and ENB represents nonbonded interactions such as electrostatic (EETS) and van der Waals (EVDW) interactions, respectively. Therefore, Eq. (3) can be written as

The term ENB can be expressed as From Eqs. (2a) and (2b), Gsolv can be calculated by means of an implicit solvent model. Therefore, this term can be expressed via the contribution of electrostatic or polar solvation energy (GP) and nonelectrostatic or non-polar solvation energy (GNP) to the solvation free energy. The g_mmpbsa script along with the APBS software, compute the above terms. To perform the above-mentioned calculation, snapshots of the last 10 ns of a MD trajectory were extracted. The extracted snapshots along with the tpr and index file were employed as input for the g_mmpbsa program to compute the binding energy.

Docking validation. The success of the molecular-docking-based virtual screening validation step is crucial for such projects. In the present work, various co-crystalline ligands were docked with the SFTS L protein.

The docking program generated 10 conformations or docking poses. Each docking pose was superimposed upon the native pose (co-crystalline orientation) of a ligand, MGP. After that, the resultant RMSD (Supplementary Fig. 1a ) was found to be ~ 1.8 Å, and the corresponding docking score was − 6.8 kcal/mol. The co-crystalline ligand MGP in its native orientation showed π-π stacking interactions with Phe1703 and Tyr1719; H-bonding interactions with Gln1707, Asp1771, and Leu1772; and hydrophobic interactions with residues Pro1706, Ile1738, and Ile1774 of the SFTS L protein (Supplementary Fig. 1b) . Protein-ligand interaction analysis of the docking results revealed that the ligand MGP engages in π-π stacking interactions with Phe1703 and Tyr1719 and H-bonding interactions with Phe1703, Gln1707, Asp1771, and Trp1725. Additionally, hydrophobic interactions with Pro1706, Ile1738, Ile1774, Leu1768, and Leu1772 were found (Supplementary Fig. 1c) . The RMSD and interaction analysis strongly supported the validity of the docking protocol 30 . Moreover, based upon the above observations, a strategy was chosen for virtual-screening-based hit identification vis-à-vis the SFTS L protein. Any small-molecule drug that manifested a Vina docking score of more than −7.0 kcal/mol and interactions with Phe1703, Tyr1719, Gln1707, Asp1771, Pro 1706, Ile1738, and Ile1774 was shortlisted and included in a rank. The π-π-type interaction was assumed to be one of the important interactions for ligand binding. The above hypothesis was tested to compile a short list of small-molecules from the DrugBank database. www.nature.com/scientificreports/

In the present study, Autodock Vina was employed as docking software. An in-house bash script was used to execute the docking program to screen out many molecules. Initially, the top 10 molecules were chosen (Table 1 ) based on a highly negative docking score. The top ligands, i.e., bromfenac, elliptinium, cinchophen, and zaltoprofen, had Vina docking scores of −8.7, −8.5, −8.2, and −8.0 kcal/mol respectively. Ligand oxcarbazepine showed a docking score −8.20 kcal/mol, but clinically, it is used as an anticonvulsant agent. Therefore, this molecule was not considered further in the present study. Other ligands such as cyproheptadine, epinastine, mianserin, midazolam, and phenytoin had docking scores greater than −8.0 kcal/mol. Bromfenac featured the lowest docking score among all the ligands subjected to the virtual screening (Table 1 ). According to the interactions with relevant amino acid residues of the SFTS L protein (Fig. 1) , the best four molecules were selected. Figure 1a shows the interaction of the co-crystalline ligand with the SFTS L protein. The ligand bromfenac shows π-π interactions with Phe1703 and Tyr1719, an H-bonding interaction with Gln1707, and hydrophobic interactions with Pro1706, Asp1771, and Leu1772 of this protein (Fig. 1b) . The ligand cinchophen engages in π-π interactions with Tyr1719 and Phe1703, H-bonding interactions with Leu1772, and hydrophobic interactions with Asp1771, Ile1774, Ile1738, and Pro1706 of the SFTS L protein (Fig. 1c) . The ligand elliptinium shows π-π interactions with Tyr1719 and Phe1703, an H-bonding interaction with Leu1772, and hydrophobic interactions with Ile1774 and Pro1706 (Fig. 1d) . The ligand zaltoprofen engages in π-π interactions with Phe1703 and Tyr1719, an H-bonding interaction with Phe1703, and hydrophobic interactions with Leu1772, Ile1774, and Pro1706 residues (Fig. 1e) .

In the docking analysis, it was found that the ligands bromfenac, cinchophen, elliptinium, and zaltoprofen were engaged in 2 or more π-π interactions. For this reason, these four ligand-protein systems were analyzed further, by MD simulations. Because these are clinically established small-molecule drugs, no theoretical absorption, distribution, metabolism, and excretion analyses were performed on them here. Table 2 . RMSD values of each frame were calculated from the total 100 ns MD trajectory and plotted against time (Fig. 2) . The RMSD parameter calculated from the MD trajectory denotes the probable changes that transpire in the geometric orientation of a protein structure. RMSD less than 3 Å indicates superior stability of a globular protein 23, 31 . From the average RMSDs in Table 2 , it can be concluded that each ligand oscillates at lower RMSD in comparison with the Apo protein (RMSD ~ 1.63 Å) except for the elliptinium-protein system.

Among the 4 protein-ligand systems, the bromfenac-protein system has the highest maximum RMSD. The zaltoprofen-protein system has the lowest RMSD profile among the four systems. Each system oscillates at an RMSD of less than 3.0 Å. An in-depth analysis of Fig. 2 indicates that none of the identified candidate systems undergo abrupt changes in the RMSD profile.

To understand the behavior of individual amino acid residues during the 100 ns MD simulation, RMSF properties were calculated next, and the average values can be found in Table 2 . Each system other than the bromfenac system possesses a lower average RMSF than the Apo protein's RMSF. The bromfenac-protein system has a 0.01 Å higher average RMSF profile than the Apo protein. Considering that MD is a stochastic process, this difference is acceptable 32, 33 . For detailed examination of the RMSF profile, RMSF of each amino acid residue was calculated from the 100 ns MD trajectory and plotted against the residue number (Fig. 3) . Phe1703, Pro1706, Gln1707, Thr1709, Tyr1719, Asp1771, and Leu1772 are the crucial amino acid residues in the ligand-binding site of the SFTS L protein, and in-depth analysis of Fig. 3 revealed that no significant changes occur when one of the 4 ligands binds to this protein.

The radius of gyration (Rg) calculated from the MD trajectory indicates the compactness or rigidity of a protein system during the simulation 25 . A consistent Rg profile indicates that the protein system may be stable and does not undergo any significant structural changes or distortions during the simulation. Average values of Rg calculated for each system can be found in Table 2 . To understand changes of Rg with time, a plot was constructed (Fig. 4) .

In Table 2 , readers can see that each system has the average Rg either similar to or less than the Apo protein's Rg of ~ 13.97 Å. Analysis of the Rg-versus-time plot constructed for each system indicates that there are no www.nature.com/scientificreports/ significant changes in the Rg profile. In accordance with the above observations, it can be stated that in each system, the SFTS L protein retains its geometric orientation and forms a stable complex with each of the identified promising ligands. SASA parameter calculated from the MD trajectory provides insight into the events in the protein-ligand contact. If, ligand binding to a protein is considered as a solvent replacement approach, then lower SASA value means that the protein or ligand is less exposed to the solvent. That is, the ligand is located deep inside the binding pocket 26 . According to Table 2 , SASA values of each ligand-protein system-as calculated from the 100 ns atomistic MD trajectory-are less than SASA of the ligand-free protein. This critical observation strongly indicates that during the 100 ns MD simulation, each ligand is situated inside the receptor-binding pocket. To understand the changes of SASA with time, a plot was constructed (Fig. 5) .

The binding energy in terms of MMPBSA was also calculated from the 100 ns MD trajectory. Highly negative binding energy implies stable binding of a ligand to a protein 15 . The average values of the binding energy in kJ/mol are shown in Table 2 . The binding energy of co-crystalline ligand MGP was calculated from molecular dynamics trajectory and compared with the binding energy of identified binders. The zaltoprofen-protein system has a highly negative binding energy −220.095 kJ/mol. The elliptinium-protein system features a moderate binding energy −101.738 kJ/mol. The bromfenac and cinchophen systems show binding energies of −183.478 and −145.806 kJ/mol, respectively. In addition, from the Table 2 it can be observed that the binding energy of each identified hits (except Elliptinium) were found to be high negative in compare to MGP bound system. The MMPBSA-based binding-energy profile of the identified ligands toward the SFTS L protein suggests that each ligand binds stably. This statement is further supported by data on other parameters, such as RMSD, RMSF, Rg, and SASA, calculated from the MD trajectory.

In conclusion, we report the identification of small-molecules (from the DrugBank database) that may stably bind to the SFTS L protein. The Autodock Vina docking software was used to implement the virtual-screening process. By the molecular-docking-based virtual screening, 10 small-molecule hits were identified on the basis In conclusion, this study presents four clinically used small-molecules exhibiting stable binding with SFTS L protein for up to 100 ns of MD simulation time. Nonetheless, proper biological evaluation is necessary to determine their efficacy against SFTS.

All data will be available upon request to the corresponding author. 

Fever with thrombocytopenia associated with a novel bunyavirus in China

Severe fever with thrombocytopenia syndrome

Prevalence of severe fever with thrombocytopenia syndrome virus in Haemaphysalis longicornis ticks in South Korea

Epidemiological and clinical characteristics of severe fever with thrombocytopenia syndrome (SFTS) in China: an integrated data analysis

Epidemiological and clinical features of laboratory-diagnosed severe fever with thrombocytopenia syndrome in China, 2011-17: a prospective observational study

Severe fever with thrombocytopenia syndrome, an emerging tick-borne zoonosis

A multicenter non-randomized, uncontrolled single arm trial for evaluation of the efficacy and the safety of the treatment with favipiravir for patients with severe fever with thrombocytopenia syndrome

Person-to-person transmission of severe fever with thrombocytopenia syndrome virus. Vector Borne Zoonot

Efficacy of T-705 (Favipiravir) in the treatment of infections with lethal severe fever with thrombocytopenia syndrome virus

Structural and functional characterization of the severe fever with thrombocytopenia syndrome virus L protein

The nucleoprotein of severe fever with thrombocytopenia syndrome virus processes a stable hexameric ring to facilitate RNA encapsidation

The Cap-Snatching Mechanism of Bunyaviruses

Molecular docking and simulation investigation: effect of beta-sesquiphellandrene with ionic integration on SARS-CoV2 and SFTS viruses

PDBsum: Structural summaries of PDB entries

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

Open Babel: An open chemical toolbox

Efficiently finding the minimum free energy path from steepest descent path

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

Identification of potential inhibitors of SARS-CoV-2 main protease and spike receptor from 10 important spices through structure-based virtual screening and molecular dynamic study

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers

CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data

SwissParam: A fast force field generation tool for small organic molecules

Comparison of simple potential functions for simulating liquid water

Insights into the structure-function relationship of brown plant hopper resistance protein, Bph14 of rice plant: a computational structural biology approach

All-atom continuous constant pH molecular dynamics with particle mesh ewald and titratable water

Polymorphic transitions in single crystals: A new molecular dynamics method

Methods of protein structure comparison

Automatic identification of mobile and rigid substructures in molecular dynamics simulations and fractional structural fluctuation analysis

Radius of gyration is indicator of compactness of protein structure

Fast pairwise approximation of solvent accessible surface area for implicit solvent simulations of proteins on CPUs and GPUs

VMD: visual molecular dynamics

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

Open Source Drug Discovery Consortium, Lynn A. g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

In silico analysis and identification of promising hits against 2019 novel coronavirus 3C-like main protease enzyme

Stochastic classical molecular dynamics coupled to functional density theory: Applications to large molecular systems

The authors declare no competing interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-00294-7.Correspondence and requests for materials should be addressed to D.-M.K.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.