Introduction

In December 2019, a novel coronavirus known as COVID-19 was initially detected in Wuhan, China. Since then, it has quickly disseminated to numerous countries globally [1, 2]. Fever, sore throat, exhaustion, diarrhea, loss of taste or smell, dry cough, runny or stuffy nose, and digestive problems are among the symptoms of the new coronavirus disease. Patients may occasionally get critically unwell and endure acute shortness of breath [3-4]. In March 2020, the World Health Organization (WHO) declared COVID-19 as a worldwide pandemic [5]. Several vaccines for COVID-19 have been introduced, and various scientific studies have been conducted to develop antiviral drugs.

Nonetheless, antiviral drugs are unavailable to target SARS-CoV-2 [6, 7]. The most effective approach to treating SARS-CoV-2 infection may be targeting viral particles, which can inhibit the infection [8]. A key player in coronavirus replication is the 3CLpro (3-chymotrypsin-like protease) enzyme, and analysis of its amino acids has shown similarities between COVID-19 and other coronaviruses [9, 10]. Computer-based drug design can be utilized to identify compounds, such as FDA-approved drugs or natural products, that have already demonstrated safety and efficacy in humans [11, 12]. Natural products, known for their low toxicity and potent biological activity, are particularly valuable in the drugs development for various diseases [13a-b]. The secondary metabolites purified from medicinal plants are considered to be highly potential. Therefore, their work is a model in drug development for new antiviral drugs against Sars-CoV-2[14]. Solanaceae is abundant in alkaloids and contains 90 genera and more than 3000 species on all continents. Datura stramonium and Hyoscyamus niger are very important medicinal plants and contain alkaloids, flavonoids, tannins, terpenes, saponins, carbohydrates, cardiac glycosides, and anthraquinones compounds. They are rich sources of tropane alkaloids. D. stramonium is widely distributed throughout the world (Europe, Asia, America, and South Africa). Likewise, H. niger is distributed in Europe and Asia [15a-b]. These plants have shown curing (anti-diabetic, antioxidant, anticancer, anti-asthmatic, anti-inflammatory, antiviral, anti-microbial, etc.) and immune-boosting properties. Parts of these plants, like leaves, fruits, stems or seeds are used to prepare drugs [16-17]. Some alkaloids have been used for various diseases such as bronchial asthma, infection of pulmonary regions and tumors, diarrhea, gastrointestinal colic, anti-tumor, mydriatic, anti-spasmodic, anti-cholinergic, anti-muscarinic, neuralgia and manic psychosis, etc. [18a-b]. Alkaloids, an extensive group of secondary metabolites characterized mainly by at least one nitrogen atom, may help seek drug treatment for COVID-19. Many studies report the antiviral effects of many compounds with alkaloid structures [19-24]. Thus, the antiviral nature of alkaloids opens a route for further research of their possible potential against Covid-19 Mpro. The objective of our research was to investigate the antiviral effects of natural substances on the primary protease of SARS-CoV-2. To achieve this, we utilized a computational approach to screen 73 alkaloid molecules as potential candidates. The impact of these selected alkaloid molecules on the structure of the main protease was assessed through molecular docking and molecular dynamics simulations.

Material and Methods

Protein/Macromolecule

To conduct molecular docking studies, we retrieved and prepared the 3-dimensional crystal structure of the 3CLpro, from the Protein Data Bank (PDB) website. We downloaded the pdb file format of the 6LU7 protein, which consists of chain A.

Ligands

We obtained the chemical structures of the ligands from the PubChem database in SDF (Structure-data file) format. We used Biovia Discovery Studio to convert the structures into pdb format. The PubChem IDs of 3D structures of the compounds obtained in this study are Hyoscyamine (CID: 154417), Moupinamide (CID: 5280537), Littorine (CID: 443005), 7-Hydroxyhyoscyamine (CID: 71587892), and 3-tigloyloxy-6-isobutyryloxy-7-hydroxytropane (CID: 129846249). These compounds have been extracted by researchers before and we used them in this study. These have exhibited anti-inflammatory, antiviral, anti-microbial, and immune-boosting properties.

