Molecular dynamics simulation approach for discovering potential inhibitors against SARS-CoV-2: A structural review

Since the commencement of the novel Coronavirus, the disease has quickly turned into a worldwide crisis so that there has been growing attention in discovering possible hit compounds for tackling this pandemic. Discovering standard treatment strategies is a serious challenge because little information is available about this emerged infectious virus. Regarding the high impact of time, applying computational procedures to choose promising drugs from a catalog of licensed medications provides a precious chance for combat against the life-threatening disorder of COVID-19. Molecular dynamics (MD) simulation is a promising approach for assessing the binding affinity of ligand-receptor as well as observing the conformational trajectory of docked complexes over time. Given that many computational studies are performed using MD along with the molecular docking on various candidates as antiviral inhibitors of COVID-19 protease, there is a demand to conduct a comprehensive review of the most important studies to reveal and compare the potential introduced agents that this study covers this defect. In this context, the present review intends to prepare an overview of these studies by considering RMSD, RMSF, radius of gyration, binding free energy, and Solvent-Accessible Surface Area (SASA) as effective parameters for evaluation. The outcomes will offer a road map for adjusting antiviral inhibitors, which can facilitate the selection and development of drug candidates for use in the medical therapy. Finally, the molecular modeling approaches rendered by this study may be valuable for future computational studies.


Introduction
Since the fall of 2019, a new respiratory infection has appeared in the world, SARS-CoV-2 was identified as a recently realized infectious virus and the main reason for the Coronavirus pandemic of 2019, which has been denominated COVID-19 by WHO (World Health Organization). A severe disease was globally classified as a disaster epidemic. The situation has worsened due to a shortage of prescription treatments to combat this novel Coronavirus. Since this is a novel virus, details on its biological activity, molecular etiology, and medical science are recently starting to appear. Being from the RNA virus of Coronaviridae group, Coronavirus is comprised of four types of alpha (a), beta (b), gamma (c), and delta (d) in which the b (b-CoV) kind includes intense respiratory syndromes of Coronavirus (SARS-CoV), (SARSCoV-2) and (MERS-CoV). SARS-CoV-2, which causes COVID-19, is much more contagious than SARS-CoV and MERS-CoV, which spreads from person to person and inflicts mortal disease [1]. Fever, nausea, dry cough, diarrhea, and shortness of breath are some of the signs of novel Coronavirus [2][3][4] that may escalate to acute respiratory diseases. The genome of SARS-CoV-2 encrypts structural molecules including S (spike glycoprotein), E (envelope), M (membrane), and N (nucleocapsid), as well as non-structural components as with M pro , papain-like protease (PLP), and RNA dependent RNA polymerase (RdRp). The transcription and replication cycle occur due to the non-structural proteins, while structural proteins are mainly liable for interactions with human organisms during the entrance process of the virus [5][6][7][8]. Viral RdRp is responsible for speeding up the replication of RNA which 97.08 % of its sequence has been partaken with SARS-CoV [9]. Also, PLP's potential to antagonize interferon (IFN) interaction and deubiquitinate viral and cellular proteins have a critical portion for infection progression [10]. On the other hand, as a crucial enzyme, the main protease (M pro )/3-CL pro is responsible for virus replication so that was considered a principal therapeutic target in the fight against SARS-CoV. The structure of SARS-CoV-2 main protease is composed of domain I (residues 8-101) and II (residues 102-184) including b-barrel, which a long loop (residues 185-200) attaches them, as well as domain III (residues 201-303), composes a-helices. The catalytic dyad in the active site of the M pro is composed of the His41-Cys145 residues, which are placed between domains I and II [11]. COVID-19 is a highly infectious virus for which there is no perfect remedy [12]. With the COVID-19 pandemic spreading at an increasing pace and the death rate rising, researchers from all over the world have attempted feverishly to devise and detect potential therapeutic agents against viral main protein targets to combat Coronavirus outbreaks. Along with many clinical trials, computational methods were considered as a spotlight for the assessment and comparison of available drugs with the ability to affect SARS-CoV2 via various pathways. Computational approaches are used to examine highly selective inhibitors against critical viral targets, providing the principally detailed foundation for anti-viral drug development. The first step in the antiviral inhibitor detection procedure is to identify possible inhibitors by logical screening. The next step entails using molecular dynamics simulation to deduce the physical and biological structure of the inhibitor binding process. The molecular dynamics method is a procedure that estimates the behavior history of atoms so that visibility of the complete system's position and velocity variations was provided by allowing atoms and molecules to interact atomically in the MD process. The molecular dynamics simulation was applied in some studies for perusing the atomic behavior of Coronavirus in different situations. Malekahmadi et al. [13] evaluated the stability and atomic behavior of Coronavirus under various conditions of temperature and pressure using MD approach. The Coronavirus was indicated by S, O, N, and C atoms and the results have been presented with calculating physical characteristics of temperature, total energy, volume variation, and atomic force. As a key factor in identifying the structural stability, the volume of Coronavirus is revealed to raise 92% and 14% by 100 K and 2 bar changes in the temperature and pressure, respectively. Moreover, as temperature and pressure rise, Coronavirus potential energy in the aqueous atmosphere reduces while the volume of the Coronavirus and consequently, instability increases. It is concluded that the variations of thermodynamic properties like temperature and pressure affect the atomic treatment of Coronavirus. Karimipour et al. [14] characterized the atomic behavior of Coronavirus during interacting with a different metallic matrix of Fe, Al, and steel using the molecular dynamics method. For simulation, the Coronavirus was demonstrated with S, O, N, and C atoms and some properties including potential energy, temperature, the volume change of simulated virus, the center of mass distance, and angle have been investigated. According to obtained results, the interaction with the steel matrix has allowed the virus to be eliminated from the surfaces to the greatest extent possible as well as the volume of Coronavirus has raised by 14.62% during interacting with the steel matrix, which causes the steel matrix could be considered as a material with antiviral properties in medical usage. Furthermore, it is reported that as the temperature rises, the repulsion between the matrixes and virus increases.
Kumar et al. [15] applied five various algorithms to prediction of high-scoring CTL epitopes. Molecular docking along with MD simulations were utilized to assessment of compatibility vaccine leads with objective receptors. The vaccine design was well-specified by physicochemical properties as remarkable result. Chatterjee et al [16] investigated the binding stabilization of CQ-TrOne complex with M pro using molecular dynamics simulation by analyzing RMSD, RMSF and R g parameters. Furthermore, the binding energy diagram displayed the effective binding attachment towards M pro accordance to essential reaction coordinate results. Their major findings revealed high capacity of considered compounds, which can be evaluated using in vitro and in vivo procedures to produce the efficient COVID-19 treatment.
Kumar et al. [17] evaluated the Noscapine Hydroxychloroquine (Nos-Hcq) conjugate as possible antiviral candidate using MD method. The RMSD, radius of gyration as well as the crucial reaction coordinate binding free energy plots, confirmed the binding stability of conjugation during simulation time. Furthermore, they demonstrated that agents with more stability to key domains of M pro could considerably boost the reaction coordinates, medicine access, and inhibitory regulation.
In another study [18], the MD simulations indicated the robust binding of vaccine for HHV-5 with MHC receptors as well as virus particular membrane receptor TLR2. Through RMSD, RMSF, and secondary structure confinement examinations, the interaction route study of the vaccine revealed steady binding with small devi-ations. Kumar et al. [19] adopted an immunoinformatic platform to develop a reliable vaccine as a future healthcare plan. They conducted MD simulation to study the stable binding with TLR-2. According to obtained results, a viable vaccination strategy was proposed to stimulate an immune replication to combat coronavirus.
Respecting the strong points of the computational methods, better knowledge of existing agents as anti-virus drugs can be used to design and find novel effective inhibitors for the treatment of COVID-19. The replication mechanism of virus can be successfully stopped by inhibiting the activity of SARS-CoV-2 main protease. In this regard, virtual screening of potential compounds derived from available databases accompanied by the MD approach was currently utilized in several types of research to characterize putative inhibitor activities. So, this review aims to cover the studies whose line of research has been concentrated on analyzing the molecular structure and stability of possible inhibitors against known targeted SARS-CoV-2 protein from the aspect of molecular dynamics simulation. Pioneer English publications that have used molecular dynamics approach for investigation of promising COVID-19 drugs with focusing on the prediction and determination of protein-ligand binding structures are included in this comparative evaluation. By considering the above point of reference, title and abstract searching have been performed according to the specific phrases and keywords contained but not restricted to the below phrases: ''Novel Coronavirus", ''Inhibitors of Coronavirus", ''MD simulation of antiviral agents", ''SARS-CoV-2", ''virtual screening of COVID-19 hit compounds", ''Drug Therapy". In general, this review has intended to provide a concise framework to help classify and identify drugs that have a high score in performance when binding to main proteases of Coronavirus.

