An understanding of coronavirus and exploring the molecular dynamics simulations to find promising candidates against the Mpro of nCoV to combat the COVID-19: A systematic review

The first infection case of new coronavirus was reported at the end of 2019 and after then, the cases are reported in all nations across the world in a very short period. Further, the regular news of mutations in the virus has made life restricted with appropriate behavior. To date, a new strain (Omicron and its new subvariant Omicron XE) has brought fear amongst us due to a higher trajectory of increase in the number of cases. The researchers thus started giving attention to this viral infection and discovering drug-like candidates to cure the infections. Finding a drug for any viral infection is not an easy task and takes plenty of time. Therefore, computational chemistry/bioinformatics is followed to get promising molecules against viral infection. Molecular dynamics (MD) simulations are being explored to get drug candidates in a short period. The molecules are screened via molecular docking, which provides preliminary information which can be further verified by molecular dynamics (MD) simulations. To understand the change in structure, MD simulations generated several trajectories such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bonding, and radius of gyration for the main protease (Mpro) of the new coronavirus (nCoV) in the presence of small molecules. Additionally, change in free energy for the formation of complex of Mpro of nCoV with the small molecule can be determined by applying molecular mechanics with generalized born and surface area solvation (MM-GBSA). Thus, the promising molecules can be further explored for clinical trials to combat coronavirus disease-19 (COVID-19).

It is surmised that several different species, including bats, civets, pangolins, and camels, are the intermediate hosts for the transmission of CoVs to humans. The discovery of a host that may act as a bridge between animals and people speeds up the process of transmission of the virus. Throughout history, there have been billions of cases of infection and millions of fatalities all over the globe. This infection causes the respiratory tract of humans and has been identified as a new strain of the coronavirus known as SARS-CoV-2 [2]. At the beginning stage of infection, it showed moderate symptoms of a typical cold, mild-fever, including coughing, a loss of the ability to smell, and issues in the throat. The level of oxygen requirement of patients is not associated with mutation of the virus. However, this may be associated with a load of viruses in host cells. The mutation may change the load of the virus by favorable or adaptation mutation for replication. Fever, cough, and weariness are the most frequent symptoms, while additional symptoms include the production of bodily fluids and others. This virus is a spherical or enclosed particle that contains positive sense single-stranded RNA. This RNA is coupled with a nucleoprotein that is embedded inside the virus's outer coating of protein matrix (capsid) [3].
This review aims to collect information on the importance of MD simulations to find promising candidates to inhibit the activity of the main protease (Mpro) of the new coronavirus (nCoV). In silico approach is important to screen the small molecules against the Mpro of nCoV. Molecular dynamics (MD) simulations are explored but keeping the importance of the MD simulations, the review manuscript has only collected information on the use of various trajectories (RMSD, RMSF, H-bond, and others) from MD simulations and has also focussed on the discussion on various thermodynamic parameters like change in free energy (FE or ΔG) for the formation of the complex between Mpro of nCoV and small molecules. Various tools like Assisted Model Building with Energy Refinement (AMBER), GROningen MAchine for Chemical Simulations (GROMACS), Nanoscale Molecular Dynamics (NAMD) are being used to perform the MD simulations and herein, the authors have collected information on the significance of AMBER and GROMACS to understand the formation of the complex. Further, it is important to mention that the authors have collected information for the inhibition of the crystal structure of Mpro of nCoV using different PDB files available on RCSB. In this review, the authors have collected information on naturally occurring chemicals and repurposed pharmaceuticals that have the potential to function as a candidate for inhibiting SARS-CoV-2. Herein, the information/literature on the inhibition of Mpro of nCoV is taken from 2020 to 2022.

Pathophysiology of SARS-CoV-2 or nCoV
Only four SARS-CoV-2 proteins, including the S protein, are required to construct the true structure of the virus. Fig. 2 shows the genome and protein of SARS-CoV-2 [4]. The assembly of new viral particles and the virus's ability to avoid detection by the host immune system are both controlled by a subset of the virus's other proteins. These non-structural proteins or nsp are produced as two enormous polyproteins, which are subsequently split up into 16 more proteins. A particularly promising therapeutic target is an enzyme known as the major protease, which is responsible for 11 of those cleavages. Spikes visible in the virus are glycoproteins and provide a way to make an entry into a host. The spikes have consisted of S1 and S2 sub-units. S1 functions in the binding of the receptor to the cell of the host, while the S2 gets fused. For getting entry of the virus into the cell, the spike protein has to be directed by an enzyme called protease. Then, the RNA genome begins to convert into structural and non-structural proteins and at the same moment, the replication of the viral genome starts. The blister of the new virus is formed with the target of the endoplasmic reticulum and after blending into the plasma membrane new viruses are released. Simultaneously, the virus enters the host cell, its antigen exposes to antigen presentation cells (APC), and gets recognized by cytotoxic T lymphocytes (CTLs). CD4 and CD8 T cells together make the majority of T lymphocytes (CTLs) which promote the immune system. In infected patients, a reduction in the number of CD4 and CD8 T cells is reported and it results in the prevention of T cell proliferation and its function (Fig. 3) [5].

Variant concern of SARS-CoV-2
The classification of a strain is based on important genetic characteristics. SARS-CoV-2 causes COVID-19 has been mutated that is resulting in different variants of CoV. Some variants are termed Variant of Concern (VOC) as they meet a comparative assessment like an increase in the transmission of the virus as well as the decreased impact of the vaccine or therapeutics. Various variants of the coronavirus have been reported in England, Brazil, California, India, and other nations of the world. A change in the genome of a virus genome due to mutation is known as a variant. Some variants are termed VOC as they meet a comparative assessment like an increase in the transmission of the virus as well as the decreased impact of the vaccine or therapeutics. Up to now, it is unclear whether Omicron is more infectious or causes more severe disease, compared to other variants. Omicron has maximum mutations found than in any other variants [6][7][8].