Molecular docking

Molecular docking is utilized to differentiate and assess protein-ligand interactions in order to observe the binding affinity and protein-ligand docking platform [25-26]. In this study, we utilized AutoDock Vina [27] to investigate the binding of 73 alkaloid compounds to the 6LU7 Mpro protein. To optimize the protein, we employed ADT to add polar hydrogen atoms and Kollman partial charges. We prepared the ligands by eliminating water, incorporating Gasteiger charges, and determining the count of torsion atoms. Both the ligand and receptor files were converted to PDBQT format. Subsequently, the grid box for each target was selected based on the x, y, and z coordinates. For each ligand, we generated 20 conformations during each docking run. We chose the top 5 ligands with the lowest binding energy for the flexible docking process from the prepared compounds. All Vina results were prepared by Auto Dock Tools, discovery studio, and LigPlot+ [28-29].

Molecular dynamics (MD) simulations

We used GROMACS 2019.3 software to perform all atom MD simulations and explore the dynamic behaviour of 3CLpro and 3CLpro: ligand complexes [30]. The MD trajectory was visualized using the VMD 2.9 program [31]. The initial coordinates for 3CLpro, ligands, and the complexes of 3CLpro with ligands were acquired based on the outcomes of the molecular docking process. Bonded (bond lengths, bond angles, dihedral, and improper angles) and nonbonded (electrostatics and Van der Waals) potential energy functions of all species were taken from AMBER 99SB-ILDN force field. The Antechamber program was used to generate the necessary parameters for MD simulations of the ligands [32]. Following this, the ACPYPE program was used to convert the Antechamber output files into topology and coordinate files compatible with GROMACS software [33]. For assigning partial charges to proteins and ligands, we adopted the restrained electrostatic potential (RESP) approach [34]. Short-range non-bonded interactions were modelled using a 12-6 Lennard-Jones potential with a cutoff of 12 Å. Periodic boundary conditions were applied in all directions, and the Particle mesh Ewald (PME) method was used to calculate long-range electrostatic interactions. The 3CLpro: ligand complex was placed in the centre of a solvation box with a cubic shape. This solvation box included water molecules modelled using the TIP3P model. The solvation box had a solvation shell of 10 Å. To balance the system, 4 sodium ions (Na⁺) were included in the simulation box as 3CLpro has a charge of -4. Furthermore, to more accurately simulate physiological conditions, a salt concentration (Na⁺Cl^-) of 0.145 M was implemented [35]. The simulation box had a length of 9.64 nm in three dimensions. To minimize thermal noise and potential energies in the solvated dendrimers, we employed the steepest descent minimization algorithm followed by the conjugate gradient method.

The equilibration simulations were conducted in two phases. In the first phase, the canonical (NVT) ensemble was used for 500 ps to bring the simulation system to 310 K. In the second stage, a process called isothermal-isobaric (NPT) equilibration was conducted for 1 ns. This was done to ensure that the pressure within the system remained consistent at 1 bar. After completing the equilibration phases, the position restraints were released, and production MD runs were carried out for 50 ns using the NPT ensemble.

Results and Discussion

Molecular Docking

The main protease 3CLpro of SARS-CoV-2 contains conserved catalytic residues, His41 and Cys145, in its active site, which are crucial for its enzymatic activity [24,25]. To inhibit virus replication and alleviate the burden on the host, it is important to block these proteases. Thus, we focused on targeting the 3CLpro to identify potential compounds that can specifically block its two catalytic residues due to the fact that all residues of the receptor were rigid except for two selected residues CYS 145 and HIS 41.

The complex between moupinamide and the Mpro enzyme showed the most favourable binding energy of -7 kcal/mol. Moupinamide formed hydrogen bonds with several amino acid residues including MET165, ARG188, THR190, GLN192, GLY143, and SER144, as well as CYS145. Moreover, a Pi-alkyl interaction was observed with CYS145 and a Pi-sulfur bond formation occurred with MET165. The ligand also formed a carbon-hydrogen bond with ASN142. Furthermore, Van der Waals attractive forces were observed between moupinamide and GLU166, GLN189, HIS163, LEU141, ARG188, and MET165 residues (Figure 1).