MD simulation
MD simulation is a computational procedure for observing the physical interactions in a biophysical environment in which the structural changes and flexibility of docked complexes can be visualized during the simulation time [20][21][22][23][24][25]. In this approach microstate parameters including a number of particles (N), the volume of the system (V), and total energy (E) which are determined as NVE ensemble [26], are affected by macroscopic variables that should be run with applying a thermostat and barostat in the simulated system. MD simulation can be exploited to visualize the real movement and structural modifications of a protein in a biological system. MD trajectories can be evaluated by calculation of Root-Mean-Square Deviation (RMSD) and the Root Mean Square Fluctuation (RMSF) of the compounds.

Selection of protein structure
Various researcher groups have already been focusing on the discovery of potential medicines or pharmaceuticals to cure infections or symptoms. Until now, by X-ray crystallography and cryo-EM tools, various crystal structures of SARS-CoV-2 proteins have been determined and stored in the Protein Data Bank (PDB) (https://www.rcsb.org/) with specified PDB IDs. For the current review, different discussed studies have utilized the initial configuration of target protein in complex with objective inhibitors from sample prepared by PDB as the 3D crystal structures of different SARS-CoV-2 proteins which are presented in discussion. Afterwards, the samples are imported into studio program so that their structural sequence can be aligned and evaluated for comparison aims and further investigations.

Docking studies
One of the significant materials in silico drug design and detection procedure, is molecular docking. It determines the interaction between the molecule and receptor using binding affinity score. Structure-based virtual screening relying on molecular docking approach is leveraged to find possible hit drugs against main proteases of COVID-19 [38]. To examine pathway and interactions between receptor and drug molecules, MD analysis of leading docked complexes are done [41]. The 3D configuration of the main protease protein is considered during the docking study. Prior proceeding the docking assay, the non-required water molecules are removed [32]. It is possible to adjust both individual and sequential docking techniques. Due to the medication repurposing study, the ligand-binding sites may differ from the protein's traditional ones, which are normally where anti-Covid medicines bind. As a result, initially a blind docking can be run in which the entire protein is placed in the grid box and potential binding sites are explored to identify best binding sites (BBS). After that, sequential docking can be done in which each ligand is docked to its appropriate BBS in order. Using this approach, the synergistic or antagonistic influence between the ligands can be prospected. Individual and sequential molecular docking combining MD simulation, RMSD, and RMSF evaluation can discover promising ligands to COVID-19 treatment [33]. Docking studies has been performed in some studies using Workflow (VSW) of Glide Schrodinger Suite [34,36] and on others with Auto dock Vina Wizard of PyRx software [37,41,42].

Protein structure energy minimization and optimization
In docking and simulation investigations, energy minimization and optimization are critical stages. After deriving the 3D structure of SARS-CoV-2 main protease associated with an selected candidates from the Protein Data Bank, the whole structure is optimized and minimized using Protein Preparation Wizard in the Schrodinger suite as a standard practice. This includes removing crystallographic waters, replacing vacant hydrogens to allocate suitable charge and protonation state. Further, the energy minimization is carried out considering a root mean square deviation (RMSD) cut-off value of 0.30 Å using OPLS3 (Optimized Potentials for Liquid Simulations) force field [34][35][36]38].

