In silico exploration of disulfide derivatives of Ferula foetida oleo-gum (Covexir®) as promising therapeutics against SARS-CoV-2.

Although vaccines have been significantly successful against coronavirus, due to the high rate of the Omicron variant spread many researchers are trying to find efficient drugs against COVID-19. Herein, we conducted a computational study to investigate the binding mechanism of four potential inhibitors (including disulfide derivatives isolated from Ferula foetida) to SARS-CoV-2 main protease. Our findings revealed that the disulfides mainly interacted with HIS41, MET49, CYS145, HIS64, MET165, and GLN189 residues of SARS-CoV-2 main protease. The binding free energy decomposition results also showed that the van der Waals (vdW) energy plays the main role in the interaction of HIS41, MET49, CYS145, HIS64, MET165, and GLN189 residues with the inhibitors. Furthermore, it is found that the Z-isomer derivatives have a stronger interaction with SARS-CoV-2, and the strongest interaction belongs to the (Z)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane (ΔG = -18.672 kcal/mol). The quantum mechanical calculations demonstrated that the second-order perturbation stabilization energy and the electron density values for MET49-ligand interactions are higher than the other residue-ligand complexes. This finding confirms the stronger interaction of this residue with the ligands.

Introduction

The first strain of human coronavirus (HCoV) was found in patients with respiratory disease in the mid-1960s [1,2]. Since then, many coronavirus (CoV) strains have been discovered that spread in human society, raising concerns about the high contagiousness and the dangerousness of the novel coronaviruses [3,4]. In December of 2019, a coronavirus infectious disease (named COVID-19) was reported by the Chinese Center for Disease Control [5]. As of Jun 1, 2020, 492,144,558 cases of COVID-19 have been confirmed globally, causing more than 6,177,552 deaths (https://www.worldometers.info/coronavirus/).

The SARS-CoV-2 is accountable for COVID-19 disease. Its structure consists mainly of four structural proteins (i.e., spike, envelope, nucleocapsid, and membrane), two polyproteins, as well as other accessory proteins [6]. The cleavage of these proteolytic into mature non-structural proteins is a crucial step for the replication and transcription of coronavirus [7,8]. The papain-like protease PLpro and the chymotrypsin-like protease 3CLpro (also known as Mpro) enzymes are involved in this conversion [7,8]. Therefore, most of the discovery projects have been launched to target specific viral proteins, mainly its main protease. The crystal structure of SARS-CoV-2 main protease was released [9]. This crystal structure gives a solid structural basis for the development of potential therapeutics that can interact with Mpro.

COVID-19 constantly changes via mutation, and sometimes these mutations result in emerging a new variant of the virus. The fast changes in COVID-19 make combating the virus even harder. Recently, World Health Organization (WHO) reported a new variant of coronavirus (named Omicron) in South Africa on November 24, 2021. Omicron is also associated with the massive number of infected people throughout the world and, consequently, the high rate of hospitalization, which again put the world on red alert.

According to WHO, currently, several vaccines are used to combat COVID-19 (https://covid19.trackvaccines.org/agency/who/). To date, the US Food and Drug Administration (FDA) has only approved remdesivir as a therapeutic antiviral treatment against COVID-19. Nevertheless, there have been many reports of the low effectiveness of this drug [10]. Therefore, numerous studies were conducted to find more efficient drugs to combat COVID-19 [11,12]. Accordingly, several computational studies have been reported to examine many potential inhibitors of SARS-CoV-2 main protease [[13], [14], [15]].

Nature is a very plentiful source of potential drugs and drug leads. Ferula foetida (F foetida) is one of the Ferula species with nutritional and medicinal applications. Traditionally the oleo-gum exudates of this plant is utilized to treat many diseases [[16], [17], [18], [19]]. F. foetida oleo-gum is a rich source of organic sulfides including thiophene, disulfides and polysulfide derivatives. The main compounds of F. foetida oleo-gum include (E)-1-sec-butyl-2-(prop-1-enyl)disulfane (E2S), (Z)-1-sec-butyl-2-(prop-1-enyl)disulfane (Z2S), (E)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane (E3S), and (Z)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane (Z3S). In addition, the oleo-gum contains a water-soluble part of oligosaccharides that typically exhibit no pharmacological activity [[20], [21], [22]]. Concerning the main disulfides of F. foetida, there is now evident that polysulfides like disulfides have the potential to release the gasotransmitter H2S at physiological levels. On the other hand, H2S releasing compounds are currently considered to have antiviral effects as well as improving viral defense of immune system. Moreover, according to an experimental and theoretical study by Wang and coworkers [23] it is was revealed that their studied disulfide compounds are excellent inhibitors of SARS CoV Mpro, It should be noted that the structural of SARS-CoV has similarities to other coronaviruses, including SARS-CoV-2, Accordingly, we performed a comprehensive in-silico study to investigate the binding mechanism of main disulfide derivatives isolated from Ferula foetida (Fig. 1 ) to SARS-CoV-2 main protease.

Fig. 1

The chemical structure of a) (E)-1-sec-butyl-2-(prop-1-enyl)disulfane, b) (Z)-1-sec-butyl-2-(prop-1-enyl)disulfane, c) (E)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane, and d) (Z)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane.

