INTRODUCTION

Since the earlier stages of this pandemic, a worldwide effort has been devoted to producing vaccines and antiviral drugs to combat this virus. The novel coronavirus primarily causes severe acute respiratory syndrome (SARS) by interacting with ACE-2 receptors and TMPPRS co-receptors for S protein priming [1–4]. In addition to the primary protease (3CL-pro) and a papain-like protease, it also encodes a cysteine protease, a serine protease (PLpro). These may participate in cleavage of viral polyproteins [5], leading to creation of active substances for replication of virus inside the host cell. One of the most intriguing molecular targets for the development of anti-SARS drugs is the major peptidase (Mpro) from the coronavirus (CoV) causing severe acute respiratory syndrome (SARS). In this situation, the development of natural medicinal therapies can aid in limiting the virus’ transmission [6]. Compared to traditional drug discovery methods which often take a long time, computational approaches provide a convenient and efficient way for discovering new compounds, especially important in the current pandemic. The pharmacological target 3CLpro/Mpro (PDB ID: 2GTB) is extensively researched for the development of drugs to com bat COVID-19, and is considered one of the most studied targets. This is indicated by the abundance of patents and potential drug candidates associated with [7]. It is a proteolytic enzyme required to cleave the viral polyprotein into multiple functionally active protein units. Its selection as a pharmacological target in this study is justified by the fact that its active site is completely conserved and immune to mutations. The major protease 2GTB found in the CoV associated with severe acute respiratory syndrome (SARS), reported that the major protease of 2019-nCov shows 96% similarity to that of SARS [8].

The interaction between a drug and its target is the most important step of pharmacological action. Drug placement within the receptor pocket is affected by hydrogen bonding, especially short contacts and van der Waals interactions. The proposed approach can be used as an alternative method to determine intermolecular interactions of metal-containing complexes in different systems (enzyme–ligand, reagent–substrate). Molecular regions are determined by all these interactions that ensure ligand-to-receptor complementarity [9]. Recently, various enzyme–ligand complexes (complexes of cyclin-dependent protein kinase, mouse acetylcholinesterase, HIV-1 protease, and EGFR) were studied using the AlteQ method [10–12]. Furthermore, Palko et al., 2021 proposed using AlteQ method to predict biological activity, molecular docking, and study mechanisms of drug action [13]. The AlteQ method uses high-resolution, high-quality, low-temperature X-ray diffraction to describe experimentally determined electron densities. AlteQ is an excellent tool for studying the electronic properties of large molecular structures [10–13]. This method calculates the electron density using Slater’s type atomic contributions in the interspace between the receptor and the ligand, and since their interactions are determined by the overlap of electron clouds, they follow the maximum complementarity principle, an equation can be obtained that describes these interactions. This method also evaluates the quality of the interaction between the receptor and the ligand, how complementary the interactions are, and due to this, it is used to reject less realistic structures obtained by docking methods.

Subsequently, de novo drug design is highly time-consuming (usually 10 to 15 years to get new drugs to market). In this regard, either the repurposing of FDA-approved drugs or the search for small molecules from natural consumable sources that are considered safe or have negligible toxicity have been considered appropriate approaches. Computational or in silico techniques represent an exciting approach to the world of drug repurposing. Drug repurposing is an effective approach to drug discovery [14–16]. Recently, it was reported that the CDK4/6 inhibitor palbociclib was one of his top-rated repurposing drugs to treat COVID-19. As a receptor for the development of severe acute respiratory syndrome coronavirus 2 (SARSCoV-2), the expression level of angiotensin-converting enzyme 2 (ACE2) is closely associated with SARS-CoV-2 infection [17]. Previously, our group tested baia and baib alkaloids and concluded that baib has novel structural features and is a highly potent CDK4 inhibitor [18]. However, to our knowledge, neither density functional theory (DFT) studies nor docking simulations have been reported for those compounds. Based on these findings, it seems important to fully understand the two molecules with a view to their potential applications.