What is the main protease (Mpro) of nCoV?
In the early 2000 s, some of the Mpro responsible for the SARS-CoV outbreak were Mpro/3CLpro, nsp5. Two overlapping polyproteins named, pp1a and pp1ab, are encoded by the replicase gene and mediate replication and transcription of the viral genome. To govern viral gene expression and replication, CoVs mature through a complicated cascade of proteolytic processing activity on polyproteins. Mpro of CoV has three-domain cysteine protease and is responsible for most of the cleavage inside the forerunner polyprotein. The 3CLpro is the most intriguing drug target for Mpro of nCoV. This is because of its critical function for processing the polyproteins that are translated from the viral RNA to generate a functional replicase complex and enable viral spread. SARS-CoV-2 S protein, on the other hand, has a 22-fold tighter binding to hACE2 than SARS-CoV S protein, which could be one reason for the significantly higher infection rate of SARS-CoV-2. The P1 and P2 locations of the substrate peptide are the most crucial for SARS-CoV 3CLpro substrate selectivity (Fig. 4) [9]. The P1 position has a conserved glutamine residue, but the P2 position contains a hydrophobic segment. 3CLpro identifies and cleaves polyprotein (Met, Leu, Val, Phe)-Gln (Ser, Ala, Gly, or Asn) sequences to produce nonstructural proteins (nsps). In place of leucine or isoleucine at position P2 in the 3CLpro of most coronaviruses, the 3CLpro of SARS-CoV may include phenylalanine, valine, or methionine. The sequences of the 3CLpro cleavage sites of SARS-CoV and SARS-CoV-2 are similar in a lot of ways. [10]. Spike protein (S), RNA-dependent RNA-polymerase (RdRp, nsp12), NTPase/helicase (nsp13), and papain-like protease (PLpro) are all potential coronaviral targets (PLpro, part of nsp3). The S protein's interaction with angiotensin-converting enzyme 2 (ACE2) has been recognized as a critical factor in SARS-CoV-2's entry into human cells. Four redundant furin cleavage sites were discovered in the S protein of the SARS-CoV-2. Furin proteases are prevalent throughout the respiratory system, leading researchers to speculate that the SARS-CoV-2 S protein's cleavage upon departing epithelial cells could result in the virus' highly contagious and deadly character. The acquisition of these furin cleavage sites  provides researchers with knowledge on how SARS-CoV-2 was able to successfully transmit from bats to humans in the first place, in addition to explaining why the virus is highly communicable between humans. A mutant SARS-CoV-2 was found recently that lacked the typical furin cleavage site. They discovered that the mutation produced lower performance in their testing. The Mpro of nCoV is considered a promising therapeutic target due to the dissimilarity to human proteases. Besides this, the genomes and structure of the virus are very closely related to other beta CoVs to help us in drug designing which is based on previous discoveries. It comprised of a long RNA strand, acts mRNA. Upon infection, the virus converted to a new replicated virus by the synthesis of two long poly proteins. These proteins are responsible for making of new RNA and virions by the replication/transcription of complex. Thus, researchers are positively using it to search novel inhibitors that block the action of the Mpro of nCoV in order to develop antiviral drugs through it. Figs. 5 and 6 shows chemical structures of some potent compounds and repurposing drugs/ natural products that showed good activity against the Mpro of nCoV.

Importance of molecular dynamics (MD) simulations over molecular docking
Molecular docking is used to filter the molecules in a library based on the binding energy against the receptor chosen. The binding energy obtained is based on the fitting of small molecules in the cavity of the receptor and depends on the different interactions like electrostatic interaction, Van der Waals interactions, hydrogen bonding, and others. Even, hydrogen bonding can be conventional and non-conventional [11]. Conventional hydrogen bonds (H-bonds) are generated in a dependable manner between molecules that include appropriate functional groups that can operate as H-bond donors and acceptors (A). In a broad sense, the number of hydrogen atoms (H) that may be donated and the number of acceptor sites that are accessible both plays a role in determining the distinct sets of H/ A pairings that are feasible. The molecular docking screening is based on the static phase of the receptor therefore we need to understand the MD trajectory over the binding of receptor and ligands interaction and change in energy throughout dynamics. The pdb file for the Mpro of nCoV has been downloaded from the RCSB (PDB ID: 6LU7 and 6Y2E) and used for the MD simulations [12].

A brief on molecular dynamics (MD) simulations
Molecular dynamics (MD) simulations are performed in the machine to investigate the changes in the atoms of the biomolecules in presence of small molecules. It generates several trajectories of proteins in the presence and absence of the molecules under different force fields and other conditions. A large variety of important processes can be observed with these stimulations such as conformational changes, ligand binding, and protein folding, and can predict the response of biomolecules at the level of atoms such as The C-terminal arrangement (the P1, P2, P3, and P4 residues) and the pro sequence (the P1′, P2′, P3′, and P4′ residues) are bound to the enzyme-substrate (D) Sequences which are sliced by SARS-CoV 3CLpro in the SARS-CoV polyproteins. (E) The consensus sequence for the cleavage through 3CLpro, was examined by the sequence logo platform [9]. mutation, phosphorylation, protonation, or ligand removal. Calculation of forces in an MD simulation may be performed using a model molecular force field. Molecular biology and drug designing can work hand in hand using MD simulations to get the benefit [13].

Basic concept of MD simulations
Molecular dynamics (MD) simulation is established in the 1970 s and it has gone from modeling of hundreds of atoms to biologically molecules like protein/ DNA with and without solvent [14,15]. Generally, MD simulations of systems of 50,000-100,000 atoms are common but they can be performed in a system having few lacs of atoms. This improvement in MD simulations is possible due to the high-performance computer (HPC) and the simplicity of the fundamental simulation method. A basic system model is created by using experimental or theoretical data. The simulation results of different biological systems have multiple degrees. The atomic model is the one that produces the most accurate replication of the real systems. Although, when huge systems or lengthy simulations are necessary, the use of coarse-grained models is becoming more common. The depiction of the solvent is an important aspect of the system definition.
Many strategies have been tried or explored and even though, it results in an increase in the overall size of the simulated systems but the explicit modeling of solvent molecules was found to be the most successful and simplest one [16]. An explicit solvent is capable of recovering the majority of the solvation effects of solvent, even those that originate from an entropic source like hydrophobicity. Once the system has been constructed, the equations that describe the force fields are derived to acquire the forces that operate on every atom, and herein, the potential energy is derived from the structure of the molecule [17].
However, despite their seeming complexity of equations, force fields are not difficult to compute. Because of the simplicity of the force-field illustration of molecular features, including springs for bond length as well as angles, energy and force calculations can be performed very quickly even for very large systems [18]. This is because of the simplicity of the force-field representation of molecular features. The atomistic molecular simulations that are now in use employ force-fields that are parameterized in a variety of different ways [19]. All the force-fields cannot represent all molecules because the recent force-field simulations are similar. Once atom forces are known, Newton's rule of motion is utilized to compute accelerations, velocities, and atom locations. To reduce instability during the integrating movement numerically, one should use a time step less than the molecule's fastest movement. This is generally between 1 to 2 femtoseconds simulations at atomic level and is the most important simulation's bottleneck. Microsecondlong simulations need iterating over this computation cycle 109 iterations. Coarse-grained techniques have this advantage. As the system is simplified, large time steps are achievable to extend the simulations. This may reduce the simulation ensemble's accuracy. Fine-tuning energy calculations, parallelization, or using GPUs have increased MD simulations' performance [20]. Molecular dynamics basic algorithms are shown in Fig. 7 [20].

Monte Carlo simulations
Monte Carlo simulations aren't limited by Newton's equations of motion, unlike molecular dynamics. This flexibility permits smart move approaches inside a statistical mechanics ensemble [21]. These non-trivial changes may speed up equilibrium sampling by 10 10 times or more. Specific Monte Carlo steps may be coupled in a simulation, enabling the modeler to approach a problem with flexibility. Monte Carlo algorithms are readily parallelizable, and some are excellent for big CPU clusters [22].