The complex formed between littorine and the Mpro exhibited a favourable binding energy of -6.8 kcal/mol. Interactions are divided into 5 categories namely: hydrogen bond, Alkyl, Pi-sulfur, Pi-donor H-bond and van der Waal bonds. CYS145, SER144, GLY143, HIS163 were seen bonding with hydrogen bond. CYS145 formed a Pi-sulfur bond, and CYS145, HIS41, and LEU27 were involved in the protein-ligand complex through alkyl bond formation. The GLU166 residue was found to interact with littorine through a Pi-donor hydrogen bond. Along with these major interactions, 11 residues were forming Van der Waals attractive interaction: THR26, THR25, MET49, HIS164, MET165, GLN189, PHE140, ASN142, LEU141, LEU27, and HIS41 (Figure 2).

3-tigloyloxy-6-isobutyryloxy-7-hydroxytropane-M^pro complex gave -6.4 kcal/mol as the least binding energy. The compound interacted 3 hydrogen bonds: CYS145, HIS41, and HIS163. Furthermore, HIS172 and HIS163 made π-interactions with alkyl group. Among the

Figure 1: The interaction between moupinamide and Mpro using LigPlot+, Auto Dock Tools, and Discovery Studio, both in 3D and 2D schematic representations

different residues, HIS41, MET49, and MET165 were seen bonding by conventional alkyl bond. GLN189 and ASN142 residues can were bound via carbon hydrogen bond. Van der Waals bond formation were seen in HIS164, TYR54, ARG188, ASP187, GLY143, SER144, PHE140, LEU141, GLU166, and MET49 (Figure 3).

The complex between hyoscyamine and Mpro, which has a binding score of -6.2 Kcal/mol, has been observed to form hydrogen bonds with SER144, CYS145, LEU141, and HIS163 residues. Likewise, HIS41 and LEU27 were interacting using an alkyl interaction. CYS145 alone was forming pi-sulfur bond. The weak Van der Waals interaction can be seen in MET49, THR25, GLN189, MET165, HIS164, PHE140, ASN142, GLY143, and THR26 (Figure 4).

The complex formed between 7-Hydroxyhyoscyamine and Mpro had minimum binding energy of -6.1 kcal/mol, and two different types of bonding were observed during docking. 4 conventional hydrogen bond were observed with THR26, GLY143, SER144, and CYS145 residues. The ligand was found to form a Van der Waals interaction with THR25, GLN189, MET49, HIS41, GLU166, MET165, HIS163, LEU141, GLN189, and ASN142 (Figure 5).

Van der Waals interactions and hydrogen bonds were significant factors in the binding process. Van der Waals interaction, the weakest intermolecular attraction between two molecules, can become very strong with multiple interactions. In protein-ligand interactions, hydrogen bonds helped stabilize the ligand, while hydrophobic and Van der Waals interactions also contributed to stabilizing nonpolar ligands. Analysing the structure of ligands and docked poses is typically used to understand the interactions and binding energy of compounds. As a result, molecular docking programs are widely used in drug discovery to identify potential drug compounds based on less affinity [36-37]. According to molecular docking results, the five alkaloid components, namely Moupinamide, Littorine, 3-tigloyloxy-6-isobutyryloxy-7-hydroxytropane, Hyoscyamine, 7-Hydroxyhyoscyamine, have showed the lowest binding energies of -7, -6.8, -6.4, -6.2, and -6.1 Kcal/mol, respectively, and the interactions indicated that the phyto-compounds under study have numerous hydrogen bonds and Van der Waals interactions, indicating that the ligands are

Figure 2: The interaction between littorine and Mpro using LigPlot+, Auto Dock Tools, and Discovery Studio, both in 3D and 2D schematic representations

