Molecular Docking and Dynamic Simulation Revealed the Potential Inhibitory Activity of Opioid Compounds Targeting the Main Protease of SARS-CoV-2

Opioids are a class of chemicals, naturally occurring in the opium poppy plant, and act on the brain to cause a range of impacts, notably analgesic and anti-inflammatory actions. Moreover, an overview was taken in consideration for SARS-CoV-2 incidence and complications, as well as the medicinal uses of opioids were discussed being a safe analgesic and anti-inflammatory drug in a specific dose. Also, our article focused on utilization of opioids in the medication of SARS-CoV-2. Therefore, the major objective of this study was to investigate the antiviral effect of opioids throughout an in silico study by molecular docking study to fifteen opioid compounds against SARS-CoV-2 main protease (PDB ID 6LU7, Mpro). The docking results revealed that opioid complexes potentially inhibit the Mpro active site and exhibiting binding energy (-11.0 kcal/mol), which is comparably higher than the ligand. Furthermore, ADMET prediction indicated that all the tested compounds have good oral absorption and bioavailability and can transport via biological membranes. Finally, Mpro-pholcodine complex was subjected to five MD (RMSD, RMSF, SASA, Rg, and hydrogen bonding) and two MM-PBSA, and conformational change studies, for 100 ns, confirmed the stability of pholcodine, as a representative example, inside the active site of Mpro.


Introduction
Multiple instances of pneumonia with an unclear origin were reported in December 2019. A WHO-designated pathogen coronavirus illness was identified in 2019 (COVID-19) [1].
Laterally, it was renamed severe acute respiratory syndrome coronavirus two (SARS-CoV-2), which belonged to a family called Coronaviridae [2]. SARS-CoV-2 is a singlestranded RNA virus that is enclosed and has a positive sense (+ssRNA) [1].

Complications of COVID-19.
Although there is no evidence of direct insult to the CNS, the virus was not detected in most CSF examinations [3]. No direct evidence of the virus in the brain has been discovered during postmortem examinations [4]. Individuals with COVID-19, in contrast, had a considerably higher rate of ischemic strokes than those with influenza [5]. There are not enough high-quality cohort studies or case series to get the whole picture.
As a complication of the neuromuscular system, death was reported in one of the 11 patients [6]. One of the fatal complications is thrombus formation due to SARSpropensity CoV-2 infecting endothelial cells via ACE-2 (angiotensin-converting enzyme 2) [7].
Morphine has been proven in both in vitro and in vivo animal studies to have various properties such as antiinflammatory, antifibrotic, anticancer, cardioprotective, and renoprotective [15,16]. Opioids can antagonize the impact of angiotensin-converting enzyme 2 (ACE2) [17], which results in the renin-angiotensin system being dysregulated (RAS). Additionally, opioids have been found to inhibit COVID-19 pathogenesis via cytokine production and inflammatory cell infiltration in the lungs during various viral infections by their immunomodulatory impact [18].