In parallel, computational strategies based on density functional theory (DFT), molecular docking, and electronic databases were employed to effectively study and study the pharmacological activities of antiviral drugs against key enzymes of SARS-CoV-2 screening [19–23]. In these studies, DFT provided fundamental insights based on frontier orbital energies and spatial distributions and optimized geometries. This method has helped scientists assess drug stability and calculate structural, electronic, and thermodynamic properties [24]. These properties have been successfully used to better understand drug behavior in biological systems [25, 26]. Frontier molecular orbital energies correlate with computationally and experimentally determined properties of organic molecules screened as drug candidates [27, 28]. Rasool et al. investigated the inhibitory potency of 19 thiazolide derivatives against Mpro and methyltransferases. Arafet   et al. used a hybrid M06-2X/6-31 + \(G(d,p)\):AM1/MM method to investigate inhibition mechanisms of a known Mpro inhibitor, N3, and created energy profiles for covalent complex formation with the Mpro. Khrenova et al. performed a detailed dynamic study of the enzyme–substrate complex of Mpro. They evaluated the electronic density features of the complex and showed that QM/MM-MD trajectories disclose substrate reactivity in Mpro, and were in good agreement with relevant experimental data. Mpro is a preferred target for SARS-CoV therapeutic design, and a wide range of drug inhibitors have been developed to efficiently target it [29–31]. Previous analyzes of genome sequences have shown that SARS-CoV-2 interacts with corresponding SARS-CoV and MERS-CoV variants with a high degree of sequence similarity [32]. Therefore, Mpro is considered a highly promising biological target for SARS-CoV-2.

This article explores DFT, molecular docking, and molecular dynamics simulations with the AlteQ method to investigate the inhibitory potential of screened baia and baib compounds against the major protease (Mpro) of the SARS-CoV-2 virus. We first characterize the title compounds in terms of several reactivity parameters using DFT calculations. We then use molecular docking, molecular dynamics simulations and AlteQ together to investigate and analyze binding affinities, binding modes, and interaction stability within the active site of SARS-CoV-2 Mpro.

MATERIALS AND METHODS

Baimantuoluoamide A (baia) and baimantuoluoamide B (baib) amide alcoid compounds have been selected from our previous reports and for further detailed theoretical calculation have been studied in the following manner [18]. Utilizing Mulliken population analysis, Fukui functions were obtained by calculating the single point energies of the N, \(N - 1\), and N + 1 species of the molecule using the same basis set 6-311 + \(G(d,p)\). The GaussView 05 program was used to display the molecule’s frontier molecular orbitals and molecular electrostatic potential surface [33]. The DOS spectrum was prepared using the Gausssum 3 algorithm [34].

RESULTS AND DISCUSSION

FMO Analysis

The chemical reactivity, kinetic stability, and chemical softness of the compounds baia and baib are all characteristics that may be explained by the molecule’s HOMO–LUMO molecular orbitals. The pictorial illustration of HOMO, LUMO, HOMO-1, LUMO+1, HOMO-2 and LUMO+2 of baia and baib molecules are shown in Fig. 1. The positive and negative phase is represented in red and green color, respectively. The plots reveal that the HOMO, LUMO, HOMO-1, LUMO+1, HOMO-2 and LUMO + 2 have predicted the variation in the electron density distribution around the whole molecule. This confirms the presence of intermolecular charge transfer in the title molecules. The global hardness, chemical potential, electrophilicity index, and softness of the baia and baib molecules were determined and are shown in along with other FMOs-related molecular parameters (see Supplementary materials Table S1). In general, HOMO–LUMO related molecular properties are used to predict a drug’s tendency to engage in drug-acceptor interactions [35]. The ionization potential values were observed larger when compared to electron affinity values describing the better electron-donating capability of the investigated baia and baib molecule. The global hardness value was found to be 2.248 and 2.214 for baia and baib ligand. The global softness value was calculated as 0.445 and 0.452. The chemical hardness values for baia and baib are greater than their softness values. Therefore, higher chemical hardness denotes lower reactivity and higher stability. The chemical potential values baia (\(\mu = - 3.840\)) and baib (\(\mu = - 3.612\)) shows that baib is more stable than baia. The values shown that the baib ligand is more reactive which had strong tendency to form complexes with different metals [33, 34].

Fig. 1.
figure 1

(Color online) Frontier molecular orbitals of the (a) baia molecule and (b) baib molecule.

Full size image

MEP Analysis

Molecular electrostatic potential (MEP) is a very useful tool for the analysis of the electronic density sites by analysing electrophilic and nucleophilic reactions sites, especially for biomolecules and drugs. It generally provides essential information about the chemical stability and reactivity of an organic molecule to understand the electrophilic and nucleophilic properties. Nucleophilic and electrophilic properties are investigated and mentioned in different colors in the region. The blue color indicates higher nucleophilic potential and the red color indicates higher electrophilic potential. The electrostatic potential usually decreases in the order with blue < green < yellow < orange < red. Electrophilic regions are mainly presented around the O atom and the nucleophilic potentials are around the H and N atoms in baia and baib compounds are noticed in Figs. 2 and 3. These reactivity centers enable the compounds to bind with charged systems in the living organism which confirms bioactivity [36].

Fig. 2.
figure 2

(Color online) Molecular electrostatic potential surface of the baia molecule.