Figure 3: The interaction between 3-tigloyloxy-6-isobutyryloxy-7-hydroxytropane and Mpro using LigPlot+, Auto Dock Tools, and Discovery Studio, both in 3D and 2D schematic representations

stabilized within the complex. These alkaloids may prepare a platform for in silico study of their antiviral potential against the main protease of COVID-19. Moupinamide has made greater number of hydrogen bonds with residues of MET165, ARG188, THR190, GLN192, GLY143, SER144, and CYS145. It has also made Pi-alkyl, Pi-sulfur and carbon-hydrogen bonds. Other alkaloids have made 4,3,4, and 4 hydrogen bonds, respectively (Table 1).

Moupinamide delineated better antiviral inhibition against SARS-CoV-2 Mpro, because more number of hydrogen bonds shows the better inhibitory and binding efficacy [38].

Molecular dynamics (MD) simulations

All-atom molecular dynamics (MD) simulations offer a dynamic representation of the behaviour of protein-ligand complexes over time. Various analysis tools, such as root mean square

Figure 4: The interaction between hyoscyamine and Mpro using LigPlot+, Auto Dock Tools, and Discovery Studio, both in 3D and 2D schematic representations

deviation (RMSD), root mean square fluctuation (RMSF), centre of mass (COM) separation distance between the protein and ligands, radius of gyration (Rg), solvent accessible surface area (SASA), radial distribution functions (RDF), protein-ligand interaction energies, and intermolecular hydrogen bonds, are used to investigate and understand the complex interactions.

Equilibration and relaxation of protein-ligand complexes

The root mean square deviation (RMSD) is determined by calculating the average distance between the corresponding atoms of a protein backbone over the course of molecular dynamics (MD) simulation. Hence, it can be used as the first benchmark to explore the stability and equilibrium state of protein in its native and complexed forms. In Figure 6, the RMSD values of the atomic positions of the protein backbone in the Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro complexes are shown. These values are calculated by comparing them to the initial structure. In addition, the RMSD plot of the native Mpro is included for comparison. From the graph, it can be observed that there are no significant changes in the protein configuration during the last 15 ns of the molecular dynamics (MD) trajectories. The well-equilibrated structure of Mpro and the sufficient simulation time indicate that the simulation can provide valuable information about the protein's microstructure, the nature of the interactions between protein and ligand molecules, the distances between protein and ligand molecules, and other relevant factors. The statistical analysis was conducted using data from the last 10 ns of the MD simulations, ensuring that the simulated systems reached an equilibrium state. The average root mean square deviation (RMSD) amounts for the protein Mpro, the complexes of Mpro with hyoscyamine, moupinamide, and hydroxytropane, were approximately 0.30 ± 0.02, 0.30 ± 0.01, 0.15 ± 0.02, and 0.29 ± 0.02 nm, respectively. This indicates that the structure of protein is considerably more stabilized after binding to Moupinamide. In the case of Hydroxytropane-M^pro complex, a peak with a deviation of 1 Å is observed at 35 ns, which may be caused by local changes in the protein microstructure; however, this amount of RMSD changes are known to be acceptable for protein containing simulated systems.

Figure 5: The interaction between 7-Hydroxyhyoscyamine and Mpro using LigPlot+, Auto Dock Tools, and Discovery Studio, both in 3D and 2D schematic representations

Table 1. Interactions and binding energies of top 5 conformations of ligands with Mpro

Figure 6: The RMSD (Root Mean Square Deviation) of the Mpro backbone was calculated for the complexes involving Mpro, Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro

Structural characterization of M^pro-ligand complexes

Figure 7 illustrates the root-mean-square fluctuation (RMSF) values of the protein residues in the trajectory. These values were calculated by aligning the residues to a reference frame that represents the average position of each residue over time. This parameter allows us to assess the flexibility of different regions within the protein structure and identify areas that display significant deviations or minimal deviations from their average conformation.