Opioid Family.
Morphine (2), as a powerful analgesic and natural narcotic compound, is considered a powerful agonist of the u-opioid receptor with high abuse potential. Compound 2 is absorbed orally with a median time to the maximum blood concentration of 0.75 h [19]. According to recent studies, heroin (15) is a more potent analgesic than its more active metabolite morphine (2) and 6acetylmorphine (3). Deacetylation of compound 15 to compounds 3 and 2 leads to questions about the exact cause responsible for the potency difference [20]. Compound 15 (diacetylmorphine) is a semisynthetic derivative of compound 2. Compound 15 is lipid-soluble and absorbed rapidly after parenteral administration [21]. Compound 15 is biotransformed to compound 3 and is much slower than compound 2 by blood and various tissue, including the brain [22,23]. Compared to compound 2, compound 15 has a greater water solubility [20,24], is the fastest onset of action [20], produces a greater degree of euphoria, and has fewer side effects. For example, compound 15 is approximately 2 to 16 times more potent than compound 2 in producing reinforcing effects in animals and subjective effects [25,26]. We have discovered that compound 2, but not compound 15, binds to opiate receptors, suggesting that compound 15 serves primarily as a lipid-soluble prodrug for compound 2's central distribution. The discovery was that compound 2 levels in the brain are higher after compound 15 administration compared to a comparable dose of compound 2. The 6-acetylmorphine possessed intrinsic activity and triggered several opiate-like effects [27,28]. Studies suggested that compound 15 is rapidly hydrolyzed to its 6-acetyl derivative, and compound 3 has opiate receptor affinity, whereas compound 15 does not. As a result, compound 15 increased potency compared to compound 2 [20]. Codeine (6), also known as 3-methylmorphine, is a mild opioid with analgesic and antitussive properties [29]. Compound 6 is conversed to compound 2 by the cytochrome P450 enzyme, CYP2D6, which is responsible for its analgesic effect. Compound 6 also has some (low) affinity for the u-opioid receptor found in the central nervous system (CNS) and peripheral tissues such as the gastrointestinal tract [30]. When used as directed, lowdose compound 6 in fixed combinations with other drugs is effective and safe [31]. In animal test systems, pholcodine has antitussive activity comparable to compound 6 [32]. Normorphine (5) is an opiate analog that was first described in the 1950s as an N-demethylated derivative of morphine. Compound 6 is a less potent analgesic than compound 2. An amount of 30 mg is indicated as a human dose that can produce less sedation, miosis, vomiting, and respiratory depression than an equal dose of compound 2 [33]. Meperidine (pethidine) (11) is a phenyl piperidine derivative that plays an opioid receptor agonist role. In the United States, meperidine is marketed under the brand name Demerol. Due to concerns regarding adverse effects, pharmacological interactions, and the neurotoxicity of normeperidine (its metabolite), several doctors are avoiding using this medication as an initial-line opioid analgesic [34]. Tramadol (12) is a racemic mixture of tramadol R (+) and tramadol S (-). Additionally, it regulates the monoaminergic system, unlike typical opioids, by reducing noradrenergic and serotoninergic reuptake [35]. As a result, tramadol is classified as an atypical opioid. Compound 12 is one of the most often recommended analgesics for moderate to severe pain owing to its special pharmacological features [36]. Compound 12 was established in Germany during the 1970s and received FDA approval in 1995 but was reclassified as a schedule IV drug in 2014 [34]. Other opioids are also used to treat pain. Other opioid analgesics and cough suppressants are also used to relieve pain. Ethylmorphine (7) is an opioid analgesic and cough suppressor. Norethylmorphine is demethylated to norethylmorphine (catalyzed by CYP3A4) and then Odeethylated to provide compound 2 (catalyzed by CYP2D6). Dihydromorphine (4) has been considered to be pharmacologically more effective because of its large selectivity for opioid receptors [19]. Levorphanol (13) is a one-of-a-type synthetic opioid due to its varied actions as an agonist for both the opioid and d-and k-opioid receptors. Compound 13 is also an antagonist of the N-methyl-D-aspartate receptor and a norepinephrine and serotonin reuptake inhibitor. The analgesic impact lasts between 6 and 15 hours [19]. 7,8-Didehydro-4,5-epoxy-17-methylmorphinan-3,6-diol (1) is derivative from epoxymorphinan which has sharp analgesic effect [37]. Compound (1) has general depressing activity consequently intended to be used in interfere with seizure threshold [38]. Based on the previously mentioned aspects, we encouraged to test the target 15 opioid compounds computationally as anti-SARS-CoV-2 candidates using molecular docking calculations as well as dynamic simulation for the most active candidates in addition to in silico ADMET prediction studies (Scheme 1).

BioMed Research International
applied to all conformers, and all minimizations were carried out using MOE with the MMFF94X force field up to an RMSD gradient of 0.01 kcal/mol and RMS (root mean square) distance of 0.1. Partial charges were then automatically generated. Molecular Database (MDB) file was used to store the collected database utilized when making docking calculations.
2.1.1. Optimization of M pro . The protein data bank provided the X-ray crystallographic structure of the binding site of M pro (PDB: 6LU7). The chemicals were docked to the target enzyme's active site.

Preparation of M pro .
Delete the cocrystallized ligand. The system was then filled with conventional geometry hydrogen atoms. Automatic correction was used to check for flaws in the atoms' connection and type. The receptor's choice and atom's potential were fixed.