Full size image
Fig. 3.
figure 3

(Color online) Molecular electrostatic potential surface of the baib molecule.

Full size image

Mulliken Atomic Charge Distribution and Local Reactivity Descriptors

The calculated Mulliken atomic charge distribution of the molecules baia and baib is shown (see Supplementary materials Tables S2 and S3). The C17 carbon atom in the baia molecule has the largest positive charge because it has connected to the electronegative atoms oxygen (O18) and nitrogen (N19). The hydroxyl group connected to the phenyl ring’s oxygen atom (O11) was projected to have the largest negative charge value. All of the hydrogen atoms have positive charge values, while the ring carbon atoms have both positive and negative charge values, indicating that the substituents have a significant impact on the carbon atoms. Calculating Fukui functions is an efficient way to determine how reactive each particular atom in a molecule is. It identifies the preferred regions of a molecule towards the specific chemical events [37, 38]. The molecular reactivity plays an important role on the designing and synthesizing of new pharmaceutical compounds and used to realize the bioactivity of a molecule [39, 40]. The most nucleophilic sites of the baia molecule are in the order of C\(20 > {\text{C}}24 > {\text{C}}37 > {\text{C}}3 = {\text{C}}30\). The most electrophilic sites of the molecule are in the order of C\(24 > {\text{C}}36 > {\text{H}}23 > {\text{C}}37 > {\text{C}}30\). The atoms O\(18 > {\text{H}}14 > {\text{C}}6 > {\text{C}}4 > {\text{H}}32\) are favorable atomic positions for radical attack.

As in the case of baib molecule, C9 carbon atom has the largest positive charge because it attached to the electronegative atoms oxygen (O22) and nitrogen (N20). The nitrogen atom was projected to have the largest negative charge value (N20). All of the hydrogen atoms have positive charge values, while the ring carbon atoms have both positive and negative charge values, indicating that the substituents have a significant impact on the carbon atoms. The most nucleophilic sites of the baib molecule are in the order of C\(24 > {\text{O}}43 > {\text{O}}42 > {\text{C}}49 > {\text{C}}52 = {\text{C}}53\). The most electrophilic sites of the molecule are in the order of C\(24 > {\text{O}}42 = {\text{H}}64 > {\text{C}}27 = {\text{C}}53 > {\text{C}}15\). The atoms O\(23 > {\text{O}}22 > {\text{H}}8 > {\text{C}}9 > {\text{C}}10\) are favourable atomic sites for radical attack.

FMO Analysis and Density of States (DOS)

Utilizing density of states, it is possible to visualize the molecular orbital (MO) composition and its impact on chemical bonding. DOS spectrum was simulated by convoluting the molecular orbitals with Gaussian curves of unit height. In the DOS spectrum, the occupied and virtual orbitals represented by the green and red lines, respectively. Figure 4 depicts the simulated DOS spectrum of baia and baib.

Fig. 4.
figure 4

(Color online) DOS spectrum of the (a) baia and (b) baib molecule.

Full size image

NBO Analysis

The intermolecular and intramolecular interaction, which is defined in terms of stabilization energy \(E(2)\), was examined using NBO analysis. The bigger \(E(2)\) value, or more electron donating propensity from electron donors to electron acceptors, demonstrates the interaction between electron donors and electron acceptors. The donor-acceptor interactions and associated second order perturbation energies are given in Supplementary materials Tables S4 and S5. The overlap between the π(C–C) bonding orbital and the π*(C–C) anti-bonding orbital produces the hyperconjugative interactions. The intramolecular charge transfer (ICT) of the baia and baib molecules is brought on by this stabilizing interaction. Two strong LP \( \to \pi {\text{*}}\) interaction in the baia molecule arises due to the lone pairs of the nitrogen and oxygen atoms which stabilize the molecule through LP (2) O\(35 \to \pi {\text{*}}\)C33-O34 and LP (1) N\(19 \to \pi {\text{*}}\)C17-O18 with stabilization energies of 46.53 kJ/mol and 39.01 kcal/mol, respectively. In baib molecule, the strong LP \( \to \sigma {\text{*}}\) interaction arises due to the lone pairs of the oxygen atoms which stabilize the molecule through LP (2) O\(38 \to \sigma {\text{*}}\)O42-C52, LP (1) O\(47 \to \sigma {\text{*}}\)O42-C52 and LP (1) O\(47 \to \sigma {\text{*}}\) C53-H60 with stabilization energies of 105.18, 107.49 and 112.04 kcal/mol, respectively. These stabilization interactions between the LP orbitals and \(\sigma {\text{*}}\) orbitals is the characteristic feature of the biological activity of a molecule.