Though few residue regions of Moupinamide-M^pro can be seen to show higher fluctuations in comparison with M^pro (0-2, 4, and 48-52), the protein residues are generally more fluctuated especially in the case of 110-150, 165-172, 191-205, and 301-306 residues. The residues 115-144, 166-172, and 190-197 of Hyoscyamine-M^pro complex show significantly higher fluctuations compared to their corresponding regions in the protein. Finally, a comparison between RMSF plots of M^pro and Hydroxytropane-M^pro reveals that except the 1-3, 154, and 297-299 residues, which are more fluctuated in Hydroxytropane-M^pro, protein structure shows higher fluctuations especially in 126-141, 166-172, 190-196, 222, 232, 244, 277, and 278 residues. The average root mean square fluctuation (RMSF) values for the free Mpro protein, the Mpro protein bound to hyoscyamine, the Mpro protein bound to moupinamide, and the Mpro protein bound to

Figure 7: During the production molecular dynamics run, the protein's root mean square fluctuation (RMSF) per residue was computed for the complexes comprising Mpro, Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro

Figure 8: These images depict the stable arrangement of the Hyoscyamine-Mpro complex and the amino acids that interact with the Hyoscyamine molecule

hydroxytropane were 0.16, 0.13, 0.13, and 0.14 nanometers, respectively. The RMSF analysis suggests that binding of alkaloids reduces fluctuation of protein building blocks.

In Figures 8-10, we can observe the complexes formed between the protein and ligand molecules, as well as the distribution of amino acid residues located within 5 Å of the ligands after 50 ns of molecular dynamics (MD) simulations. The snapshots reveal the specific locations on the protein where ligands can interact. Based on Figure 8, Hyoscyamine molecule is surrounded by residues HIS41, THR45, SER46, CYS145, MET165, PHE181, MET49, THR25, VAL186, GLN189, ASP187, ARG188, CYS44, HIS164, and GLU166.

Figure 9 shows that in the equilibrated structure of Moupinamide-M^pro, amino acids LEU27, HIS172, HIS41, GLY143, MET165, LEU50, PRO52, GLN189, LEU141, SER144, HIS163, ARG188,

Figure 9: These images depict the stable arrangement of the Moupinamide-Mpro complex and the amino acids that interact with the Moupinamide molecule

HIS164, ASN142, MET49, PHE140, PHE181, ASP187, ASN51, CYS145, PRO39, and GLU166 surround the ligand.

Finally, Hydroxytropane is interacting with amino acids GLU166, SER46, HIS164, ASN142, GLU166, MET49, HIS163, PHE140, MET165, HIS41, GLN189, ASP187, CYS145, HIS172, ARG188, SER144, and LEU141. From the equilibrium snapshots, one can conclude that all the ligands adopt a similar preferential binding site.

To examine the distances of interaction between Mpro and the ligand molecules throughout the MD simulations, we generated radial distribution function (RDF) plots for HIS41 and CYS145 in relation to the ligands. These plots are depicted in Figure 11. The RDF plots are used to calculate the likelihood of B atoms being distributed around A atoms as the reference atoms. The equation below represents the function gA-B(r) that describes this relationship:

Equation 1 illustrates the distribution of B atoms around A atoms within a spherical shell with a thickness of Δr. The nB value represents the number of B atoms present in this shell, while NB indicates the total number of B atoms in the amorphous cell. The radial distribution function (RDF) analysis was performed to investigate the interaction between ligands and the HIS41 residue of Mpro. The RDF profiles of Ligand-HIS41, Moupinamide-HIS41, and Hydroxytropane-HIS41 were obtained from the molecular dynamics (MD) simulation. The RDF profile of Ligand-HIS41 showed that Hyoscyamine is predominantly located at a distance of 0.50-0.75 nm from the HIS41 residue. In the case of Moupinamide, there was a sharp peak indicating that the ligand maintains a distance of 0.35-0.45 nm from HIS41 throughout the MD simulation. The RDF profile of Hydroxytropane exhibited a constant distribution zone starting from 0.40 nm and extending up to 0.70 nm, suggesting that the Hydroxytropane-HIS41 interaction distance fluctuates at different MD frames. The sharpest Ligand-CYS145 RDF peak is observed to appear at 0.50-0.70 nm for Hydroxytropane, while Hyoscyamine shows a relatively broad peak distributed in the 0.50-0.90 nm. In the case of Moupinamide, a bimodal distribution is seen to take place including 0.30-0.40 nm and 0.50-0.70 nm. In general, radial pair distribution functions revealed a relatively constant interaction distance between ligands and M^pro amino acids.