Computer-aided drug design methods emerged as efficient tools in modern drug discovery and development. Generally, these methods, according to scoring protocols, can identify and optimize the drug leads. The docking calculations score the protein-ligand complexes based on obtained molecular mechanics energy. This method can be used to explore big molecular sets, but it is not very accurate. Molecular dynamics (MD) simulation can be employed to characterize the dynamic behavior, nature of the interaction, and stability of the best-docked protein-ligand complexes. Furthermore, free energy calculation methods, such as the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method, have been widely used for evaluating the thermodynamic properties of complexes during an MD simulation.

In this work, we conducted a computational study to investigate the binding affinity of four disulfide derivatives isolated from Ferula foetida with SARS-CoV-2 main protease. The molecular docking and molecular dynamics simulation is performed to find high-scored configurations and evaluate the stability of the four drug-receptor complexes. Furthermore, free energy calculation using MM- PBSA method is employed to identify the role of hotspot residues of SARS-CoV-2 in interaction with the studied potential drugs. Finally, the quantum mechanics (QM) calculations are used to provide deeper insights into the nature of drug interactions with SARS-CoV-2 main protease at the electronic level.

Computational details

Ligand and receptor preparation

The starting crystal structure for SARS-CoV-2 main protease was taken from RCSB Protein Data Bank (PDB) (https://www.rcsb.org/, PDB ID: 6LU7) [24]. The SARS-CoV-2 structure has two chains (A and C) and chain A was selected as the target receptor. The protein was prepared for docking by removing the water molecules, the inhibitor, and the other heteroatoms from its structure (Fig. 2 ).

Fig. 2

Three-dimensional structure of main-protease (PDB ID: 6LU7).

In this work, four disulfide drugs including, (E)-1-sec-butyl-2-(prop-1-enyl)disulfane, (Z)-1-sec-butyl-2-(prop-1-enyl)disulfane, (E)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane, and (Z)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane, were selected as the drug candidates. The spectral data of these compounds are provided in supplementary material (Figures S1-S6). All of the ligand molecules were fully optimized using Gaussian 03 software [25] at the M062X/6-311G** level of theory [26]. The optimized structures of these molecules are shown in Fig. 3 .

Fig. 3

The optimized structure of a) (E)-1-sec-butyl-2-(prop-1-enyl)disulfane, b) (Z)-1-sec-butyl-2-(prop-1-enyl)disulfane, c) (E)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane, and d) (Z)-1-(1-(methylthio)propyl)-2-(prop-1-enyl)disulfane.

Molecular docking protocol

Before beginning the docking investigation, the co-crystalized ligand (N3) on its receptor is re-docked to validate the docking calculations. It is observed the root-mean-square deviation (RMSD) is about 1.9 Å which is the acceptable value. The optimized structures of the studied inhibitors are selected as ligands. The protein was prepared by adding the hydrogen atoms to its structure using Chimera software [27]. Then by the steepest descent algorithm [28], energy minimization was performed, which terminated after 1000 steps. The SARS-CoV-2 main protease binding site was determined from the previous studies to ensure the ligands exactly bind to the receptor active site [13,14]. Accordingly, the binding region is specified by defining a grid box with 22.0 × 23.0 × 22.0 Aͦ 3 dimension, which its center was in the coordinates of x = −11.14 Aͦ, y = 19.90 Aͦ, and z = 63.06 Aͦ. Four different systems are designed in which each of the inhibitors is located at the center of the grid box. The dock screening was executed using the AutoDock [29] executable plugin in chimera software. The docking tool was generated ten poses per ligand in which the ligands attached to the protein binding site. The best-docked pose for each drug candidate was extracted for further analyses.