Theoretical basis for Monte Carlo
A Monte Carlo simulation method is the exact estimation of equilibrium for thermodynamic as well as physical parameters of a system. Consider the possibility of estimating the average value of a property denoted by 'A' Eq. (1). This might be computed by taking into account the following factors: , U = potential energy and r N = the positions of all N particles.
Probability density of finding the system in configuration r N is given by Eq. (2).
where the configurational integral serves as the Eq. (1) denominator. If NMC points can be generated at random in configuration space using Eq. (2), then Eq. (3) may be written as A set of Monte Carlo steps that produce a Markov chain of states make up a Monte Carlo algorithm. Because the chance of generating new configurations relies solely on the present configuration and not on any prior configurations, this Markov process has no history dependency. When our system is in state m, the chance that it will transition to state n is given by the expression mn , where is the transition matrix [23]. Let's establish a probability vector ρ, that indicates the likelihood that a certain system state exists [24]. If someone is simulating in the canonical ensemble, then ρ * is reached when ρ i is equal to the Boltzmann factor (Eq. (4)) for all states. r=r*p (4) The net flow between two states must be zero at equilibrium to satisfy the more stringent requirement of detailed balance (Eq. The matrix representing the transition between states m and n may be expressed as Eq. (6).
= p mn mn mn (6) The Metropolis acceptance criteria is as follows (Eq. (7)) [25]: Monte Carlo simulations may be run in a variety of statistical mechanics ensembles, and the distribution from which we sample depends on the ensemble. The canonical ensemble, in which the number of particles N, volume V, and temperature T are all held constant, the isobaric-isothermal ensemble, in which N, pressure P, and temperature T are all held constant, the grand canonical ensemble, in which, V, and T are all constant, and the semi grand canonical ensemble are all examples of such ensembles [26,27]. For protein simulations, the modeler is unlikely to deviate from the typical ensemble where the partition function is (Eq. (8)):

Ab initio molecular dynamics
It is necessary to find a solution to the electronic properties of various molecules, one should solve the Schrodinger equation. This results in the getting energies and forces, which are referred to as potentials, which make it possible to characterize the dynamics of a reaction following its explicit electronic structure. This strategy is called "ab initio molecular dynamics" (AIMD) [28,29].

Entropy from unconstrained MD
MD accurately models actual systems to illustrate the role of temperatures with finite temperature effects. AIMD has become a standard method for calculating the change in free-energy as well as their deconvoluted enthalpic and entropic contribution in heterogeneously and homogeneously catalyzed systems [30]. FE change can be speedily determined from the distribution of states derived from the AIMD trajectory by using the ΔG = -k B T ln(p/p o ), where p and p o are the probabilities of the state of concern and the reference, respectively. This allows for an accurate estimation of the changes in FE [31,32]. Employing the quasi-harmonic approximation (QHA) with AIMD data is one of the strategies that has garnered the most positive feedback for its effectiveness in compensating for anharmonicity [33,34]. Optionally, improved sampling methods with constrained MD simulations (meta-dynamics approach) can be used to determine entropies with the relationships, ΔS = (ΔU-ΔG)/T, using the FE (G) and internal energy (U) obtained from the simulation. These methodologies can be applied to calculate entropies with relationships [35]. It is important to note that unconstrained AIMD only empowers us to determine ΔS, whereas amplified sampling methods additionally enable us to evaluate ΔS of transition states and it should be taken into consideration [36,37]. By employing the QHA method, the vibrational density of states (VDOS) is calculated by taking the Fourier transform of the velocity autocorrelation function and applying it to the extracted velocities derived from the equilibrated AIMD data [38]. This is done in the following way (Eq. (9)): where v is the velocity and the angular brackets are the average with time. Using this equation, we can accurately define the degree to which systems, at a particular temperature, are anharmonic [39,40].
The VDOS that is produced as a consequence of solving Eq. (9) may then be used to appropriately weight the harmonic partition function v (Eq. (10)).
where N is the total number of atoms and 3 N-n is the total number of degrees of freedom available in the system. When dealing with molecules in the gas phase, both the translational and rotational modes are projected out and then processed independently [41]. After that, the following formula may be used to compute the translational and rotational entropy of AIMD trajectories (Eqs. (11,12)). Determining enthalpy and entropy computationally is critical for comprehending complicated catalytic systems' FE landscape. Entropic factors to the FE change cannot be easily derived from static DFT calculations, unlike enthalpic contributions. Harmonic approximation (HA) is used to examine firmly bonded adsorbates at relatively low temperatures when surface-adsorbate bonding is barely affected. HA is the 2 nd gradient of the equilibrium Born-Oppenheimer energy surface. The ensuing vibrational modes have high frequencies specified by HA. Anharmonicity and cumulative dynamics in catalytic systems make the typical HA technique inaccurate for determining entropy [42]. Atoms and molecules access anharmonic areas of the surface of potential energy due to atomic oscillations. Many weakly bound, solvated, restricted, or high-temperature adsorbates exhibit such phenomena. Determining entropy through anharmonic consideration is critical for thermodynamic and kinetic views of complicated catalytic systems at limiting temperatures [30].
Several computer programs are predominantly used for MD simulations, viz. AMBER, Abalone, CHARMM, CHEMKIN, Discovery Studio, GROMACS, and GROMOS [43]. AMBER software suite has several programs to apply the force field of AMBER to perform the MD simulations of proteins in the presence or absence of small molecules or drug-like candidates. Hamiltonian has been integrated to study or investigate the forces as well as the velocities of the protease with and without drug-like candidates.

Computational Tools (GROMACS, AMBER, NAMD) for MD simulations
GROMACS is open-source software and can be used while AMBER is commercial software. GROMACS can be used easily as the command line and most of the analysis codes are inbuilt. GROMACS is one of the most extensively employed and well-liked open-source bioinformatics programs. It is often used to simulate macromolecules in presence and absence of small molecules. GROMACS is used to do MD simulations of non-biological systems such as polymers as well as biological macromolecules including proteins, nucleic acids, and lipids [44].
Another set of preparation and analysis techniques included in GROMACS is free-energy computation. In comparison, AMBER offers a variety of built-in and tutorial alternatives, such as guided MD simulations and QM/MM simulations [45]. GROMACS requires the patch external codes to function. NAMD is a freeware for MD simulation program and it was developed using the Charm+ + parallel programming language [46]. It is well-known for its efficient and is often used to simulate huge systems. Since then, NAMD has expanded, gaining many new systems and being able to scale to thousands of processors. The most current stable version is 2.9. NAMD is accessible at no cost for non-commercial usage by individuals, academic institutions, and businesses for internal operations [47]. AMBER tool is a model construction and parameter preparation tool, while GROMACS lacks these capabilities for nonstandard residue [43,48].