Figure 10: These images depict the stable arrangement of the Hydroxytropane-Mpro complex and the amino acids that interact with the Hydroxytropane molecule

Figure 11: The analysis of ligand molecules for radial pair distribution profiles with the HIS41 and CYS145 amino

The radius of gyration (R_g) for a protein with N atom is calculated as follows:

To assess protein compactness, the radius of gyration (Rg) was chosen as an analysis tool. Rg is computed by using a formula that involves various variables. The protein's total mass is represented by M, while N represents the total number of atoms in the protein. The mass of the kth atom is denoted as mk, and rmean represents the central position of the protein. The distance between the kth atom and the center is represented as (rk-rmean). One benefit of utilizing Rg is that it enables us to assess the stability of protein folding. A protein that is folded in a stable manner will consistently display a similar Rg value during the MD trajectory. Conversely, an unfolded protein will show variations in Rg [26]. Figure 12 shows the radius of gyration (Rg) for unloaded Mpro and its complexes with Hyoscyamine, Moupinamide, and Hydroxytropane as a function of simulation time. The average Rg values during the last 10 ns of the Molecular dynamics trajectory are approximately 2.21 ± 0.01, 2.21 ± 0.01, 2.24 ± 0.01, and 2.23 ± 0.01 nm for free Mpro, Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro, respectively. The Rg value of Mpro remains relatively constant when interacting with Hyoscyamine, while Moupinamide and Hydroxytropane slightly decrease the compactness of Mpro.

Figure 12: The Rg values of Mpro without any load, as well as the complexes formed with Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro

Hydrogen bond analysis

Hydrogen bonds play a crucial role in biological systems as significant noncovalent interactions. These hydrogen bonds form when a hydrogen atom from a group X―H, where X represents an electronegative atom, interacts with one or more electronegative atoms (Y) that possess a pair of non-bonded electrons [39]. Hydrogen bonds between the amine proton atoms and carbonyl oxygen atoms of the peptide backbone are responsible for the formation of secondary structures in proteins, such as α-helices and β-sheets. In an α-helix, a hydrogen bond is established between the carbonyl oxygen and the hydrogen of N―H that is four residues away on the same strand. On the other hand, in a β-sheet, hydrogen bonds are formed between the relevant functional groups of distinct protein strands. Thus, hydrogen bonds are crucial interactions in biological systems, playing a significant role in determining the stability of host-guest complexes. To investigate the arrangement of protein-ligand complexes and comprehend the influence of hydrogen bonds on their stability, we employed the gmx h bond analysis tool for the examination of hydrogen bonding. Figure 13 indicates the quantity of hydrogen bonds formed between Mpro and substrate molecules during a 50 ns MD simulation. From Figure 13, it is evident that Moupinamide consistently forms hydrogen bonds with Mpro, whereas the intermolecular hydrogen bonds of Hyoscyamine-Mpro and Hydroxytropane-Mpro frequently break many times during the simulation run.