Molecular dynamics (MD) simulation protocol

MD simulations were performed to investigate the relative stability of the ligand-receptor and generate a set of conformations for Molecular mechanics/Poisson−Boltzmann (Generalized Born) surface area (MM/PBGBSA) calculation. In MD simulation, to construct the ligand-receptor complexes (i.e., E2S-6LU7, Z2S-6LU7, E3S-6LU7, and Z3S-6LU7), we extracted the best-scored docking poses. Then, these complexes were placed separately in the center of simulation boxes. Furthermore, for comparison, another system is designed in which the apoprotein is placed in the center of a simulation box. It should be noted that the simulation boxes are large enough to avoid their interaction with neighbor boxes. Three-site transferrable intermolecular potential (TIP3P) water model [30] added to systems to mimic the solvent molecules. Na+ and Cl− ions were added for neutralizing the systems. The initial snapshot from one of the systems is given as an example in Error! Reference source not found. Based on this protocol, five systems, hereafter named as 6LU7, E2S, Z2S, E3S, and Z3S, were setup which more details about them are provided in Table S1. For all of the components of systems, the force field parameters are applied by the CHARMM36 force field [31].

In this work, the MD simulations and post-MD analyses were performed using the GROMACS package version 2021.3 [32]. Through the steepest descent algorithm, 1000 steps energy minimization was applied. After that, systems have been subjected to two separate runs in NVT and NPT ensembles for 200 ps. The temperature of systems was slowly raised to 310 K using the V-rescale (V stand for velocity) thermostat [33], and pressure reached 1 bar by Berendsen algorithm [34]. MD productions for the studied systems were run under periodic boundary conditions for 90 ns. All of the bonds constrain at their equilibrium position with the LINear Constraint Solver (LINCS) algorithm [35]. Long-range electrostatic interactions were treated by the particle-Mesh Ewald (PME) method [36] and nonbonded interactions were calculated with a 1.2 nm range cutoff.

Free energy calculation

MM/PBGBSA method [37] was used to provide more details about the interaction of the ligand with the receptor. The MM-GBSA method using molecular mechanics (MM), the Generalized Born (GB) electrostatics, and solvent accessibility (SA) models estimate the binding energy. This method is a reliable and efficient model, which is widely utilized to calculate binding free energy binding free energies of noncovalently bound complexes [13,[38], [39], [40]] Here, the gmx_MMPBSA python script is used for compute binding free energy from the MD productions. The binding free energy for ligand-receptor complexes can be computed as follows:where the , , and are the Gibbs free energy of ligand-receptor complex, receptor, and ligand, respectively. In the MMPBSA method, each of the Gibbs free energy term can be obtained from the following equation:where the correspond to the molecular mechanical energy in the gas phase, the solvation energy, temperature, and entropy, respectively. On the other hands, G is equal to the change in the enthalpy (H) minus the change in the product of the temperature times the entropy:

More details about this method can be found in Tresanco et al. work [41,42]. Furthermore, gmx_MMPBSA_ana (a graphical user interface tool) is used to analyze gmx_MMPBSA results.

Quantum mechanical calculations

Quantum mechanical calculations are performed to provide an insight into the electronic properties of the ligand-residue complexes. To build the systems for these calculations, the hotspot residues that are involved in interactions with the ligands from the outcome of binding free energy decomposition were selected. Well-established density functionals such as the B3LYP functional are a good candidate for systems in which hydrogen-bonding interactions govern the min role. However, many of the gas-phase biological systems contain significant dispersion interactions. For evaluating these systems, many functional levels such as M06–2X have been specifically developed for representing the dispersion interactions that offer a potentially powerful alternative to accurate ab initio methods. Then, single-point energy calculations are performed on the ligand-residues complexes at the M06–2X/6-311G** level of theory. The M06–2X functional reported by Truhlar's group that used in the investigation of a wide range of molecular systems [43,44]. This level is a hybrid meta functional with 54% Hartree–Fock (HF) exchange, which was designed to handle weak interactions by including dispersion correction [26]. All of the density functional theory (DFT) calculations are performed using the Gaussian 03 software. More details about DFT calculations will be discussed in ‘Results and discussion section.

Results and discussion

Docking screenings

In this study, the interaction of SARS-CoV-2 main protease with four disulfide drugs (see Fig. 3) is investigated. At first, the obtained results from the docking calculation are evaluated. The energy and configuration of the best-ranked of each ligand are given in Table 1 and Fig. 4 , respectively.