Analysis Data Molecular Dynamic Simulation

3.6.1. Conformational dynamics of each system: stability, flexibility, and compactness. The conformational dynamics of each system are evaluated through several variables, such as stability, flexibility, and compactness [41, 42]. The system stability analysis can be done through the convergence of total energy in each system. The results show that there is no excessive energy fluctuation in each system (see Supplementary materials, Fig. S3). Additionally, the calculation of root-mean-square deviation (RMSD) aims to determine the stability of each system in the form of complex, backbone, and ligands (see Supplementary materials, Fig. S4). In particular, the RMSD complex plays a crucial role in stabilizing the ligand and protein [42]. The results showed that the RMSD complex did not significantly fluctuate during the 40 ns simulation. It indicates that over time the interaction between the ligand and the receptor has a stable interaction. This assumption is supported by the relatively similar average value of the RMSD-complex baia: (\(0.24 \pm 0.02{\kern 1pt} \)) nm and the baib: (\(0.27 \pm 0.02{\kern 1pt} \)) nm.

In general, RMSF results showed that the ligand target is sufficiently stable, and the conformation has low flexibility. The system flexibility analysis aims to see the flexibility of the protein structure during the simulation time [43]. This variable can be calculated by the root-mean-square fluctuation (RMSF) (Fig. 5). The results show the flexibility of the protein structure of the baib > baia system, it can be seen in the average RMSF value of each system (Table 1). These results identify that the baia system has a more rigid structure than the baib because it has low fluctuations. More specifically, the fluctuations occurred in the loops (216–221) and \(\alpha \)-helix (256–276) regions. Meanwhile, the residues on the active site (radius of 5 Å from inhibitor coordinate) did not fluctuate excessively. It indicates that the region has a good interaction between inhibitor-amino acids and amino acids-amino acids. Additionally, the results show that each system has a relatively compact structure with insignificant fluctuations (~2.21 nm) by calculating the radius of gyration (RoG) (Fig. 5). These results indicate that during the 40 ns simulation time, each system has a well-folded structure [44]. It is confirmed by the average structure over a simulation time of 40 ns (Fig. 5). The results show that both systems have well-folded structures over time compared with crystal structures. The conformational dynamics of each system show good requirements during the simulation process. This becomes a crucial parameter in the evaluation of other variables, such as binding free energy and hydrogen bonding, which will be discussed in the next discussion.

Fig. 5.
figure 5

Root-mean-square of fluctuation and the radius of gyration plotted along 40 ns of trajectories. The superposition showed by the average structure of each system: initial coordinate (crystal) and after simulation along 40 ns (baia and baib).

Full size image
Table 1. Average parameters of trajectory analysis during the simulation time
Full size table

3.6.2. Hydrogen bonding analysis. Hydrogen bonding (H-bond) plays a crucial role in inhibitor-protein interactions [45–47]. H-bond analysis was calculated at 40 ns trajectories during the simulation time for each system (Table 2). The analysis of the H-bond occupation plays a major role in the evaluation of inhibitor-protein interactions. The results show that the baia system has four H-bonds and the baib system has only one H-bond is recorded during the simulation time. Unfortunately, the H-bond category shows a weak category because it has a very small H-bond occupancy, which is indicated by the percentage of fraction < 70% [42]. It is because during the simulation time the H-bond was recorded with the highest occupancy presentation only in the low-cost system (Fraction: 16.33%) of 3266 frames out of 20 000 total frames. However, the H-bond analysis provides a fairly clear picture in looking at the type of interaction between inhibitors and proteins.

Table 2. Hydrogen bonding analysis using 40 ns trajectories (cutoff value: distance < 3.5 Å and angle > 120°)
Full size table

3.6.3. Prediction of binding affinity. Determination of the binding affinity of each inhibitor to the target protein was calculated based on 4000 frames extracted from 40 ns trajectories. The binding affinity calculation process uses the MM–GBSA approach [48, 49]. The calculation process was performed on gas and solvation terms by considering several energy components (Table 3). The calculation process of each energy component becomes the main focus in looking at the contribution of the energy component to the free energy binding (\(\Delta {{G}_{{{\text{bind}}}}}\)) [50]. The energy component in the gas phase (\(\Delta {{G}_{{{\text{Gas}}}}}\)) shows that the van der Waals energy contribution (\({{E}_{{{\text{vdW}}}}}\)) has a good contribution compared to electrostatic energy (\({{E}_{{{\text{ele}}}}}\)). Besides, the energy contribution in the solvent phase (\(\Delta {{G}_{{{\text{sol}}}}}\)) shows an unfavorable contribution to the polar energy (\(\Delta G_{{{\text{solv}}}}^{{{\text{eve}}}}\)) for the Generalized Born model. However, overall the \(\Delta {{G}_{{{\text{bind}}}}}\) of each complex showed a good interaction. It is indicated by a higher negative value in thermodynamics terms. It is hoped that inhibitors that have a good binding affinity (higher negative value) can bind strongly to the protein active site. Each inhibitor that binds to two different active sites shows a relatively similar value (difference: ~ 1.06 kcal/mol). These results provide a fairly clear picture of the inhibitor’s ability to bind to the target protein.