Table 2 presents a summary of the key donor and acceptor groups involved in the interaction between different species. As shown in Table 1, the N―H groups of residues GLN189, ASN142, SER144, CYS145, HIS163, GLU166, and GLY143 act as hydrogen bond donors, forming interactions with the oxygen and nitrogen atoms of Hyoscyamine. Conversely, the O―H group of Hyoscyamine establishes intermolecular hydrogen bonds with the oxygen atoms of GLN189, GLU166, and HIS164 residues. In the Moupinamide-Mpro complex, the most significant donor groups of the protein are the N―H groups of residues GLN192, ASN142, GLY143, GLN189, HIS163, HIS172, and HIS41, which establish intermolecular hydrogen bonds with the ligand's oxygen atoms. In addition, the N―H and O―H groups of Moupinamide participate in hydrogen bond interactions with amino acids ASN142, HIS164, GLU166, THR190, PHE140, LEU141, SER144, and GLU166. Finally, the amine groups of HIS41, GLY143, HIS163, GLU166, and HIS172 residues create hydrogen bonds with the oxygen atoms of Hydroxytropane. The only hydrogen bond in the interaction between Hydroxytropane and Mpro where the ligand's O―H group acts as a donor is formed with the amino acid HIS164. On average, there were approximately 0.8 hydrogen bonds formed between Mpro and Hyoscyamine, 1.5 hydrogen bonds formed between Mpro and Moupinamide, and 0.8 hydrogen bonds formed between Mpro and Hydroxytropane. The results of this study are consistent with the binding score values calculated by molecular docking. These results suggest that Moupinamide interacts more strongly with Mpro than Hyoscyamine or Hydroxytropane.

Solvent accessible surface area

The solvation behaviour of a protein is a relevant consideration in evaluation of protein-ligand complex microstructure. The solvent accessible surface area (SASA) will tell to what extent protein building blocks have access to aqueous environment. Thus, SASA analysis was performed to evaluate the ability of M^pro to perform chemistry with solvent molecules and to interact with Hyoscyamine, Moupinamide, and Hydroxytropane (Figure 14). The average SASA values of protein during the last 10 ns of MD simulations are calculated to be about 152 ± 2.5, 148 ± 2.1, 150 ± 2.5, and 155 ± 2.8 nm² for M^pro, Hyoscyamine-M^pro, Moupinamide-M^pro, and Hydroxytropane-M^pro, respectively, showing that the complexation process does not affect the densification of M^pro, significantly.

Table 2. The primary hydrogen bond donor and acceptor groups in the complexes of Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro

Figure 13: The number of intermolecular hydrogen bonds formed between bound (a) Hyoscyamine, (b) Moupinamide, and (c) Hydroxytropane molecules and M^pro taken during the MD simulations

Figure 14: The surface area that is accessible to solvent molecules is measured over time in simulations of free Mpro and Mpro bound to Hyoscyamine, Moupinamide, and Hydroxytropane

In terms of solvation free energies, the estimated values for free Mpro and the complexes of Mpro with Hyoscyamine, Moupinamide, and Hydroxytropane were -25.0 ± 4.8 kJ.mol^-1.nm^-2, -28.6 ± 4.6 kJ.mol^-1.nm^-2, -26.1 ± 4.5 kJ.mol^-1.nm^-2, and -25.2 ± 5.1 kJ.mol^-1.nm^-2, respectively. Figure 15 presents the average area per residue throughout the simulation for the different systems studied. It demonstrates that all systems display a comparable pattern. A closer look at the charts shows that PRO9, TYR118, ASN119, TYR154, GLU178, LEU286, VAL297, ARG298, and SER301 of free M^pro are more exposed to solvent molecules compared with the Hyoscyamine-M^pro complex. In contrast, the residues ARG4, LYS137, ILE213, ASN214, PHE294, and VAL303 have a less access to aqueous solvent in free M^pro. In the case of Moupinamide-M^pro complex, ARG4, MET49, LYS137, ASN142, and ARG298 are more exposed to solvent molecules compared to the free M^pro, while LYS12, THR26, LEU50, ASN51, TYR118, ASN119, ARG217, and CYS300 have less access to solvent molecules. The solvent accessible surface area per residues of M^pro in Hydroxytropane-M^pro shows the least change compared to the free M^pro. Finally, the TYR118, ARG217, and LEU286 amino acids of free M^pro are more exposed to water molecules. Also, ARG4, LYS137, ASN214, PHE294, and VAL303 are the residues that are less exposed to solvent in free M^pro.