Table 1

The docking results of the selected poses for MD simulations.

Ligand	Energy (kcal/mol)	Amino acid	Type of interaction	Interaction distance (Å)
E2S	−3.6	HIS41	Pi-alkyl	3.80
			Pi-alkyl	4.06
		MET165	H-bond	2.98
			Alkyl	4.25
Z2S	−3.7	HIS41	Pi-alkyl	4.09
			Pi-alkyl	4.37
			Pi-alkyl	5.16
		MET49	H-bond	3.72
		CYS145	Alkyl	3.09
			Alkyl	5.22
		MET165	H-bond	3.16
E3S	−3.7	HIS41	Pi-alkyl	4.00
			Pi-alkyl	5.34
		MET49	H-bond	3.89
		CYS145	H-bond	2.84
		HIS163	Pi-alkyl	4.70
		MET165	H-bond	5.05
Z3S	−3.7	HIS41	Pi-alkyl	3.93
			Pi-alkyl	4.65
			Pi-alkyl	4.73
		MET49	H-bond	3.42
			Alkyl	5.20
		CYS145	H-bond	3.27
		MET165	H-bond	3.04

Fig. 4

Predicted binding modes of four SARSCoV-2-ligand systems:A) E2S, B)Z2S, C)E3S, and D) Z3S. (I) and (II) stands for 3D structure and 2D structure, respectively.

The obtained energy results show that the interactions of the studied inhibitors are weaker than the reported data for Remdesivir [45]. It can be related to the ability of Remdesivir to form hydrogen bonding and its bigger size that interact with more residues. Keep in mind the disulfide derivatives interactions are strong enough to form stable complexes with Mpro. The interacting residues of Mpro involved in the binding of the ligands to the protein active site can be observed in Fig. 4. As can be seen in this Figure panel A, E2S molecules interact with HIS41 and MET165. In the Z2S case, the ligand-receptor interaction becomes stronger in which four intermolecular interactions between the Z2S molecule and HIS41, MET49, CYS145, and MET65 are formed (Fig. 4B). In two other systems, E3S and Z3S molecules form more intermolecular interactions with the binding site (E3S interact with HIS41, MET49, CYS145, HIS163, and MET165 sites and Z3S interact with HIS41, MET49, CYS145, and MET65). In these interactions, amino acids that have a sulfur group or an imidazole ring are involved in the formation of interactions with the ligands. According to the obtained results, it can be concluded that the hydrophobic interactions have the main role in the binding of the studied ligands to the active site of SARS-CoV-2 main protease.

MD simulation

The stability receptor-ligand complexes are evaluated by calculating the RMSD of the protein backbone during the MD simulation time. For four investigated systems, the RMSD values are calculated during 90 ns time and plotted as a function of the simulated time (Fig. 5 ).

Fig. 5

The calculated RSMD for protein during simulation time.

As can be seen in Fig. 5, the RMSD plots for all of the investigated systems reached their equilibrium state after ∼5 ns and fluctuated around 0.15–0.25 nm. This finding shows all of the studied complexes were almost stable during MD simulations time. It should be noted that most fluctuation belongs to the E2S system, where the RMSD value after about 25 ns increased from ∼0.22 nm to ∼0.27 nm. This observation can be related to the lowest number of the formed interaction of E2S with SARS-CoV-2 main protease. Furthermore, the comparison of docking and MD outcomes indicated that during MD simulation, the change in interacting residues with E2S and Z3S is more significant than two other inhibitors. This behavior leads to more fluctuation in the RMSD curves of E2S and Z3S systems.

The radius of gyration (Rg) can be used for measuring structure compactness and understanding the folding properties of the protein and receptor-ligand complexes. Furthermore, it is helpful to discover the effect of a drug molecule on conformational changes of the protein. The Rg is changed according to the folding state of the protein complex. A low Rg value shows the protein structure is more compact and vice versa. The variation of Rg as a function of simulation time for the studied receptor-ligand complexes are depicted in Fig. 6 .

Fig. 6

Radius of gyration (Rg) of the investigated complexes as function of simulation time.

It is observed that the Rg values for all of the studied systems (except the E2S system) fluctuated around 2.22–2.23 nm. Nevertheless, in the E2S system, the Rg curve decreased to ∼0.18 after 25 ns, while after 40 ns, the Rg returns to its average value. These obtained results indicated that the SARS-CoV-2 main protease protein in complex with Z2S, E3S, and Z3S structure retains its compactness. The E2S-6LU7 complex shows an initially loose packaging system for the protein after ∼25 ns, but after 40 ns regains its structural compactness.