AMBER
Molecular dynamics (MD) simulations is an effective method that can be used to generate atomic trajectories for a system that contains N particles. This is accomplished by employing integration to Newton's equation of motion, as shown in Eq. (13), and doing so in a way that satisfies macroscopic constraints, such as the evaluation of the passage of time for a set of N interacting particles by way of the solution of Newton's equations.
F=ma (13) Where F represents the force that is acting on an atom, m represents the mass of the particular atom being considered, and a represents the position vector. The use of a molecular force field helps to explain the molecular interactions that occur during the execution of an MD simulation in the program AMBER. The chosen molecule will always have the same response to the force field that was applied with settings. A command line interface (CLI) on a computer running Linux serves as the foundation for MD simulationw in the AMBER program [49].

Force field
AMBER software suite offers a collection of tools that may be used to implement the AMBER force fields into MD simulations of protein. In simulations of molecular mechanics, the force field is an energy function known as a Hamiltonian, and the parameter settings are used to calculate the potential energy of the system. As shown in Eq. (14), the Hamiltonian is integrated with MD simulations to get an understanding of the system's underlying forces and velocities [50].  (14) It is necessary to choose a force field before beginning an MD simulations to comprehend and compute the potential energy of systems. AMBER's FF14SB force field was used for the protein that served as the target [51]. After the force field has been applied, a sophisticated mechanism is used to solvate the molecules with specified water molecules. The system has been shown to have a periodic boundary condition by this simulation. It is essential to emphasize that the periodic box needs to be of sufficient size so that the water may round the particular protein of interest. Therefore, there is no interaction between the periodic pictures and the molecules that are present in the system. To perform MD simulations, TIP3P 8.0 water model is being selected. This particular kind of water box, known as a TIP3P 8.0 box, is designed to solvate within 8.0 Angstroms, which direct that molecules have a buffer of at least 8.0 A° with the wall. The parm7 and rst7 files were then stored in the working directory that was now active. The ff14SB force field was used to determine the values for the parameters [52].

Preparation of input file
It is carried out to minimize the target's energy with and without the drug. After then, the temperature of the system steadily increased. In addition, the MD simulations were carried out on a system with a certain temperature and pressure. However, the manufacture of MD requires three phases, which are as follows: • Heating with constant volume and temperature (NVT) for 20 picoseconds; • Minimization of release strain; • Heating for 20 picoseconds; and • Perform MD simulations for 100, 500, and 2000 nanoseconds with no change in pressure or temperature.
The original system is simplified to reduce the likelihood of steric collisions, and after that, the system is equilibrated. It brings about a state of equilibrium in which the energy, temperature, and pressure average values are maintained. The subsequent stage is a production phase, during which one must establish the qualities of the interest. The trajectory was figured out, and the output file was created and saved after every 500 and 5000 steps, respectively. After that, the parameter and topology file parm7, the coordinate file rst7, and the input files Min.in, Heat.in, and Prod.in were created and utilized to prepare and perform the actual minimizing, heating, and MD simulations [49,52].

Analysis of MD output file
After this, the cpptraj software was employed in AMBER to analyze the data sets and trajectories. Using a variety of various trajectories, the data that was obtained from the MD was evaluated so that further information about the structure, flexibility, and stability could be obtained. The RMSD, RMSF, H-bonds, and change in binding FE calculation are used to do this. The RMSD plot is used to get an understanding of the binding and stability of the complex, as shown in Eq. (15), and the RMSF plot is utilized to gain an understanding of the flexibility of the protein, both with and without the drug, as shown in Eq. (16) [53]. The following formulae were used in the calculation of these values.
Where N = number of atoms, x m , y m , z m = Cartesian coordinates of the initial structure, and x i , y i , z i are the Cartesian coordinates of trajectory at frame t.
Where T is the number of trajectory frames and X is the timeaveraged position.

Molecular mechanics with generalized Born and surface area solvation (MM-GBSA)
MM-GBSA was executed by using the script that was included in the AMBER package. The MM-GBSA methodology is the most effective force-field-based method for determining the binding FE. The binding energy is comprised of the sum of the following components: the Coulomb energy, the covalent binding energy, the Vander Waals energy, the lipophilic energy, the Generalized Born electrostatic solvation energy, the total energy, the H-bonding energy, the Pi packing energy, and the self-contact correction energy [54].
Changes in free energies are employed to discover probable compounds and predict outstanding inhibitors based on the strength of their binding to the target. This is done so to determine which compounds have the potential to be excellent inhibitors. The redistribution of potential energy models, solvent models, and protein flexibility was helped along by docking energy. On the other hand, there are several approaches to use when estimating absolute or relative free energies, the length of the simulation, accuracy, and so on [54].
Using the MM-GBSA approach, the relative binding energy of a complex system with and without drug molecules was determined. For calculating the change in FE, the MM-GBSA approach is applied and it is an efficient, acceptable, as well as in usage [55]. This technique incorporates Generalized Born (GB) electrostatics, molecular mechanics (MM), and solvent accessibility (SA) models (Eqs. (17)-(21)) [56]. Using the provided equations, the relative binding energy terms of complex, target, and drug are estimated using MD simulations trajectories for a different time.

GROMACS
GROMACS is an open-source software program that is employed most often in the field of chemistry. Its basic function is to do dynamical simulations of biomolecules [58]. It offers a comprehensive collection of calculation types, as well as tools for preparation and analysis. There are a few more advanced methods for calculating FE that are supported. It was first established in the Department of Biophysical Chemistry of University of Groningen, and it is currently maintained by contributors in universities and research institutions all over the globe [44].

Force field
Chemistry at Harvard Macromolecular Mechanics (CHARMM) is a prominent library of force fields for MD, as well as the software package used for the simulation and analysis of MD [59]. Various force fields for proteins are CHARMM19, CHARMM22, CHARMM27, CHARMM36, CHARMM36m, and CHARMM36IDPSFF [60][61][62][63]. From quantum chemistry calculations of model compound-water interactions, the atomic partial charges in the CHARMM22 protein force field were derived. Additionally, CHARMM22 is parameterized for the explicit water model TIP3P. Nevertheless, it is often employed with implicit solvents (Eq. (22)). In 2006, a modified version of CHARMM22/CMAP was developed for usage with implicit solvent GBSW [64].
6.6.2. Preparation of input file 6.6.2.1. Generation of the GROMACS input files. To begin, one get the pdb file generally from the RCSB Protein Data Bank. This peptidespecific file format provides the aminoacidic sequence of esculent as well as extra structural information. This format, however, varies from GROMACS structure inputs. To convert the file to GROMACSreadable format, we will utilize the pdb2gmx GROMACS program [65].
6.6.2.2. Set up the simulation box. The next step is to construct the simulation box, after which we will introduce water molecules into the system. The GROMACS application.editconf is used.

Configure the simulation.
To this point, researchers performed a successful configuration of the molecular structures inside our system, as well as the simulation box and the interaction parameters that will ultimately define the dynamics of the system. Configuring the simulation itself, including its duration, step time, temperature/pressure coupling, and so on, is the last step that must be completed before conducting the simulation. This is accomplished via the use of yet another input file known as the.mdp file.