Figure 15: The surface area of the protein that is accessible to the solvent over time in simulations for the Mpro protein alone and in complexes with Hyoscyamine, Moupinamide, and Hydroxytropane.

Table 3: During the MD simulations, the average protein-ligand interaction energies (expressed in kJ.mol-1) were calculated

Protein-ligand interaction energies

The strength of the interaction between the protein and ligand was determined by calculating the short-range non-bonded interaction energies. Table 3 presents the average values for the total interaction energies, Lennard-Jones (LJ) and short-range Coulombic. Both types of short-range interaction energies follow the trend: Moupinamide > Hydroxytropane > Hyoscyamine. The calculated interaction energies suggest that Moupinamide forms a significantly more stable complex with Mpro compared to Hyoscyamine and Hydroxytropane. This finding is consistent with the analysis of hydrogen bonds.

Conclusion

Currently, there is not confirmed antiviral drug for cure of Sars-CoV-2 Virus. Antiviral potency of alkaloid molecules against main protease (3CLpro) of COVID-19 were investigated by molecular docking analysis. Among the five compounds studied, Moupinamide demonstrated the highest potency against COVID-19 Mpro with a binding energy of -7 Kcal/mol, which can be attributed to its ability to form a larger number of hydrogen bonds. To investigate the dynamic behavior of Mpro: ligand complexes, the researchers utilized atomistic molecular dynamics (MD) simulations. The calculated average RMSD values for Mpro, Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro were approximately 0.30 ± 0.02, 0.30 ± 0.01, 0.15 ± 0.02, and 0.29 ± 0.02 nm, respectively. These results suggest that the protein structure becomes significantly more stable upon binding to Moupinamide. Additionally, the mean value of RMSF was determined to be around 0.16, 0.13, 0.13, and 0.14 nm for free Mpro, Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro, respectively. The RMSF analysis suggests that binding of alkaloids reduces fluctuation of protein building blocks. Radial pair distribution functions (RDF) analysis revealed a relatively constant interaction distance between ligands and M^pro amino acids. Moupinamide forms persistent hydrogen bonds with M^pro, while Hyoscyamine-M^pro and Hydroxytropane-M^pro intermolecular hydrogen bonds break many times during the simulation run. Regarding the formation of hydrogen bonds between Mpro and the ligands, the average numbers were calculated to be 0.8, 1.5, and 0.8 for Hyoscyamine, Moupinamide, and Hydroxytropane, respectively. This finding is in good agreement with the binding score values obtained from molecular docking calculations and indicates the stronger interactions between Moupinamide and M^pro compared to Hyoscyamine and Hydroxytropane. The data reveals that the average short-range Coulombic, Lennard-Jones (LJ), and total interaction energies for Hyoscyamine-Mpro, Moupinamide-Mpro, and Hydroxytropane-Mpro are -28.5, -120.2, and -148.7 kJ.mol^-1, -55.1, -169.3, and -224.4 kJ.mol^-1, and -33.2, -138.3, and -171.5 kJ.mol^-1, respectively. These results indicate that the complex formed between Moupinamide and Mpro is significantly more stable compared to those formed between Hyoscyamine and Mpro, as well as Hydroxytropane and Mpro.

Disclosure Statement

No potential conflict of interest was reported by the authors.

Funding

This study did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Authors' contributions

All authors contributed toward data analysis, drafting, and revising the paper and agreed to responsible for all the aspects of this work.

Conflict of interest

The authors declare that they have no conflicts of interest in this article.

ORCID

Ali Ramazani

https://orcid.org/0000-0003-3072-7924

HOW TO CITE THIS ARTICLE

Neda Tadayon, Ali Ramazani*. Molecular Docking and Dynamics Analysis of COVID-19 Main Protease Interactions with Alkaloids from Hyoscyamus Niger and Datura Stramonium. Chem. Methodol., 2023, 7(11) 883-903

DOI: https://doi.org/10.48309/chemm.2023.422059.1732

URL: https://www.chemmethod.com/article_183350.html