RMSD analyses
The RMSD parameter which is the rate of the mean distance between the atoms is required for examining the equilibration and the structural stability of the protease in the presence of a docked ligand. The RMSD value of a solitary protein and a ligand-docked one can be compared to evaluate changes in protein molecular dynamics as well as the conformational stability of the protein-ligand complex. The RMSD of theatomic coordinates between the initial position at (t = 0) and the new atomic sites at specified time t can be employed to analyze system conformations through determining the resemblance in three-dimensional structure [46]. The structural stabilization of the complex is demonstrated by a lower RMSD value of ligand-docked protein than solitary one. This boosts the efficacy of chosen ligands as possible COVID-19 therapeutics [33]. The variations in the RMSD values over time suggests a possible conformational deviation in the enzyme structure as a result of inhibitor attachment. The higher the RMSD for one or a group of atoms during the simulation, the greater their structural [28].
The following equation has been presented to compute RMSD quantity: where x m , y m ; z m and x 1 ; y 1 ; z 1 represent the initial coordinates and trajectory coordinates at frame t, respectively. N is the number of atoms.

RMSF analyses
The RMSF can be used to survey the flexibility of the protease in the involvement a compound as well as for identifying local variations in the protein chain [27]. The variation of flexibility in terms of RMSF parameter can be utilized to investigate inhibitor binding to the target. RMSF values for agent atoms would be obtained from trajectories recorded during MD simulation runs to determine the stiffness and flexibility of residues in COVID-19 protease following binding of the identified drugs [34]. The RMSF plot would be a typical representation for residues that have suffered significant changes throughout the MD simulation operation [33]. The peaks in the RMSF diagram depict the residues with greatest oscillation during simulation [37]. Also, higher RMSF values mean that the protein has more flexible domains. This parameter is calculated through the given equation as follow: where T demonstrates trajectory frame numbers and x À is the timeaveraged position.

H-bond interaction analyses
The H bond is constructed by covalent bonding of the hydrogen atom to an electronegative atom with another one. Hydrogen bonds between the main chain and the amide nitrogen stabilize secondary structures. Hydrogen bonding are also linked to the protein's structural stiffness [46]. HBs have a critical role in maintaining protein stability. The hydrogen bond is critical in providing a stable foundation for biological systems. As a result, the MD simulation trajectories may be used to compute the hydrogen bonds in drug-protein combinations [39]. This analysis delves more descriptions into the enzyme-ligand binding process with particular consideration [34]. By analyzing positions at each time step and measuring the proportion of available H-bond throughout the simulation, the molecular predictors of H-bond are determined [41]. Using MD scheme, the average values for Hydrogen bonds between the main protease and respective inhibitors are estimated for whole trajectories.

Binding free energy analysis
The analysis of binding energy reveals significant residues for inhibitor design. The structural stabilization of inhibitor in the catalytic region is reflected by thermodynamic energy proportion to the total binding free energy of the complex. The binding affinity, stability, and selectivity of the inhibitor are all influenced by interactions in the residues. To determine the influence of the evaluated agents, it is required to analyze the binding affinity of all of the chosen hits towards main protease [34]. The ligand-receptor complex binding free energy can be applied to assess diverse conformations. To grade conformations derived by molecular docking, the scoring tool is based on free energy calculations [37]. In fact, higher rate of free energy corresponds to more desirable binding affinity with the protein higher energy [39]. Although the virtual screening approach present highest-scoring hit compounds, dock-ing is unable to quantify the ligand-receptor complex's binding free energy. As a result, accurate techniques are required to estimate the receptor-ligand complex's binding free energy. In this regard, through molecular dynamics simulation the conformational variation according to free energy distribution assessment upon ligand binding can be anticipated [32]. To investigate the binding free energy of the top scored agents, Molecular Mechanics with Generalized Born and surface area solvation (MM-GBSA) and Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) methods can be used, which are recognized strategies for predicting binding energy and identified to be more reliable than other methods [29,32]. The procedure is based on MD simulation of complexes that has been tried on a variety of systems with varied degrees of success. Using routs of MD processes the comparative enthalpy variation for the generation of the complex over MD simulation times is obtained [43]. In MM-GBSA the binding free energy of a ligand (L) to the receptor (R) organizing a complex (RL) is calculated as below: In above equationsDG binding ,DE vdW and DE ele are related to bonded energy, Van der Waals, and electrostatic portions, respectively. Also, DG GB and DG SA are involved in the free energy of solvation which is computed from generalized Born equation and Solvent-Accessible Surface Area (SASA), respectively.

Radius of gyration (Rg) analysis
Another parameter is radius of gyration Rg which is considered as a fundamental indicator of the total size of chain molecule. It can be employed to measure how much a protein's structure varies during MD simulations. The Rg evaluates the compactness manner and flexibility of the protein inside a biological environment so that compares the structure of the protein per time to the hydrodynamic radius that can be monitored experimentally. Respect to the MD approach for evaluated SARS-CoV-2 main protease, Rg demonstrate resembling structural compactness of comparable candidates [44]. The lower values of Rg describe a more rigid structure during the simulation. In this regard, the structural change of complexes can be quantified using the Rg that is stated according to the below formula [30]: where rðiÞ displays the coordinates of the atom i and r center is the center of mass. N refers to the number of protein atoms.

Protein Solvent Accessible Surface Area (SASA)
Protein Solvent Accessible Surface Area (SASA) has been regarded as a critical element in protein folding and stability research, which is known as the surface identified around a protein by the Van der Waals contact surface of the molecule and a hypothetical center of a solvent sphere. Examining the compactness treatment of backbone atoms can be assessed by the SASA parameter. Over time, the changes in SASA quantities illustrate the folding and unfolding of complexes. It's mentioned that a lower solvent accessible surface area value indicates greater compactness. From the atomic coordinates, the SASA of a native protein can be calculated numerically. Also, there are many simulation approaches for calculating the unfolded form of a SASA for a protein. A variety of unfolded state layouts have been presented in an attempt to approximate the variations in SASA relevant to protein folding [31].