Radius of gyration
The equation for the radius of gyration, also known as Rg, may be found Eq. (23). This term refers to the radial distance of a point from the rotating axis where it is believed that the whole mass of the body is concentrated. In the fields of molecular biology and bioinformatics, the folding and unfolding of proteins is an important process that contributes to the performance of a specific biological function. When it comes to the folding and unfolding of proteins, the topology of the atoms plays a crucial role. However, for proteins with the same chain length, topology on its own is not sufficient to describe the folding. In situations like this, the Rg performs an important and distinctive function in defining the compactness of the system.
Where r i the position of atom i for the center and m i is the mass of atom i.

Root mean square deviation
The ligand is simply docked into the active binding cavity when the molecular docking procedure is carried out. As a result of the ligand being fitted to the protein, the coordinates of some of the atoms in the protein's backbone were moved slightly from their usual location. A comparison of the shift in coordinates of the atoms that make up the protein may be carried out by measuring the root mean square deviation (RMSD). Estimating the atomic coordinates is done via the use of computational modeling so that structural similarities may be compared. RMSD may be used to compare two structures using a single number. The RMSD, is a form of average distance that is measured between the atoms of stacked proteins (Eq. (24)). The equation that is used to compute RMSD in GROMACS.
Where M = = m i N i 1 and r i (t) = position of the atom I at time t.

Molecular mechanics Poisson-Boltzmann surface area analysis
The interaction of low molecular weight compounds with biomolecular systems is mostly investigated in terms of binding FE change. Many computational methods are employed to estimate the FE change of binding. Linear interaction energy, FE perturbation, thermodynamic integration, molecular mechanics, molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) [66], and Generalized Born surface area are the most prevalent methodologies (MM-GBSA). The majority of computational methods alter molecules along an irreversible route and are costly to use. These methods are likewise confined to tiny compounds with a small number of atoms. Linear interaction energy and MM-PBSA are route-independent methods. The FE of binding is obtained using MM-PBSA by evaluating ensembles of the starting and final states. Due to this, the efficiency of MM-PBSA is better [67].
The post-processing computation is performed using the g mmpbsa program in conjunction with GROMACS [68]. This postprocessing following MD simulations is highly effective for determining a more precise binding energy estimate. Because the method is based on quantum physics, the outcome is predictable. Electrostatic energy, polar solvation energy, Van der Waals energy, Weeks-Chandler-Andersen (WCA) energy, solvent accessible surface area (SASA) energy, solvent accessible volume (SAV) energy, and MM-PBSA binding energy were all estimated (Eqs. (25)- (27)) [69].

NAMD
To efficiently simulate complex biomolecular systems, the parallel MD algorithm was created known as NAMD. NAMD can run on low-end commodity clusters with tens of processors, high-end parallel systems with hundreds of processors, and individual desktop and laptop computers [70]. In terms of available features, parameters, and file formats, NAMD is compatible with both AMBER and CHARMM [71].

Force field
When doing an MD simulation on all of the system's atoms, it is standard practice to assume that each atom is subject to a force that is defined by a model force field and takes into consideration how that atom interacts with the other parts of the system. At the moment, these types of force fields give an adequate balance between the accuracy and the computing efficiency of their predictions. The NAMD algorithm makes use of a common potential energy function, which may be broken down into the following components. The first three terms explain the stretching, bending, and torsional bonded interactions, where bonds are the number of covalent bonds in the system, angles is the angle between two covalent bonds that share only a single atom at a vertex, and dihedral is the description of atom pairs that are distinguished by three covalent bonds, having the central bond being subject to the torsion angle [72]. The last two terms in the Eq. (28) explains the nature of the interactions between unbound atom pairs [71]. U total = U bond + U angle + U dihedral + U vdW +U coulomb (28) The parameters ki, bond, r0i, and so on for the interactions provided in Eqs. (28)-(31) are put down in force field parameter files for each particle in a specific context of bonds. These files are used to simulate the force field. Next, the force field is tested to see how effectively it reproduces bulk characteristics and how accurately it reproduces the structure, kinetics, and thermodynamics of small molecules with strong experimental definitions. Parameters are often determined using all the combinations of quantum mechanical calculations and empirical methods. The parameterizations from the CHARMM7 and AMBER [73] force field standards may be used with NAMD without any problems [74].