To further study the effect of the investigated ligands on the Mpro, secondary structure changes upon ligand binding as a function of the simulation time are analyzed. The Dictionary Secondary Structure of Proteins (DSSP) analysis [46] is performed to determine the secondary structures of each residue of SARS-CoV-2 main protease. Fig. 7 panel A represents the percentage of protein secondary structure for 6LU7 and the ligand-6LU7 complexes.

Fig. 7

A) The percentage of protein secondary for 6LU7 and ligand-6LU7 complexes. Number of residues in various secondary structures for B) 6LU7, C)E2S-6LU7, D) Z2S-6LU7, E) E3S-6LU7, and F) Z3S-6LU7.

As can be seen in this Figure, the distribution of secondary structure for alpha-helix, beta-sheet, and beta-bridge elements are intact in free and complex forms of 6LU7. While, in comparison to free 6LU7, the coil and turn percentage decrease and increase in the complex forms, respectively. It should be noted that the bend is slightly decreased in the case of the 6LU7-Z2S. This change in the secondary structure of 6LU7-Z2S may be responsible for fluctuation in the previous plots (i.e., RMSD and Rg). The variation of the average numbers of residues for the 6LU7-E2S complex as a function of time indicated that the beta-bridge trend is almost constant over time. The coil and bend fluctuated during the simulation and decreased at the end of 90 ns. The rest elements also fluctuate over time but eventually have an upward trend. The other ligand-receptor complexes behave similarly in which the beta-sheet, alpha-helix, and bend slightly increased while the coil and turn decreased (see Fig. 7).

Binding free energy calculation

The endpoint MM-PBSA method is used to evaluate the binding affinity of the studied ligand to the SARS-Cov-2 main protease. It is known that due to the high flexibility of ligand-6LU7 complexes, the entropy term must be included in binding free energy calculations. The obtained binding free energy results for the investigated system are given in Table 2 and Figure S8. The negative values of enthalpy and Gibbs free energy reveal that the ligand-binding is exothermic, and the formed complexes are thermodynamically stable. As can be seen in Table 2, Z-isomers have a stronger bond to the 6LU7, and the highest Gibbs free energy belongs to Z3S (−18.672 kcal/mol). It is observed that the entropy contribution in Gibbs free energy is considerable so that neglect of it can change the order of complexes stability.

Table 2

The obtained results from MM-PBSA binding free energy calculations (all energy term in kcal/mol).

	E2S	Z2S	E3S	Z3S
ΔH	−19.624 (±3.01)	−17.140 (±2.04)	−21.513 (±2.30)	−21.185 (±1.42)
-TΔS	8.723 (±0.56)	2.500 (±0.03)	2.870 (±0.04)	2.512 (±0.08)
ΔG	−10.900 (±3.01)	−14.639 (±2.04)	−18.643 (±2.30)	−18.672 (±1.42)

Interestingly, the most change in entropy value belongs to the E2S complex (8.723 kcal/mol). Keep in mind that most of the fluctuation in the MD results has been related to this complex. This finding can be attributed to the reorientation of the E2S ligand in the binding site. The MM-GBSA binding free energy decomposition is performed to recognize the hotspot residues of the SARS-CoV-2 main protease ligand binding. This calculation is helpful to design potent and selective inhibitors for specific purposes. The obtained results from the per residue binding energy decomposition analysis suggested HIS41, MET49, CYS145, HIS64, MET165, and GLN189 have major contributions in the binding of the studied ligands. The binding energy values for all the interacting residues with investigated complexes are given in Figure S9. The three most important interacting hotspot residues for each complex are shown in Fig. 8, Fig. 9 . Furthermore, the decomposition of energy terms for these residues is illustrated in Figure S10.

Fig. 8

Close-up snapshot from interaction region obtained from the MD simulation where the ligand interacts with residues in A) E2S and B) Z2S systems.

Fig. 9

Close-up snapshot from interaction region obtained from the MD simulation where the ligand interacts with residues in A)E3S and B) Z3S systems.