Docking of the 15
Molecules to M pro Active Site. The MOE-Dock software was used to dock the target molecules. In general, the following methods were used.
The Dock tool was launched after loading the enzyme active site file. The program's requirements were changed to include the following: (i) Dummy atoms as the docking site (ii) Triangle Matcher will be utilized as the placement approach (iii) London dG was chosen as the scoring mechanism, and its default parameters were set (a) Dock calculations were automatically performed after loading the MDB file of the ligand that needed to be docked (b) After studying the acquired poses, the poses that best represented the interactions between the ligand and the enzyme were chosen and saved for energy calculations 2.2. MD Simulations. CHARMM-GUI web interface and CHARMM36 force field were used to prepare the M propholcodine complex. The NAMD 2.13 package was used for all of the simulations. The periodic boundary conditions were set with a dimension of certain dimensions in x, y, and z, respectively, and the TIP3P explicit solvation model was utilized. The CHARMM general force field was used to produce the parameters for the best docking findings. After that, (Cl-/Na+) ions were used to neutralise the system. Production, equilibration, and minimization were all part of the MD protocols. All MD simulations used a 2 fs time step of integration, with the canonical (NVT) ensemble used for equilibration and the isothermal-isobaric (NPT) ensemble used for production. The pressure was maintained at 1 atm using a Nose-Hoover Langevin piston barostat with a Langevin piston decay of 0.05 ps and a period of 0.1 ps throughout the 100 ns of MD generation. The Langevin thermostat was used to set the temperature at 298.15 K. The Lennard-Jones interactions were smoothly trimmed at 8.0, and a distance cut-off of 12.0 was applied to short-range nonbonded interactions with a pair list distance of 16. The particle-mesh Ewald (PME) approach was utilized to handle long-range electrostatic interactions, with a grid spacing of 1.0 being applied to every simulation cell. The SHAKE method was used to restrict all hydrogen atom covalent bonds. We used the same protocol for all MD simulations in order to maintain consistency.

Physicochemical Properties and Lipophilicity.
A wide range of cheminformatics tools supporting the manipulation and processing of molecules are available from SwissADME and Molinspiration, including the conversion of SMILES and SD files, normalisation of molecules, production of tautomers, molecule fragmentation, calculation of various molecular properties required in QSAR and drug design, high-quality molecule depiction, and molecular database tools supporting substructure and similarity searches. Additionally, these tools allow data visualisation, bioactivity prediction, and fragment-based virtual screening. Because Molinspiration tools are designed in Java, they can essentially run on any platform. In order to eliminate structures with unsuitable properties for drugs and choose promising drug candidates, calculated molecular descriptors may be utilized for property-based virtual screening of vast collections of molecules. The following molecular characteristics were determined using Molinspiration and SwissADME. The Lipinski rule of five, which states that any compound considered to be a drug should have a partition coefficient less than 5, a polar surface area within 140 2, an H-bond acceptor less than 10, an H-bond donor less than 5, and a molecular weight within 500 dalton, is the foundation for the calculation of these properties.
2.5. ADME Data of Tested Compounds. Using the SwissADME software, the ADMET descriptors (absorption, distribution, metabolism, excretion, and toxicity) of the opioid derivatives 1-15 were identified. The studied compounds were first produced and minimized in accordance with the synthesis of small molecule methodology, after which the CHARMM force field was applied. Models for human intestinal absorption, aqueous solubility, blood-brain barrier penetration, plasma protein binding, cytochrome P450 (CYP1A2, CYP2C19, CYP2C9, and CYP3A4) inhibition, and skin permeability were among the ADMET descriptors that were included in the application.  Table 1, respectively. Most of the opioid derivatives 1-15 exerted high binding affinity to the M pro as their ΔG values range from -0.5 to -5.3 kcal/mol, compared to the reference (ΔG = −1:0 to -1.2 kcal/mol).