Table 3. Determination of the energy component (kcal/mol) of each complex using the MM–GBSA approach. Data are shown as mean ± standard error of the mean (SEM)
Full size table

Integration Over Atomic Basins of Inhibitor-Protein Complexes

Integration over atomic basins of inhibitor-protein complexes was performed using quantum-chemical orbital-free AlteQ approach. This method has recently been shown to describe 3D electron density maps obtained for organic and inorganic compounds using high resolution low temperature X-ray diffraction data with very good quality, superior to other widely known quantum chemical methods [51, 52]. Then various integral molecular and atomic characteristics using AlteQ 3D maps and the Quantum Theory of Atoms in Molecules, or QTAIM, proposed by Bader [53] were computed. We used the software “Electron properties calculation: Integration over atomic basins” available on the site www.chemosophia.com. Radial integration over atomic basins was applied, starting from the atomic center and ending when the electron density surface reaches 0.001 a.u. or when the electron density is equal to the electron density of the neighboring atom. This method is described in more detail in [54]. The most important obtained integral characteristics including molecular solvent-accessible area \((S)\), volume \((V)\), charge \((Q)\), the number of electrons that atoms take from their neighbors (\({{N}_{{{\text{neigh}}}}}\)), the number of electrons which atoms donate to neighbors (\({{N}_{{{\text{on}}{\kern 1pt} {\text{neigh}}}}}\)) of the ligands in the gas phase and in the receptor–ligand complex as well as their differences \((\Delta )\) are given in Table 4.

Table 4. Integral characteristics (S, V, Q, \({{N}_{{{\text{neigh}}}}}\), \({{N}_{{{\text{on}}{\kern 1pt} {\text{neigh}}}}}\)) of the ligands in the gas phase and in the receptor–ligand complex as well as their differences \((\Delta )\)
Full size table

Integral characteristics show that baia insignificantly transfers electrons to the enzyme in the process of their interaction \((Q > 0)\), while baib insignificantly accepts electrons from the enzyme \((Q < 0)\). Moreover, the interactions of baib with the enzyme are more effective, which manifests itself in more significant overlaps of the electron clouds of baib with the enzyme: with a constant value of covalent overlaps, intermolecular overlaps are significantly enhanced, which is reflected in a serious increase in \({{N}_{{{\text{neigh}}}}}\) and \({{N}_{{{\text{on}}{\kern 1pt} {\text{neigh}}}}}\) numbers of electrons, as well as a serious reduction in volume and surface area.

In both cases, interactions are formed both with the transfer of electron density from the ligand to the enzyme, and vice versa. The atomic basins of the ligands in the complexes with the enzyme are shown in Fig. 6. Indeed, the basins of both oxygen atoms and hydrogen atoms of the ligand are truncated. The atoms of the enzyme in contact with them are shown schematically in Fig. 6.

Fig. 6.
figure 6

(Color online) Atomic basins of the ligands in the complexes with the enzyme: (a) baia and (b) baib.

Full size image

CONCLUSIONS

An alternative strategy that enables the quick identification of prospective therapeutic leads to treat quickly evolving viral infections is the repurposing of pharmaceutical products that have already received approval, including drug prospects. The HOMO and LUMO orbitals and MEP to study intermolecular interaction through the charge transfer within the molecule have been made obvious by DFT. The compounds baia and baib were discovered to be promising leads for the creation of drugs that are selective, safe, and effective against COVID-19, according to DFT studies. Integral characteristic investigations using the AlteQ method indicated that interactions were formed in baia and baib with the transfer of electron density from the ligand to the enzyme and vice versa. Overall, our results propose that baia and baib bearing good binding potency are components of Datura metel L that suggest the biologically safe profile of these compounds further supporting the potential of these compounds as starting points for therapeutics against COVID-19. Our studies were performed for identifying bioactive compounds that can inhibit COVID-19 Mpro effectively, as well as providing useful information for future studies. However, further studies should be conducted for the validation of these compounds using in vitro and in vivo models to pave a way for these compounds in drug discovery.