As can be seen in Fig. 8 panel A, E2S mainly interacts with HIS41, MET49, and GLN189 residues. It is found that the MET49 residue has the highest contribution in ligand binding. The van der Waals (vdW) energy term plays the main role in the interaction of hot spot residues with E2S (Figure S10). It should be noted that a small amount of electrostatic energy also participates in the interaction of HIS41 with E2S. Z2S bond to MET49, MET165, and GLN189 residues (Fig. 8B). Close inspection of Fig. 8B indicated that the interaction hot spot residues with Z2S are weaker than those interactions with E2S. It is noteworthy that the van der Waals energy still dominates the interaction.

In the case of E3S, the interaction of the ligand with CYS145, MET165, and HIS164 is stronger than the other residues (see Fig. 9A and S10). Furthermore, it is found that the interactions of E3S with CYS145 and MET165 have a hydrophobic nature and are mainly due to vdW interaction. However, in E3S–HIS164 interaction, the electrostatic energy term is predominant (Figure S10). As mentioned above, the Z3S molecule has the highest binding affinity to SARS-CoV-2 main protease. This high affinity is mainly because of the formation of strong interaction between the ligand and HIS41, MET49, and GLN189 residues. MET49 has the highest Gibbs free energy, where the vdW term governs the main role.

QM calculations

According to obtained results from the binding free energy decomposition (see Figure S9) for each ligand, a complex is designed in which the ligand interacted with hotspot residues. It should be noted that these complexes are extracted from the last frame of MD simulation. The designed complexes for the studied systems are given in Figure S11. A reduced density gradient (RDG) analysis is carried out on these systems using the Multiwfn program [47], and obtained results are given in Fig. 10 . In the RDG plot, the blue, green, and red regions imply strong attractive interactions, weak attractive interactions, and strong steric effects, respectively. In this study, more intermolecular interactions are observed in the green area. It is found that the intermolecular interactions in Z-isomers more the E isomers, as well as the Z3S system has the most interactions with the selected residues. This result is in line with the binding free energy outcomes.

Fig. 10

The reduced density gradient (RDG) vs sign (λ2)ρ plots for A) E2S, B) Z2S, C) E3S, and D) Z3S system.

Furthermore, natural bond orbital (NBO) and atom in molecule (AIM) theories are employed to obtain the orbital interactions and electronics properties of the investigated complexes. The obtained results of NBO and AIM analyses are summarized in Table 3 . For each residue-ligand interaction, the donor and acceptor orbitals with the highest E(2) (The second-order perturbation stabilization energy) are listed. It is found that in most cases the charge is transferred from the residues to ligands. Furthermore, when the MET49 is involved in the interaction, the energy values of E(2) and the electron density are higher, indicating the stronger interaction of this residue with the ligand.

Table 3

The second-order perturbation energy (E(2)) and electron density corresponds to the charge transfer between the hotspot residues and the studied ligands.

Ligand (L)	Residue (R)	Donor	Acceptor	E(2) (kcal/mol)	ρ(a.u.)
E2S	HIS41	π N - C (R)	σ* C – H (L)	0.24	0.006237
	MET49	LP S (R)	σ* C – H (L)	0.58	0.007071
Z2S	MET49	LP S (L)	σ* C – H (R)	1.98	0.010934
	MET165	LP S (R)	σ* C – H (L)	0.72	0.006008
	GLN189	LP O (R)	σ* C – H (L)	1.14	0.009420
E3S	CYS145	LP S (R)	σ* C – C (L)	0.17	0.005322
	MET165	σ C – H (R)	σ* C – H (L)	0.08	0.005259
Z3S	MET49	π C - C (L)	σ* C – H (R)	0.25	0.005973
	GLN189	σ C – H (R)	σ* C – H (L)	0.25	0.005913

Conclusion

In summary, the binding mechanism of four disulfide compounds from Ferula foetida to SARS-CoV-2 main protease is investigated. By using molecular docking, the high-scored configurations for all drugs are identified. The MD simulation results indicated the receptor-ligand complexes are almost stable. The HIS41, MET49, CYS145, HIS64, MET165, and GLN189 residues of SARS-CoV-2 main protease have the most interactions with the studied inhibitors. Furthermore, it is found the interaction of Z3S with SARS-CoV-2 main protease is relatively stronger than the other inhibitors. This inhibitor formed a stable complex with SARS-CoV-2 main protease through the vdW interactions. MET49, due to higher values of E(2) and the electron density, has a stronger interaction with the inhibitors.

Data availability

All of Computational models are available from the corresponding author upon request.

Note

The authors declare no competing financial interest.

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.