Results and Discussion
The docking result of reference compound 1 is completely consistent with that obtained for opioid derivatives ( Figure 4). The 2D diagram showed a crucial binding with hydrogen bonding interaction with MET165 and hydrophobic interaction with GLU166 amino acid residues.
Docking results with the M pro of opioid derivatives 1-15 revealed that most of the opioid compounds showed good binding with the M pro making several vital interactions comparing the reference and the most active compounds 3, 7, and 8 (Figures 1-3). Compounds 2, 4, 7, 8, 10, 12, 13, and 14 exhibited hydrophobic interaction with the conserved amino acid GLU166, typically as the reference. On the other hand, none of the tested compounds interacted with MET165 amino acid residue. Additionally, compounds 1, 3, 10, and 15 showed hydrogen bonding as (an H-acceptor) interaction with the HIS163 amino acid residue (Figure 1 for compound 3). Moreover, extra binding more than the reference, such as compounds 5 and 10, showed interactions within the binding site of the M pro , revealing hydrogen bonding and hydrophobic pi-H interaction with SER144 and ASN142 amino acid residues, respectively, for compounds 5, 6, and 9; however,  Figure 5(a)). The M pro -pholcodine complex expressed a minor level of fluctuation till~60 ns and stabilized later till the end of simulations (100 ns). The amino acids' flexibility of M pro was examined in terms of RMSF to figure out the protein's region that fluctuated through the process of simulation. As Figure 5(b) demonstrates, the binding of pholcodine does not make M pro flexible. The solidity and stability of the M pro -pholcodine complex were investigated by the  computation of the radius of gyration (Rg) that is inversely proportional to both solidity and stability. Figure 5(c) indicates that the Rg of the M pro -pholcodine complex at 100 ns was almost similar to that at 1 ns. The solvent accessible surface area (SASA) was computed over 100 ns to estimate the interactions between the M pro -pholcodine complex and the solvents in the media. As illustrated in Figure 5(d), the M pro -pholcodine complex featured a decrease in the surface area as the SASA values were computed to be lower at the end of the study than at the start. The hydrogen bonding 3.3. MM-PBSA. By calculating the average free binding energy from MD trajectories with a time interval of 100 ps, we were able to determine that the pholcodine has an extremely low binding free energy of -243 KJ/mol with the M pro and that the binding energy was stable over the whole period of our analysis, showing accurate binding (Figure 6(a)).
Next, we analyzed the total binding free energy of the M pro -pholcodine complex to elucidate the different parts of the binding energy and to reveal which amino acid residues played a major role in binding the ligand to the M pro to identify which amino acids had the most favorable impact in binding. Figure 6(b) illustrates that five amino acid residues of the M pro (GLU-47, ASP-48, GLU-166, ASP-187, and ASP-197) contributed more than -30 kJ/mol binding energy and are therefore considered hotspot residues in binding. Figure 7 illustrates the conformational changes that occurred because of the binding of the M pro -pholcodine complex during the 1st and 100th ns of the MD production run, indicating the incidence of some conformational changes. Additionally, the binding stability and the integrity of the complex were confirmed as pholcodine kept binding firmly to the M pro throughout the study.

Conformational Changes.
3.5. In Silico Prediction of Drug-Likeness Profiles. The design and applied new drugs are complicated because of the unacceptable ADMET parameters (distribution, excretion, absorption, metabolism, and toxicity) and the high costs for new drug development. Accordingly, it is critical to estimate the ADMET properties of a new drug [39]. Recently, the in silico ADMET analysis was applied vastly, decreasing the degradation in late production stages [40,41]. Many parameters, such as the aqueous solubility, polar surface PSA, partition coefficients, cell permeability, and intestinal absorption, have been investigated in several virtual screening studies. Lipinski's rule links the good level of oral bioavailability of a certain drug with different parameters. The molecular weight (M Wt.), Log P, hydrogen bond (HB) acceptor atoms, and HB donor atoms should be >500, >5, >10, and >5, respectively [42]. The rotatable bond's number indicates molecular flexibility that is essential in oral bioavailability. The percentage absorption (% ABS) was found to be inversely proportional to the polar surface area (PSA) measured by the equation %ABS = 109 − 0:345 tPSA [39].
Herein, we employed the software of Molinspiration [43], Molsoft [44], and SwissADME [39] to predict the ADMET characteristics of the examined opioid candidates. Table 2 shows that all compounds except compound 10 obey Lipinski's rule with Log P range values 1.09-3.36 (<5), MW range 247.33-398.50 (<500), HBD from 0 to 3 (≤5), and HBA from 1 to 6 (<10) ( Table 3). The examined compounds would have a high-level oral absorption. Also, the values of  Additionally, some other pharmacokinetic parameters were evaluated using the SwissADME software as follows:   GIT absorption level, P-gp substrate, and the inhibitory potential against several cytochrome P-450 targets. The results are listed in Table 4. The investigated opioid candidates showed medium to low ability in the skin permeability model with Log Kp ranging from -5.51 to 8. 16. Also, all compounds exerted high level of human intestinal absorption. Most of the examined compounds 2-7, 10, 11, 14, and 15 were highly bound to human glycoproteins.

Conclusion
In the global pandemic of the SARS-CoV-2, a huge need for a new treatment modality has been emerged. Our molecular docking studies revealed that the tested opioid (1-15) had a better binding affinity with the SARS-CoV-2 M pro active site (PDB ID 6LU7) than their corresponding reference and might be better alternatives to prevent SARS-CoV-2 and