Discussion
One of the difficulties in the controlling of COVID-19 is the development of inhibitors as effective antiviral drugs. Considering the above challenge, many studies support the use of MD simulation to investigate conformational variations and flexibility of different combinations with the main protease of SARS-CoV-2. For this effect, an MD simulation has been executed to peruse conformational changes of the docked complexes including lopinavir, darunavir, hydroxychloroquine, interferon-alpha, remdesivir, niclosamide, ribavirin, umifenovir, ritonavir, ivermectin, and as well as phytochemicals such as isoobtusitin, ellagic acid, apigenin, carnosic acid, drymaritin, morin, scutellarein and triterpenoids with the main protease of SARS-CoV-2 [32]. To examine the dynamical characteristics and binding effectiveness of the mentioned complexes, MD simulations were done using the Gromacs software. For the 100 ns trajectory, the simulation process was run with the all-atom force field CHARMM 36. AutoDock Vina was employed to analyze the phytochemicals and medicines affinity for the main target. Due to the great importance of binding energy between host and receptor, an evaluation was performed on the bonding behavior of the mentioned complexes. The results of MD simulation for 100 NS represent the values of À98.858 and À47.326 kJ/mol for the binding energy of remdesivir and carnosic acid, respectively. Among investigated inhibitors, the binding tendency of remdesivir and carnosic acid into the main protease has demonstrated superior results. A survey of interaction energy contribution suggests that Van der Waals energy is more prominent in comparison to others, which detects effective hydrophobic contact of remdesivir and carnosic acid with the receptor. Assessment of RMSD outputs for SARS-CoV-2 main protease and listed complexes demonstrates the stability of the combinations. Based on the obtained trend of oscillations for the main protease, carnosic acid and remdesivir docked complexes; it was observed that all three components reach stability after 60 ns. On the other hand, the results of RMSF describe the stability behavior of protein so that fewer oscillations have been perceived for SARS-CoV-2 main protease backbone atoms compared to the combination with remdesivir and carnosic acid. From the radius of gyration calculations for determining the treatment of protein compactness, the higher value of Rg is comprehended for remdesivir-docked complex and as a result, has more flexibility than carnosic acid docked complex and main protease. Prediction of the number of hydrogen bonds during the simulation time explains that up to 40 ns the number of hydrogen bonds for a complex of carnosic acid is more than two others, while in the rest of the simulation time the highest number of hydrogen bonds was related to the main protease.
From the sight of another strategy, composed drug therapy was investigated as an anti-covid therapeutic facility by molecular dynamics simulation. Dextromethorphan, Prednisolone, and Dexamethasone were utilized as ligands in combination with receptor of COVID-19 main protease (M pro ) [33]. In the first step, a complex of SARS-CoV-2 protease (PDB Id: 6LU7) was made with Dexamethasone named 6LU7-D1. Next step, Prednisolone was docked with 6LU7-D1 and finally, the created complex was combined with Dextromethorphan and the docked complex of 6LU7-D1-P-D2 was formed for analysis. The conformational behavior of considered ligand-receptor docking was inspected using MD simulation for 100 ns time scale so that the atomic trajectory for surveying structural stability has been assessed by computing RMSD and RMSF properties. GROMACS software was used to create the topology files and ran MD simulations on the 6LU7-D1-P-D2 complex using the Gromacs96 53a6 force field. It has been found that the synergistic impact of these combined agents produces greater anti-Covid property than their separate effects. The MD results of RMSD computation for main protein (6LU7) as well as the protein-ligand complex (6LU7-D1-P-D2) present the numerical range of 1.5 to 3.1 for 6LU7 and 1.0 to 2.1 for 6LU7-D1-P-D2. Since lower RMSD values suggest greater structural stability, it is concluded that the proposed combined ligands with target receptors can be considered as anti-covid therapies. Measurement of RMSF versus residue number for 6LU7 and 6LU7-D1-P-D2 has displayed a similar fluctuation scheme so that the MD simulation has confirmed that there is no significant change in the 6LU7-D1-P-D2 complex pattern. Given this, RMSF examination revealed the conformational stability of simultaneous docking of propounded ligands to the M pro protease.
Another approach to combat COVID-19 has been examined in which peptide-like inhibitors that can be impressible against SARS-CoV-2 infection have been analyzed using MD simulation [34]. For this purpose, the main protease of SARS-CoV-2 (M pro ) was complexed with five peptide-like N3-ILP compounds as potential elements with PubChem IDs 54456426, 54152887, 54035018, 91366909, 57076946 considering 80% resemblance to basic N3ligand. Choosing these compounds has been done relying on the binding scores, lowest RMSD, and binding interactions parameters. The 100 ns MD simulation has been carried out for all proposed ligand-protein complexes. The energy minimization was accomplished using the OPLS-2005 force-field. The AMBER 18 package was utilized to run Partial Mesh Ewald Molecular Dynamics (PMEMD) modeling on all of the produced complexes. The general AMBER force fields (GAFF) were used to parametrize the enzyme and the proposed ligands. Since changes in the structure of the enzyme affect its behavior and performance, it is essential to study the structural dynamics and conformational changes during the simulation. In this regard, RMSD calculation was performed for determining the structural deviation of the enzyme in coupling with ligands. It is revealed that the lowest average RMSD with the value of 2.03 Å has been associated with 54035018-M pro , while this parameter has shown the values of 2.45 Å and 2.33 Å for 54152887-M pro and 54456426-M pro , respectively. Considering the outcomes, it is proposed that binding of 54035018-M pro , 54152887-M pro , and 54456426-M pro creates more stability in comparison with N3-ILP-M pro linkage. In addition, the calculations of Rg to confirm the results of conformational changes represent the values of 38.87 Å, 40.39 Å, and 43.68 Å for 54035018-M pro , 54152887-M pro , and 54456426-M pro , respectively, which indicate small variation in the compactness of these protein-ligand complexes. In addition, among them, the most stability has belonged to the binding of 54035018-M pro . From RMSF computations to examine the rigidity and flexibility of protein-ligand bindings, it is reported that the minimum oscillations in the residues are seen in the 54035018-M pro complex with the value of 11.12 Å. Furthermore, the RMSF result for N3-ILP-M pro indicates a slightly higher value than 54035018-M pro , which has suggested 54035018 as a preferable inhibitor of M pro enzyme. Also, the binding energy assessment of the considered compounds has affirmed that the binding energy of 54035018-M pro (À37.40 kcal/mol) is greater than N3-ILP-M pro (À30.89 kcal/mol), and finally, 54035018 has been recommended as a superior lead for the COVID-19 therapeutic.
As an absorbing inhibitor against Coronavirus, guanine-N7 methyltransferase has been investigated with molecular dynamics simulation to detect the molecular binding behavior of Guanosine-P3-Adenosine-5 0 ,5 0 -Triphosphate (G3A) with C-terminal N7-MTase domain of nsp14 from SARS-CoV-2 [35]. The simulation has been carried out for G3A complex nsp14 and five ligand-bound complexes including TCM 57025 (5a), TCM 3495 (5b), TCM 20111 (5c), TCM 31007 (5d), and TCM 5376 (5e) over a timescale of 50 ns to analyze the interaction comportment. The GROMACS 5.1.4 is employed to execute molecular dynamics. The dynamic behavior and interaction process of substrate-bound complexes as well as ligand-bound ones are modelled for 50 ns so that OPLS-AA force field is used to create the top ligand-bound structures. Computing the RMSD values has exposed the stability of these examined compounds during the simulation. In particular, among them, the value of RMSD for TCM 3495 and TCM 20111 has been reported higher than 0.94 nm which was the highest value compared to other complexes. In complementary for RMSD oscillations the RMSF values have been analyzed for substrate and ligand complexes. It is indicated that in all complexes, the deviation has taken place only in the N-terminal and C-terminal residues which proves the stability of RMSD. On the other hand, investigation of hydrogen-bonding interactions has disclosed the mean hydrogen bond of 3.6 for substrate G3A bound complex with ARS-CoV-2N7-MTase while this parameter for TCM 57025, TCM 3495, TCM 5376, TCM 20111, and TCM 31007 was 3.40, 3.35, 4.00, 3.97 and 3.34, respectively. The PCA outputs for probing the binding of the substrate-ligand complex have demonstrated that all complexes except the TCM 20111 have had acceptable compactness with high consistency of protein-ligand interactions. On the whole, it has been claimed that these five TCM components have a promising possibility as antiviral phytochemicals that could suppress SARS-CoV-2N7-MTase.
In light of the key protease of M pro /3CL pro found in CoVs as appealing COVID-19 targeted therapies and the drug pool's tracking, some potential drugs have been screened in which conceivable molecules including Pepstatin A, Leupeptin Hemisulphate, Nelfinavir, Lypression, Octreotide and Birinapant, have been examined as possible M pro impediments via MD simulation [36]. In this case, firstly the simulation has been performed for finding a competent site in which APO of M pro (PDB-ID: 6M03), as well as COM1 and COM2, have been assessed during 200 ns of MD simulation. The SiteMap tool from Schrodinger Suite was applied to compute binding sites. The OPLS-2005 force field was utilized, with 15 site points per reported site cropped. The OPLS3 force field was also employed to create the systems with the Systems builder. The site of bound peptide 'N3' has been considered as the main area of ligand coupling. From RMSD observations, it has been declared that both APO and COM proteins have been stable across the trajectory with only slight deviations. Also, The RMSF correlations have reported that N3 binding creates M pro protein more pliable, mainly in the region III relative to13b. Furthermore, the fluctuations in this sector in the existence of ligands are due to the involvement of region III in the creation of a homodimer. In comparison to the cocrystallized ligands N3 and 13b, it is found that six factored substances have smaller binding free energy as well as higher specific therapeutics value. The MD evaluation of mentioned compounds confirms that Nelfinavir and Birinapant have provided the most protein stability during simulation in comparison to the other four chosen items. On the other hand, as it was evident in the outcomes the other hits have been discovered to be stable. Performing binding energy analysis on the trajectory planning of the hits and support molecules has confirmed that apart from Leupeptin Hemisulfate all other hits have had significantly higher binding energy than the support molecules. Based on the total investigation, the discovered groups of molecules were picked up as primary hit molecules with appearing optimistic antiviral effect.
In terms of attention to find possible small molecules to treat the 2019-nCoV virus, in the study of Tripathi et al. [37] the efficacy of the best five ligands against SARS-CoV2 including MolPort-002-530-156, MolPort-002-701-723, MolPort-000-410-348, MolPort-027-852-450, MolPort-039-338-091 has been evaluated using the MD approach. 3CL pro has been considered as a potential target for inhibitor developments during 5 ns molecular dynamics simulation to measure the binding free energy. The pmemd module in Amber 20 was employed to conduct the MD assessment. The parameter and coordinate files for target complexes were generated with tleap through amber ff14SB force field. Regarding this, analyzing two ligands of MolPort-002-530-156 and MolPort-000-410-348 have provided premier binding free energy than 3CL Pro which is known as an inhibitor of Coronavirus protease with sufficient reliability. These two ligands exhibit more acceptable electrostatic and Van der Waals interaction energies, which lead to the creation of stable complexes. The RMSD of the ligands has not been remarkably greater than the RMSD of the protein, reflecting that the ligand has not dispersed away from the primary binding site.
The mean values of RMSD for MolPort-000-410-348 and MolPort-002-530-156 protein complexes have been reported 1.38 and 1.46, respectively. Determining RMSF for assessing regional variations in the protein chain has represented the average value of 0.54 Å and 0.51 Å for MolPort-002-530-156 and MolPort-000-410-348, respectively, which confirms that during the simulation, the residues involved in ligand interactions have kept extremely stable. For confirmation of complex stability, the results of binding-free energy have been evaluated and presented to be À63.34 ± 2.03 and À61.52 ± 2.24 kcal.mol À1 for MolPort-000-410-348 and MolPort-002-530-156, respectively which suggest that the complexes are quite stable. In total, the MolPort-000-410-348 and MolPort002-530-156 components were identified as virtual agents, expecting to show the predicted finding in experimental trials.
Virtual screening technique prediction results to find possible candidates have indicated a list of possible inhibitors against M pro protease and the five top-scored identified hits were assessed by MD simulation [38]. The best hits including AG-690/11060013, AG-690/11203374-2, AG-690/11203374-1, AH-034/04857012 and AG-690/11203374-3 were chosen based on glide score, MMGBSA energy, and ADMET characteristics for computational experiments. Applying the OPLS 2005 force field, the protein structure was refined to constrained molecular mechanics. The Desmond molecular dynamics program was conducted to MD simulation of the detected candidates. The force field for the system was OPLS3e. The analysis data of RMSD for all labeled ligands has suggested a stable ligand-protein complex during the simulation in which the values were estimated ranging from 0.87 to 2.7 for both the protein and the ligand. The 50 ns MD trajectory results have revealed that attaching the mentioned ligands to the receptor leads to potent stability with minor conformational changes. Excluding the terminal amino acids, which have demonstrated mild movement in the 50 ns run, the RMSF values of the whole protein have changed from 0.4 to 2.0 Å, 0.3 to 1.7 Å, 0.4 to 2.0 Å, and 0.4 to 2.0 Å for AG-690/11203374-1, AG-690/11203374-2, AG-690/11203374-3 and AH-034/04857012, respectively. From the outcomes, AG-690/11203374-1 and AG-690/11203374-2 were suggested as the best productive hit compounds, with the greatest interaction number and acceptable binding energy. As a general finding, these two components were proposed among the five high-rated hits as the impressive leads for developing novel antiviral drugs to tackle the CoVs pandemic.
In another virtual screening review, leading inhibitors of novel Coronavirus were sieved, and following that MD simulation has been applied for comparing seven attractive candidates to repurpose against main protease (M pro ) of SARS-CoV-2, involving sapanisertib, napabucasin, ornidazole, daniquidone, lenalidomide, indoximod, and salicylamide could be vital for the treatment of COVID-19 [39]. The MM/GBSA and docking scores of these poten-tial candidates were used to choose them. Control 1 and control 2 have been spotted in a combination of Z element and main protease complex (PDB ID: 5R7Y) with N3 agent. The measured values of 1. 13 1.08, 1.19, 1.194, 1.14, 1.11, and 1.17 Å are reported for control 1, control 2, daniquidone, indoximod, lenalidomide, napabucasin, ornidazole, salicylamide, and sapanisertib complexes, respectively. The mobility and structural stability of the complexes were evaluated using molecular dynamics simulation in the YASARA software package. First, the complexes were purified in the program and the system's hydrogen bond was modified. Considering a cubic simulation box, the AMBER14 force field was assigned. The structural stability of the lenalidomide and napabucasin was maintained for the majority of the simulation period so that for both of them RMSD fluctuation was not considerable. In the case of other complexes, smaller values of RMSD were detected which is the reason for lower flexibility over time. From the standpoint of SASA parameter evaluation, up to 20 ns the complex of salicylamide and main protease has remained stable, while by going to the end of the simulation time a higher enhancement was detected compared to the other complexes. The results of the gyration radius for complexes of salicylamide and sapanisertib have demonstrated higher values, which indicate extra labile structure, as well as the loose packaging mechanism along with the simulation. The overall conclusion has conveyed the stability of all investigated ligand-receptor complexes with the exception that daniquidone, oridinazole, and sapanisertib have demonstrated more optimized efficiency due to higher hydrogen bonds, lower RMSD, and higher binding free energy.
Based on the stream of thought for repurposing licensed medications of other diseases to attach SARS-CoV-2 glycoprotein RBD, 2421 molecules were sieved and docking results of six topranked compounds of cefsulodin, cromoglycate, nafamostat, nilotinib, penfuridol, and radotinib have been evaluated using MD simulation [40]. Spike proteins on the surface of the virus and the human enzyme angiotensin (ACE2) have interceded the binding between the virus and the receptor cell. In this regard, inhibitors with the ability to attach to the Receptor-Binding Domain (RBD) to prevent ACE2 binding would be an efficient way to block the infection. In silico, cefsulodin and nilotinib have been the most successful binders. By building a TIP3P solvated box, the consequence system was provided for modeling with the Amber14SB force field.
The equilibration as well as a 200-ns MD simulation were performed with ACEMD. For a more thorough evaluation of the complexes' stability, the 13 simulations were expanded to 500 ns with strong stability. The computation of binding energy has presented the value of À53.2 ± 4.1 and À41.3 ± 6.7 kcal/mol for nilotinib and cefsulodin, respectively, pointing out that both compounds were known to be the most stable ligands. It is found that accompaniment of cefsulodin with RBD for advancing towards ACE2 has improved the binding manner between the two proteins, which bode well for drug-repurposing techniques aimed at the primary steps of SARS-CoV-2 infection.
Considering observed antiviral properties of several phytocompounds, some promising agents were examined using MD approach as feasible therapeutics against major objectives called RNA dependent RNA polymerase (RdRp) and main protease (M pro ) of SARS-CoV-2 [41]. The MD simulations were accomplished via Patel et al approach with GROMACS ver.2016.4 in which Amber99SB force field was operated. For producing inhibitor parameters and topology, the ACPYPE server was employed. For 50 ns MD simulations, the best-docked component of Mulberroside E, Darunavir were considered as candidate of Mpro inhibi-tors, and Emblicanin A and Remdesivir have been chosen for RdRp. The RMSD results have rendered the mean value of 0.258 ± 0.06 nm, 0.255 ± 0.03 nm, 0.879 ± 0.11 nm, 0.286 ± 0.02 nm and 0. 109 ± 0.01 nm for complex of Mulberroside E with M pro , Mulberroside E, RdRp protein, Remdesivir and Emblicanin, respectively, which obtained values have been in the acceptable range. Also, the significant lower value of RMSD for RdRp with respect to Emblicanin A in comparison to the reference complex has proved that the docked complex remained stable through simulations. Based on the findings of RMSD deviations, Emblicanin A performs fine inside the active site of the RdRp protein. Measuring the binding free energy (MM-PBSA) has led to the values of À111.62 ± 6. 788, À141.443 ± 9.313, À30.782 ± 5.85, À89.424 ± 3.130 kJ/mol for M pro _Darunavir, M pro _Mulberoside E, RdRp_Remdesivir and RdRp_Emblicanin, respectively, that indicates M pro _Mulberroside E and RdRp_Emblicanin A are more preferred compared to other studied inhibitors. In general, the ultimate outcomes confirm that Mulberroside E and Emblicanin A have represented more improved interaction and stability than Remdesivir drug for COVID-19. So, these compounds can be produced as a single or combined treatment for progressing the strategy of multi-pronged therapy against SARS-CoV-2.
In another study, molecular dynamics simulations have been developed to peruse the structural stability of unearthed inhibitors from the screening of Asinex Focused Covalent Library (AFCL) as well as anti-hepatitis-C virus (HCV) remedies of paritaprevir and simeprevir from FDA-approved protease inhibitor library for tackling SARS-CoV-2 3CL pro [42]. Close examination of candidate inhibitors led to choose of three highest rated options of 51, 78, and 223 for detailed 50-ns MD simulation. All of the simulations were conducted by the AMBER 18 simulation package. Also, the Antechamber package in AmberTools was employed to generate ligand parameters using the AMBER force field (GAFF). The RMSD study reveals that involvement of bound inhibitors with main protease has improved the stability in comparison with its apo (without ligand) composition which has showed more oscillation in its backbone. In addition, the smallest averaged RMSD has belong to the complex of M pro /cmp78 with value of 0.5 Å. Based on the binding free energy measurements, all investigated complexes have ascertained that the portion of Van der Waals, electrostatic and non-polar solvation energies have been desirable for interactions of ligand-receptor. Furthermore, about two other ligands, Cmp 78 has had the preferable overall binding free energy with value of À60.05 kcal/mol. In addition, the complexes of anti-hepatitis C virus (HCV) drugs of paritaprevir and simeprevir with SARS-CoV-2 main protease have demonstrated acceptable stability during MD simulation. The mean values of RMSD during the last 40 ns of simulation have been 3.2 and 3.5 Å for paritaprevir and simeprevir, respectively. In addition, evaluation of both inhibitors through binding free energy confirms the appropriate value of À47.15 and À51.84 kcal/mol for paritaprevir and simeprevir. In conclusion, anti-HCV medicines of paritaprevir and simeprevir can help speed up the drug development procedure for medical studies as a COVID-19 therapy.
Kumar et al. [43] have analyzed the performance of two potential molecules of ZINC20601870 and ZINC00793735 come from screening the ZINC database as high scored inhibitors against main protease of the SARS-CoV-2 using MD scheme. MD assessment of target protease with and without selected combinations were carried out through pmemd modules in AMBER18 suite so that the Amber ff14SB force field was applied. The top screened ZINC molecules were created in three dimensions using Marvin sketch, and next evaluated using Gaussian 09 with the B3LYP/6-31G. The binding energies of ZINC20601870, ZINC00793735 in combination with n-Coronavirus main protease have been reported À3.96 and À6.20 kcal/mol, respectively so that the ligand-protein interaction for ZINC20601870 agent has been occurred cos of hydrogen bonds. The stability evaluation with RMSD measurement has illustrated the values of 0.89-2.95 Å and 0.96-3.21 for ZINC20601870 and ZINC00793735, respectively, which affirms that combination of the ZINC20601870 with main target leads to more stable complex rather than ZINC00793735. Besides, from structural stability investigation the RMSF parameter has been detected 1.3-11.80 for ZINC20601870 and 1.38-12.10 for ZINC00793735, which has imparted the interaction of these agents with main protein residues. Also, for providing participants of main protease with ZINC20601870 and ZINC00793735 the free energy has been determined to be À2.66 and À4.55 kcal/mol, respectively. Altogether, it is stated that ZINC00793735 is a possible active compound for combating the novel Coronavirus.
The use of eight phytochemicals sieved from Withania somnifera and Azadirachta indica as well as two repurposed medicines with targeting the spike glycoprotein and the main protease of SARS-CoV-2 have been tested using MD approach [44]. Gromacs 2020.2 software package with AMBER99SB-ILDN force field was adopted for performing the simulations. 100 ns MD simulations have been performed on the effective docked NSP5 complexes containing four phytochemicals (P1-P4) and one repurposed agent (PR5). Withanolide R (-141.96 kJ/mol) and 2,3 Dihydrowithaferin A (-87.60 kJ/mol) have been found with the minimum binding energies for the main protease and spike proteins, respectively. Obtaining the average RMSD for complexes with phytochemicals of P1, P2, P3, P4 as well as the repurposed drug-PR5 and the apo-protein without ligands with values of 0.29, 0.27, 0.25, 0.22, 0.23 and 0.22 nm, respectively, has mean that the ligands do not disperse away from their original binding situations with favorable stable binding rendering them promising inhibitor candidates. Furthermore, the estimated values of average RMSD for complexes with phytochemical ligands in equal order to S1:0.09 nm, S2:0.13 nm, S3: 0.05 nm, S4: 0.05 nm and SR5: 0.17 nm, have con- Investigation of using the combination strategy of two HIV-1 proteinase inhibitors containing lopinavir and ritonavir as favorable candidate against SARS-CoV 3CL pro has been examined using MD simulation [45]. The simulated system of SARS-CoV 3CL pro free enzyme (free SARS), as well as its docking with lopinavir (SARS-LPV) and ritonavir (SARS-RTV), were performed to develop MD trajectories. The AMBER 7 calculation package was exploited to run MD simulations. The SANDER module of AMBER 7 by considering the Cornell force field was utilized to conduct energy minimization and MD calculations. The parm99 force field was employed to determine structure and conformation of combinations. After 600 ps simulation time, the total RMSD values of the three compo-nents have achieved equilibrium. Adding up all of the portions has represented the value of À47.2 and À45.3 kcal/mol for SARS-LPV and SARS-RTV, respectively, which the structural correspondences of both complexes account for their near resemblance. From the results, six and seven hydrogen bonds were discovered for SARS-LPV and SARS-RTV complexes, respectively. Summarizing, it is expressed that the function of LPV and RTV against the 3CL pro enzyme of SARS-CoV have not displayed a remarkable distinction.
The idea of applying HIV protease inhibitors (HPIs) as potential medicine against main SARS-CoV-2 protease (M pro ) was perused in another study by molecular dynamics simulation [46]. For this purpose, simulations were carried out for 10 feasible HPIs including Lopinavir, Amprenavir, Fosamprenavir, Atazanavir, Ritonavir, Darunavir, Indinavir, Saquinavir, Nelfinavir, and Tipranavir with medicinal properties for tackling the target protein. The GROMACS package was leveraged to perform molecular dynamics simula-  tions. The topology was built using the CHARMM36 force field and the tip3p water model. The average values of RMSD for complexes of main protease with HIV inhibitors have determined that the docked complexes with Darunavir and Fosamprenavir candidates create higher RMSD value as well as more oscillations. In addition, the findings of RMSD for the complexes of Amprenavir, Lopinavir, Nelfinavir, Ritonavir, and Tipranavir wereassigned with lower values than that for N3 inhibitor. On the other hand, the radius of gyration analyzing has eventuated the highest oscillation for Atazanavir, while the lowest value of Rg has belonged to Indinavir ligand. A more detailed review of Indinavir and Nelfinavir has revealed the potent RMSD fluctuation for complex with Indinavir, which could indicate that the ligand has undergone a structural change and was put at a disadvantage during the bond formation. In return, the similarity of RMSD value of Nelfinavir with N3 inhibitor was apperceived, confirming that it can maintain a stable structure. Furthermore, Darunavir, Fosamprenavir, and Saquinavir have the highest mean number of hydrogen bonds in comparison to other studied HIV agents. Altogether, Nelfinavir outperformed the N3 inhibitor, especially in terms of the RMSD and Rg dynamical criteria as well as the free energy related to structural biomolecular reactions. Altogether, Nelfinavir outperformed the N3 inhibitor, especially in terms of the RMSD and Rg dynamical criteria as well as the free energy related to structural bio-molecular reactions. Although experimental studies results have shown that Lopinavir and Ritonavir have not been effective in the remedy of SARS-CoV-2 disease, the computational simulation has presented that they are identified as comparable drugs for anti-SARS applications to other studied candidates. To sum up, the examined candidates against targeted protease as well as the best-scored inhibitors in the reviewed studies are summarized in Table 1.
The mean values of RMSD for high score agents corresponding to their SARS-CoV-2 main protease is demonstrated in Fig. 1 [34,35,37,41,44] in which the MolPort-002-530-156 depicted the lowest average RMSD of 1.16 Å. The RMSF values for docked complex of MolPort-002-530-156, Remdesivir and 6LU7-D1-P-D2 are  plotted in Fig. 2 [32,33,37]. As can be seen MolPort-002-530-156 complex represented the least oscilliations in the residues with average value of 0.54 Å. The time-averaged values of the radius of gyration for some selected agents is illustrate in Fig. 3 [32, 41,46]. The compound Mulberroside E presented the lowest Rg in comparison with the other complexes. The Fig. 4 depicts the number of hydrogen bonds for combining some ligand with corresponding main protease, which were time averaged [41,46]. As shown, among them the Tipranavir exhibited the minimum number of hydrogen bonds while RdRp-Emblicanin A indicated the maximum value. Along with various MD trajectory data analyses for ligandreceptor, binding free energy analysis has been carried out to assess the binding capacity between the main protease of SARS-COV-2 and docked complex of identified top molecules using the MM-PBSA and MM-GBSA calculations from the equilibrated MD trajectories. The comparison of binding energies of excellent candidates is shown in Table 2. From the table, it can be perceived that the binding energy values have been in the range of À6.20 (for ZINC00793735) to À105.15 (for Birinapant) kcal/mol. The binding energies have confirmed that overall, all the studied complexes are stable, but the ligand Birinapant have shown a better affinity toward inhibition of main protease in comparison with other complexes [47][48][49]. In addition, after that, candidates Nelfinavir, MolPort-000-410-348, MolPort002-530-156 and Cmp 78 have almost the same binding energies with slight differences, which indicates that they have almost the same affinity toward inhibition of main protease.

Conclusion
Essential requirements to restrain the development of novel Coronavirus pandemic has forced researchers to discover functional agents that can impact SARSCoV2 via various inhibitory effects. In this regard, computational methods such as MD simulation have been able to accelerate the speed of researches, which can be the spotlight on the drug discovery path. Although several researches were done to identify potential candidates, there is no comparison or evaluation between the proposed drugs. Due to the lack of a general assessment of the obtained results and the comparison of the discovered virtual hits in the literatures, the present review has supplied a list of available drugs with potential impressions on the main protease of SARSCoV2 that have been assessed using the MD approach. To address these objects, the intended candidates as well as their bound targets were appraised with the parameters of RMSD, RMSF, binding energy, hydrogen bonds, the radius of gyration, and SASA. In conclusion, some examined candidates were identified as lead that is more effective compounds for fighting against COVID19 from the point of view of the molecular dynamics simulation criteria and are suggested for the forthcoming clinical tests. Table 2 MM-GBSA/MM-PBSA binding energy analysis of best candidates.