Free energy calculation
MD may also be employed to build an ensemble of configurations by which thermodynamic variables such as FE differences, and F, can be determined. In summary, there are three alternative methods for calculating ΔF: (1) Calculate the appropriate probability distribution, ρ[U(r 1 ,., r N )], from which a FE difference may be derived using the formula − 1/β ln [U(r 1 ,., r N )]/ρ ο , where ρ ο represents a normalization term: (2) calculate the FE differential directly; and (3) calculate the FE derivative, dΔF(ε)/dε, along some order parameter (collective coordinate), ε compatible with a mean force, and integrate it to yield ΔF [75]. The common umbrella sampling approach, which seeks the probability of locating the system along a particular response coordinate, belongs to the first group [76]. The necessity to anticipate the external potential or bias required to overcome the obstacles of the FE landscape is a glaring flaw of this technique; this issue may quickly become complex in the event of qualitatively novel situations [77].
Based on the literature MD simulations of the natural compound with nCoV were reviewed and reported in Table 1 which are proficient in shackling nCoV [78]. MD simulations of ZINC02123811 by GROMACS against Mpro (PDB ID: 6lu7) indicate equilibration in RMSD over 100 ns, confirming the stability of the best poses that came after docking. Stable and reduced RMSF with ZINC02123811 binding show protein-ligand stability. Rg plot showed Mpro folded stably with ZINC02123811 and acted like Mpro of nCoV free. H-bond plot shows no notable change in intramolecular H-bonds between Mpro of nCoV with and without ZINC02123811 complex. This research found Mpro of nCoV with ZINC02123811 complex is a steady overall simulation [79]. RMSD plots predicted by MD simulations of eribulin mesylate (eri) and Soblidotin (sob) against Mpro of nCoV (6Y2F) are below 0.3 nm, averaging 0.17 and 0.18 nm for eri and sob, respectively. The RMSD for protein-eri and protein-sob complexes, respectively, was 0.23 and 0.24 nm. Although the RMSD profiles of the complexes increased somewhat, they remained within an acceptable range. The Rg values of the protein varied between 2.16 and 2.28 nm. Eri had high Rg values, whereas sob peaked at 93,920 ps, plateaued at 94,910 ps, and then levelled out. The average number of H-bonds with eri and sob was 1.9 and 2.1. This indicates that the ligand is in the required capsule and interacts with the target protein throughout the simulation [80]. MD simulations of Lucenin-2, Saponarin and Isoschaftoside with Mpro of nCoV (6LU7) was done using AMBER. RMSD plots of isoorientin and saponarin imply to be steady within 21-25 ns [81]. RMSD and RMSF plots of Lucenin-2 and Isoshaphoside with 6LU7 are shows in Figs. 8 and 9.
The top MD-generated complex configurations were put into MD simulation by the Desmond software. The kazinol T-protease complexes were steady in a dynamic surrounding at 300 K and 1.01325 bar, according to a simulation of 100 ns. The RMSD value ranged between 1.7 Å. The relative binding energies of the ligand with protein ranged from 30 kcal/mol to 75 kcal/mol from 0 to 100 ns. The small compounds discovered to suppress nCoV protease might potentially be used to inhibit other viral proteases [82]. AMBER 18 was used to conduct MD simulations to investigate the potential of bioactive substances found amla, bhumi amla, and giloy to suppress the enzymatic action of nCoV Mpro. The top 10 were selected employing computational filtering methods. According to simulation results, amritoside, pectolinarin, astragalin, apigenin-6-C-glucosyl7-O-glucoside, 7-ketositosterol, and quercetin engage the binding of substrate cleft of Mpro of nCoV effectively [83]. Ellagic acid (Ela) and kievitone (Kie) interactions with Mpro of SARS-CoV-2 (6LU7) were studied by MD simulations via GROMACS. Mpro of SARS-CoV-2, Mpro of SARS-CoV-2 -Ela, and Mpro of SARS-CoV-2-Kie had average RMSD values of 0.2981 ± 0.0552 nm, 0.2171 ± 0.0267 nm, and 0.3230 ± 0.1366 respectively. RMSD is in the allowed value within 0-10 ns in a situation of Kievitone dynamics with Mpro of nCoV receptor. However, afterward 10 ns, the complex's stability was disturbed. We infer that the catalytic dyad residues stay stable and there were no adequate differences in the RMSF plots of these complexes when correlated to apoprotein. Furthermore, the RMSF values of Mpro of nCoV-Ela (0.1147 ± 0.0501 nm) are smaller as compared to that of the apoprotein, Mpro of nCoV complex with Kie (0.10301 ± 0.1012 nm). For the compete 20 ns simulation, Mpro of nCoV with Ela had averaged intermolecular H-bond creation of 1.81, that is greater than those of Mpro of nCoV with Kie (1.04) and the basic (1.501) [84]. From Fig. 10 we observe the RMSD reading for Mpro of nCoV, Mpro of nCoV -Ela and Mpro of nCoV-Kie respectively 0.2981 ± 0.0552 nm, 0.2171, and 0.2692 ± 0.0433. AMBER 18 was used to perform MD simulations of N1, N10-dicoumaroylspermidine (1), (R)− 3′,7-Dihydroxy-2′,4′-dimethoxyisoflavan (2), and (E)− 2′-Geranyl-3′,4′,7-trihydroxyflavanone (3) against Mpro of nCoV (6LU7). The range of average values is between 1.61 ± 0.02 Å and 2.92 ± 0.06 Å. The complex1 had the greatest divergence (2.92 ± 0.06 Å), whereas the complex2 exhibited the least (1.62 ± 0.02 Å). In the last 125 ns of runs, complex2 and complex3 are more stable in comparison to complex1 based on the average RMSD deviations. Higher RMSD of complex 1 are mostly attributable to changes in the conformational change of extended loops (residues 185-200) and domain III (residues 200-306) RMSF plot, we discovered that domain III of complex 1 is more whippy. Domain III is primarily engaged in the formation of a homodimer, which contributes to greater fluctuations in dietary compounds, while domains I and II are responsible for ligand binding. All three complexes have a comparable degree of compactness. After the binding of a new ligand, the protein's structure and solvent-accessible surface area (SASA) are often considerably altered [85]. Using GROMACS, MD simulations of Table 1 Promising molecules based on natural products to inhibit the activity of Mpro of nCoV studied by molecular docking and molecular dynamics simulations.
Target Kadsurenine and methysticine against Mpro of SARS-CoV-2 (6LU7) for 30 ns had performed to provide dynamic data. Kadsurenin and methysticin RMSD demonstrated that 5 ns and 10 ns respectively reached simulated equilibrium for ligand-protein interactions. Kadsurenine's RMSF was 0.90, 0.14, and 0.52 nm; Methysticin's was 0.64, 0.12, and 0.38 nm. During MD simulations, certain amino acid residues showed steadiness in the protein-complex state owing to phenolic chemicals' less RMSF. During a 30 ns MD simulation, the average distance between two complexes is still 2.2 nm [86]. To confirm the accuracy of the MD simulation by AMBER 18, the RMSD of the produced plots was evaluated at 70 ns. All simulated systems converged within 10 ns and 30 ns, a sign of well-balanced systems; therefore, the insights drawn from the evaluated trajectories are reliable. Hypoxoside binding is characterized by an increase in the deviation of the c-atoms of nCoV, with the attached conformation exhibiting a mean RMSD of 6.12 and the free conformation exhibiting an average RMSD of 5.80. The mean Rg of the bound Mpro of nCoV, and the average Rg of the unbound structure was 19.91 [87]. After the interaction of MQC, the variation of every residue inside the protein were compared. In the instance of MQC complexes, variations were found to be minimal at many residues. Nonetheless, there is a modest rise in the RMSF of MQC, which may be due to the increased size of the attached blocker (Fig. 11). MD simulations by GROMACS for 95 ns for Mpro of nCoV in combination with terpenoid (T3) indicated binding stability. Protein backbone atoms fluctuated more after 20 ns, with a mean RMSD of 0.46 nm. Internal proteinligand T3 fluctuations were smaller than 0.35 nm. The mean RMSF is 0.3 nm and reflects the protein's steadiness following ligandbinding. The smallest number of H-bonds was 0 while the highest was 5 within 10 ns. MD simulations show that T3's binding to SARC single bondCoV-2 Mpro of nCoV is stable. More hydrogen bonds indicate greater ligand-protein interactions [90]. Using the GROMOS, MD simulations of Mpro of nCoV and Mpro of nCoV-complexed with anisotine are run for 100 ns. The RMSD and RMSF plots (Fig. 12) of        Fig. 13.
To access configuration space for testing molecular geometries, 50 ns simulations of the nCoV PR protein both with and without the introduction of the blocker were performed. As a typical model for the MD simulations research utilizing AMBER 18, the docked complexes of the herbal drugs (hesperidin and sesamin) with the lowest binding energy were chosen. RMSD was used to monitor protein structure stability and the placement of ligands in the binding site cleft compared to their most fundamental structure. In comparison to the uncovered nCoV PR protein, steady oscillation and minimal fluctuations of RMSD were found in every complex model, suggesting that the complex forms were more robust and had fewer structural changes during MD simulations [93]. Rutin and Dieckol's MD simulation investigations against Mpro of SARS-CoV-2 (6Y2E) have further confirmed the stability of the selected ligands over time in the catalytic site of the chosen targets, as well as their tightness, the arrangement of conformations in three-dimensional space, and the lack of major variations in binding poses. It has been shown that the dual-targeting inhibitor dieckol has a significant blocking effect on both RdRp and Mpro of nCoV. Moreover, rutin produced from Ruta graveolens, Tephrosia purpurea, and Eucalyptus sp. has a strong   capacity to block Mpro of nCoV with appropriate steadiness, indicating its candidacy as prospective nCoV inhibitors [94]. Amentoflavone, isorhoifolin, nicotiflorin, and rutin persist in equilibrium after 75 ns of MD simulation with Mpro of nCoV (6LU7). No complex experienced structural instability throughout the simulation, excluding the naringin complex with Mpro of nCoV. From 42 ns, RMSD oscillates by greater than 3.0 Å indicating instability. Maximum RMSD deviations of 3.0 Å suggest equilibrium. A mentoflavone had more stable RMSD curves than isorhoifolin, nicotiflorin, and rutin by roughly 1.0 Å. The amentoflavone Rg trajectories reveal uniform variations of Rg as a function of the molecule trajectory. Isorhoifolin, rutin, nicotiflorin, and naringin have large Interquartile ranges values, implying their Rg values change considerably based on molecular trajectory. 5.31, 4.51, 4.25, and 4.84 Å. Although the interquartile range of Rg is greater for these situations, it never exceeds the amentoflavone values by as 0.6 Å, demonstrating that the ligands remain very stable over the 75 ns simulation for Rg, with only modest conformational changes [95]. The dynamic binding interactions of five inhibitory drugs were studied, and AMBER 14 simulations of ligand-protein complexes were run for 100 ns. The complexes achieved equilibrium within 5 ns of the simulation, according to the calculation of the RMSD of the ligand trajectory, which showed values between 1.5 and 3.5. Significant alterations in vicenin-2 and narcissoside suggest a flexible entrance to the active site of 3CLpro. For example, molecules with higher stability include rutin, isoschaftioside, and kaempferol-3-O-gentiobioside. The degree of configuration change for all of these compounds within the binding pocket is fairly stable with an RMSD of less than 0.5. (Fig. 14). To assess each ligand's ability to bind to 3CLpro, the binding FE was calculated using MM-PBSA. [96].
RMSD for rutin is 0.39 nm. The RMSD for (R)-amygdalin is somewhat greater than for rutin, indicating lesser stability. Rutin and amentoflavone docked better with SARS-CoV-2's protease. The flavonoid amentoflavone shows few changes from the beginning conformation with a mean RMSD of 0.13 nm, indicating its combination with the major protease is stable (Fig. 15). In SARS-CoV-2 complexes, rutin and amentoflavone create 12 and 7 hydrogen bonds, respectively. During simulation, rutin and amentoflavone produced 5.95 and 2.05 hydrogen bonds, respectively. Rutin forms more hydrogen bonds than amentoflavone. Rutin may bind better than the other two chemicals owing to hydrogen bond formation [97].
MD analysis was done between Mpro of nCoV (6LU7) complex and its apo form by using GROMACS software. The mean RMSD of Mpro of nCoV complex and apo form was 0.37 and 0.29 nm, respectively. Mpro of nCoV had lower backbone atom deviations with 5GDE than apo. Average Rg values were 2.14 and 2.15 nm for apo and complex form, indicating comparable simulated compactness. RMSF studies indicated that the apo and complex of Mpro of nCoV, residues had mean RMSFs of 0.21 and 0.15 nm, respectively. In complex form, Mpro of nCoV backbone atoms fluctuated less than in apo form. The MD simulations revealed that the configuration of the Mpro of nCoV complex is compact like its apo configuration, but with improved structural stability and less residue fluctuation [98]. Desmond was used to conducting 100 ns MD simulations of Mpro of nCoV (6LU7) complexed with Amentoflavone and Morusin. During the whole simulation, the ligand-protein combination remained nearly constant, with RMSD values ranging from 0.8 to 3.0 Å for both morusin and the major protease complex. By measuring the RMSF of each major protease residue, the flexibility of the protein was assessed further. The RMSF value of the whole main protease, excluding the terminal amino acid residues, was between 0.5 and 1.5 Å. This may be because the N and C terminals move more quickly than the rest of the protein. The morusin associated well within the primary protease binding site, with Glu166 and Gln192, in particular, demonstrating effective ligand interaction throughout the simulation period [99]. The MD simulations of C00006660, C00014803, and ANLT0001 were done via GROMACS. C00006660 and ANLT0001 exhibited greater RMSD swings throughout the MD run. Throughout 200 ns, RMSF had also been measured to determine the atomic mobility of structurally changing backbone atoms (Figs. [16][17][18]. In general, it was observed that molecule-specific C00014803 docked protein residues exhibited less variation than other top docked hits [100]. MD simulations results of repurposing drugs against Mpro of nCoV were reported in Table 2 which is efficient in inhibiting the action of Mpro of nCoV. GROMACS is employed to simulate the linkage of hydroxychloroquine (HCQ) with the main protease of nCoV (6LU7). RMSD data indicates that the protease moiety undergoes a dramatic structural shift at 2 ns, just as the medication starts to engage more closely with the protease. In around 7 ns, the drugprotease system reaches equilibrium. In addition, the RMSD figure demonstrates that the primary protease backbone becomes destabilized upon contact with the drug molecule. It is intriguing that when the drug-protease complex develops, the number of H-bonds within the drug and protease grows but the average H-bonding distance decreases. The current MD and MD simulations experiments reveal that the medicine HCQ may have a major influence on the structural main protease of nCoV, which may result in a strong inhibitory effect [101]. MD simulations of the inhibitor (118098670)protease system demonstrate that the system diverges little from its x-ray crystal structure, with 3 Å for the protease and 1.5 Å for the blocker from their original structures. The MD simulations between inhibitor (104161460) and Mpro of nCoV reveal that the protease has an RMSD of 3 Å. MD simulations of the blocker protease system  Interestingly, the RMSD of the complex system was more stable than the Mpro of nCoV. In addition, the residues 130-299, which correspond to the interaction region of the Mpro of nCoV-conjugate complex based on MD simulations, had substantially smaller RMSF statistical disparities. Intriguingly, the binding area for lead conjugation was also within residues 137-298 before MD simulations, and this region remained unchanged after simulation, indicating an effective binding of Nos-HCQ with the Mpro of nCoV that gives stability for the Mpro of nCoV-conjugate. The resulting plots was then estimated to determine the Rg of Mpro of nCoV with and without noscapine. The lead combination with the Mpro of nCoV exhibited Rg fluctuations between 2.12 and 2.26 nm and for Mpro of nCoV between 2.14 and 2.26 nm, demonstrating the steady and strong interaction of Nos-HCQ with Mpro of nCoV [103]. Theaflavin digallate, amodiquine, and lopinavir were all simulated using MD simulations against the Mpro of SARS-CoV-2 (6LU7), and the structural variation was identified using the RMSD data of the protein with ligand complexes from 0 to 20 ns. In the range of 0-5 ns, RMSD data gradually climbed until they reached a steady level. Amodiquine, lopinavir, and theaflavin digallate were all found to have mean RMSD data of 0.25, 0.23, and 0.22 nm, respectively. The fluctuation of each atom during the simulation is given via RMSF. The binding site residue changes are observed less when the RMSF for a protease including 306 amino acids and 3 potential drug compounds was estimated. The RMSF readings for amodiquine, lopinavir, and theaflavin digallate, on the mean, were 0.15, 0.17, and 0.2 nm, respectively. First, the Rg of protein-ligand complexes was found to be 2.13-2.20 nm. In the case of lopinavir and theaflavin digallate, input values were reduced and then stabilised between 5 and 20 ns, indicating a stable binding conformation. After 5 ns, the Rg values of the Mpro of nCoV with amaodiquine stabilised [104]. The apo-form of the Mpro showed the greatest RMSD values, but docked complexes had equivalent or lower value systems, indicating how docked ligands stabilized the structure of the protease. Once the locations of these ligands were set, they did not alter much. CID3010243 and CID44271958 had the greatest RMSD scores, indicating the greatest position change, whereas CID10009410 and CID44271905 had the least. As predicted, major protease terminal residues revealed greater RMSF values, demonstrating their mobility. While attached to CID3010243, residues 49, 52, 140,169, and 189 were more mobile than in apo-protease. Similar to residue 120, residues 186-190 were more flexible when bound to CID44271958. Different ligands modify some residues. Rg of the Mpro of nCoV was determined for every complex to determine the influence of the ligands. Rg reading of    [105].
The use of plitidepsin in the fight against COVID-19 may be acceptable throughout the development of antiviral medications. It has been stated that plitidepsin is an even more promising option than remdesivir in the fight against the SARS-CoV-2 infection. In the study of Vishvakarma et al., the ability of plitidepsin to suppress the activity of Mpro of nCoV was evaluated. At 300 K, MD simulations were performed to find the binding through the complex of plitidepsin and the Mpro of SARS-CoV-2 (Fig. 19). Additionally, the binding was examined at a greater temperature (325 K). At 300 K, plitidepsin and the Mpro of nCoV form a stable complex [107]. Authors have performed the MD simulations repurposed drugs acyclovir (A), ganciclovir (G), and a ganciclovir derivative (D) with Mpro of nCoV using WebGro at 290, 300, 310, and 320 K to achieve more reliable findings (Fig. 20). The results of MD simulations showed effective binding and inhibition [2]. In the beginning, it was concluded that HCQ and remdesivir are somewhat active, and now both  GROMACS using GROMOS96 43a1 force field N-substituted isatin derivatives and pyrazolone compounds could be used as a potent inhibitor and may possess an antiviral activity against SARS-CoV-2 [106] of these substances are being used as drugs for the management of COVID-19. Keeping this information in mind, Deepak et. al. have designed various compounds/ derivatives (Fig. 21) and study how they inhibit the Mpro of nCoV. Based on the results of the docking investigations, they have concluded that DK4, DK7, DK10, DK16, DK17, and DK19 have the potential to serve as inhibitors for Mpro of nCoV and therefore, effectively regulate COVID-19. In addition, MD simulations were done on the molecules HCQ, DK7, THC, and DK16 since these were the molecules with the best score results obtained from the docking investigations. To confirm the stability of the protein-ligand complexes, RMSD, RMSF, and H-bonds have been investigated. The results indicated that DK16 displayed better stability in comparison to both the parent molecule and DK7. The stability of the DK7 molecule with the Mpro protein of SARS-CoV-2 was shown to be equivalent to that of HCQ. It is hoped that both of these compounds will be viable candidates for medication treatment and that they will be able to undergo further experimental validation [108]. Patients who suffer from persistent congestion and sinus problems may find relief from taking amino acids. In addition, it may lessen the severity of infection by focusing on the reproduction of the virus; it could also be useful in cleaving the virus, etc. The efficacy of the amino acids to bind to the Mpro of SARS-CoV-2 was evaluated. With binding energy of − 0.94 kcal/mol, Arginine had the greatest affinity for the enzyme. The thermodynamic parameters change in enthalpy (ΔH) and relative change in entropy (ΔS) were derived from MD simulations (Fig. 22) conducted on the SARS-CoV-2 protease with and without arginine. This made it possible to calculate the change in FE (ΔG) for the building of the complex, which was determined to be 2.74 kcal/mol [109].

Challenges
Molecular dynamics (MD) simulations of biological systems have a great difficulty overcoming perception problems among the wider scientific community as the results obtained from computational tools may differ from the wet laboratory experiments. Although, the importance of MD simulations is still constrained by two primary obstacles: the force fields that are used require further refinement, and high computational demands prohibit routine simulations of longer than a microsecond time scale, which in many instances leads to insufficient sampling of conformational states. It takes a huge  time to complete a one-microsecond simulation of a very modest system.

Future prospectives
MD simulations provide a strategy to screen millions of molecules in a short period in comparison to wet laboratory experiments. Using MD simulations and MM-GBSA calculations, one can have an idea of the promising nature of the molecule against the Mpro of nCoV. This is a time-saving approach to finding promising drug candidates to combat COVID-19. One important thing has to be done to corroborate the wet laboratory experiments by extending the MD simulations for several microseconds at different temperatures. Various research groups performed MD simulations in a few microseconds but still need to further extend it. Further, MM-GBSA or MM-PBSA calculations can provide useful information like change in FE from the change in enthalpy and change in entropy to know the feasibility of the formation of the complex at a particular temperature [110].

Conclusion
The efficient implementation of MD simulations for fast drug discovery would result in considerable time and money savings. New methods, software, and hardware have made it easier for the pharmaceutical industry to use MD simulations approach in a big way. GROMACS, AMBER, and NAMD are the tools used for MD simulations. The authors focussed on Mpro of nCoV as a target because it plays a vital role in facilitating viral replication and transcription. The exponential increase in the cases due to coronavirus has attracted researchers to find effective molecules to inhibit the activity of this virus. But, finding a molecule or drug is not an easy task, as it is very difficult to claim the biological potency of a molecule without clinical trials which need a longer timeframe. Therefore, the researchers have started exploring the computational tools (MD simulations, molecular mechanics with MM-GBSA or MM-PBSA calculations) to get a prominent candidate against the Mpro of nCoV. The reason to target the Mpro of nCoV is due to its involvement in the replication of virus. In this context, researchers have used repurposing drugs, natural products, and new candidates to inhibit the activity of Mpro of nCoV. Various trajectories (RMSD, RMSF, Rg, hydrogen bond) are obtained to acknowledge the change in the structure of the main protease of nCoV. MM-GBSA or MM-GBSA calculations are carried out to determine the change in free energy through the change in entropy and change in enthalpy for the complex formation. This approach enables us to filter the compounds for the inhibition of this virus and researchers have suggested it for their clinical trials.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.