ORIGINAL RESEARCH article
An Analysis Based on Molecular Docking and Molecular Dynamics Simulation Study of Bromelain as Anti-SARS-CoV-2 Variants
- 1Department of Biology, Faculty of Mathematics and Natural Sciences, Sam Ratulangi University, Manado, Indonesia
- 2The University Centre of Excellence for Biotechnology and Conservation of Wallacea, Institute for Research and Community Services, Sam Ratulangi University, Manado, Indonesia
- 3Pharmacy Study Program, Faculty of Mathematics and Natural Sciences, Sam Ratulangi University, Manado, Indonesia
- 4Department of Animal Production, Faculty of Animal Husbandry, Sam Ratulangi University, Manado, Indonesia
- 5Department of Pharmacy, Faculty of Mathematics and Natural Sciences, Universitas Syiah Kuala, Banda Aceh, Indonesia
- 6Department of Biology, Faculty of Mathematics and Natural Sciences Education, Universitas Pendidikan Indonesia, Bandung, Indonesia
- 7Department of Pharmacy, BGC Trust University Bangladesh, Chittagong, Bangladesh
- 8Institute of Pharmacy, Martin-Luther University of Halle-Wittenberg, Halle, Germany
- 9Microbiology Laboratory, Department of Genetic Engineering and Biotechnology, University of Rajshahi, Rajshahi, Bangladesh
- 10Department of Pharmacology, College of Pharmacy, King Khalid University, Abha, Saudi Arabia
- 11Department of Clinical Laboratory Sciences, College of Applied Medical Sciences, Najran University, Najran, Saudi Arabia
- 12Department of Biotechnology and Genetic Engineering, University of Development Alternative, Dhaka, Bangladesh
- 13Department of Pharmaceutical Chemistry, Faculty of Pharmacy, Erciyes University, Kayseri, Turkey
The rapid spread of a novel coronavirus known as SARS-CoV-2 has compelled the entire world to seek ways to weaken this virus, prevent its spread and also eliminate it. However, no drug has been approved to treat COVID-19. Furthermore, the receptor-binding domain (RBD) on this viral spike protein, as well as several other important parts of this virus, have recently undergone mutations, resulting in new virus variants. While no treatment is currently available, a naturally derived molecule with known antiviral properties could be used as a potential treatment. Bromelain is an enzyme found in the fruit and stem of pineapples. This substance has been shown to have a broad antiviral activity. In this article, we analyse the ability of bromelain to counteract various variants of the SARS-CoV-2 by targeting bromelain binding on the side of this viral interaction with human angiotensin-converting enzyme 2 (hACE2) using molecular docking and molecular dynamics simulation approaches. We have succeeded in making three-dimensional configurations of various RBD variants using protein modelling. Bromelain exhibited good binding affinity toward various variants of RBDs and binds right at the binding site between RBDs and hACE2. This result is also presented in the modelling between Bromelain, RBD, and hACE2. The molecular dynamics (MD) simulations study revealed significant stability of the bromelain and RBD proteins separately up to 100 ns with an RMSD value of 2 Å. Furthermore, despite increases in RMSD and changes in Rog values of complexes, which are likely due to some destabilized interactions between bromelain and RBD proteins, two proteins in each complex remained bonded, and the site where the two proteins bind remained unchanged. This finding indicated that bromelain could have an inhibitory effect on different SARS-CoV-2 variants, paving the way for a new SARS-CoV-2 inhibitor drug. However, more in vitro and in vivo research on this potential mechanism of action is required.
Nature has an abundance of biological compounds that can be used as drug candidates to treat a wide range of diseases, both infectious and degenerative. Today, research is focused on discovering therapeutic agents that are derived from natural ingredients, such as plant materials. Pineapple, where the fruit is frequently consumed by the community, is one of the plants with nutraceutical properties. This plant has been used as a medicinal agent for a long time, and some of its benefits have been scientifically established. Bromelain, the key component of pineapple fruit and roots, is one of many active compounds found in pineapple that possesses therapeutic properties. Bromelain is made up of a variety of proteases, as well as phosphatase, glucosidase, peroxidase, cellulases, and glycoprotein (Bhattacharyya, 2008). Minor thiol endopeptidase, ananain, comosain, protease inhibitors, and organically bound calcium are all present in pineapple bromelain (Gautam et al., 2010; Bala et al., 2013).
Fresh pineapple is one of the favorite fruits in tropical countries. Apart from its nutritional value, bromelain’s active ingredients have been shown to have therapeutic properties, such as antiviral (e.g. anti-SARS-CoV-2), anti-inflammatory, anticancer, antimicrobial, and immunomodulatory function (Dopheide and Ward, 1981; Chandler and Mynott, 1998; Rathnavelu et al., 2016; Sagar et al., 2020; Chakraborty et al., 2021). Bromelain has a half-life of 6–9 h (Castell et al., 1997) and a plasma concentration of 2.5–4 ng/ml after an oral dose of 8.6 g per day (Kumakura et al., 1988). In an artificial stomach juice, the bromelain concentration of 3.66 mg/ml was detected after 4 h and it remained in artificial blood at a concentration of 2.44 mg/ml after 4 h (Pavan et al., 2012). According to recent research, bromelain can inhibit SARS-CoV-2 infection in VeroE6 cells by lowering the expression of the host cell receptor angiotensin-converting enzyme 2 (ACE-2) and the primed SARS-CoV-2 spike (S) protein. Bromelain has thus been suggested as an antiviral agent for COVID-19 therapy (Sagar et al., 2020). In another study, bromelain at a concentration of 100 g/ml was found to be effective in dissolving SARS-COV-2 spikes and envelope proteins. When combined with 20 mg/ml of acetylcysteine, the SARS-COV-2 spike and envelope proteins are fully disintegrated (Akhter et al., 2020). Multiple SARS-CoV-2 variants have recently been reported as circulating globally, according to the Center for Disease Control and Prevention (CDC). To date, at least four variants have been identified, the most notable of which are the B.1.1.7 lineage (the United Kingdom/United Kingdom variant/Alpha variant), B.1.351 lineage (South Africa/SA variant/Beta variant), P.1 lineage (Brazil/BR variant/Gamma variant), and B.1.429 lineage (California/US variant/Epsilon variant). The B.1.1.7, B.1.351, and P.1 lineages are considered variants of concern (VOC) because they have demonstrated a clear impact on disease transmission and severity, as well as immunity that may influence disease epidemiological situations. Meanwhile, the B.1.429 lineage is classified as a variant of interest because there is evidence that this variant possesses a mutation that alters the mode of transmission, the sensitivity of the test kit, the severity of symptoms, and the virus’s ability to evade the immune system. However, there is currently insufficient evidence, necessitating additional research.
The B.1.1.7 lineage has a mutation in the spike glycoprotein’s receptor-binding domain (RBD) at position 501, where the amino acid asparagine (N) has been replaced with tyrosine (Y) (abbreviated as N50Y). The B.1.351 lineage has multiple mutations in the spike protein which include K417N, E484K, and N501Y. The P.1 lineage has three RBD spike mutations: K417T, E484K, and N501Y. The B.1.429 lineage has the L452R mutation. Thus, these variants have mutations in the spike glycoprotein Xie et al. (2021), where this protein is responsible for the entry of the virus into the host cell through its attachment to the angiotensin-converting enzyme 2 (ACE2) receptor (Ni et al., 2020; Shang et al., 2020). A structural study analysis revealed that the RBD of the SARS-CoV-2 spike glycoprotein contains residues that are required for ACE2 binding (Lan et al., 2020). This receptor is a type I membrane protein (single transmembrane) protein found in several organs, including the oral and nasal mucosa, nasopharynx, lung, arteries, heart, kidneys, stomach, and intestines (Gheblawi et al., 2020; Bian & Li, 2021). As a result, several researchers have proposed spike protein as one of the targets for anti-SARS-CoV-2 drugs (Huang et al., 2020; Prajapat et al., 2020; Tallei et al., 2020; Khairan et al., 2021).
In the management of coronavirus disease 2019 (COVID-19), patients are treated with antiviral drugs (remdesivir, lopinavir, umifenovir, favipiravir, and oseltamivir), viral protease inhibitors (lopinavir and darunavir), anti-inflammatory agents (tocilizumab), and antimalarials (chloroquine and hydroxychloroquine) (Singh et al., 2020; Wu et al., 2020). However, there are still concerns about the efficacy of these drugs due to a lack of valid clinical trial data. As a result, natural compounds remain a promising alternative. Although bromelain has been shown to inhibit SARS-CoV-2 replication in vitro, with the development of new variants of this virus, computer modelling using molecular docking approaches and molecular dynamics simulation is required to evaluate bromelain’s ability to inhibit various variants of this virus. The purpose of this study was to investigate the interaction between bromelain and various mutated SARS-CoV-2 receptor-binding domains.
Materials and Methods
Preparation of Fruit Bromelain 3D Structure
The crystal structure of fruit bromelain (Bro) was unavailable in the PDB data-bank, thus homology modelling was utilized to generate its 3D structure using SWISS-MODEL (Schwede et al., 2003). The SWISS-MODEL structure assessment online server was used for validation (https://swissmodel.expasy.org/assess/) (Waterhouse et al., 2018). The protein sequence of fruit bromelain was obtained from GenBank with accession number QIM61761.1 (https://www.ncbi.nlm.nih.gov/protein/QIM61761.1).
The crystal structure of the SARS-CoV-2 spike receptor-binding domain (RBD) was retrieved from the Protein Data Bank with PDB ID 6YLA (https://www.rcsb.org/structure/6YLA), assigned as wild-type (WT). The 3D structure WT and its variants were created using homology modelling on the SWISS Model web-server and assigned as BR (P.1 lineage), SA (B.1.351 lineage), United Kingdom (B.1.1.7 lineage), and United States (B.1.427 lineage). Energy minimization was employed to overcome minor structural distortions, unfavorable interactions, and collisions introduced during the modelling phase.
Multiple Alignment Analysis
The multiple alignment analysis of the amino acid sequences of RBD WT and its variants was performed using the UCSF Chimera package release 1.15 (Pettersen et al., 2004).
Molecular Docking Analysis
The simulation of molecular docking was performed on the ClusPro protein-protein docking server (https://cluspro.bu.edu) (Kozakov et al., 2013; Kozakov et al., 2017; Vajda et al., 2017; Desta et al., 2020). The docking results were visualized in the LigPlot to confirm the binding position of bromelain and the receptors. As a comparison, other interaction graphics were generated using the EMBL-EBI tool PDBsum (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html). The binding affinity (ΔG) and dissociation constant (KD) predicted values were obtained from the Prodigy server (https://wenmr.science.uu.nl/prodigy/) (Vangone & Bonvin, 2015; Xue et al., 2016). The UCSF Chimera package release 1.15 was used to display the position of the interaction between RBD, hACE2, and bromelain.
Binding Free Energy Calculation of the Complexes
Estimations of the binding energies for complexes of Bro with RBD of wild type, B.1.427, B.1.1.7, B.1.351, and P.1 lineages were performed by using the MM-GBSA method. The method is a combination of molecular mechanics approaches with the generalized Born method and the surface area continuum solvation model.
Eint (internal) refers to bond, angle, and dihedral energies and Eele (electrostatic) refers to Coulomb force, while Evdw (van der Waals) refers to van der Waals interactions values. In the equation, GGB is the electrostatic solvation energy and, lastly, GSA is the non-electrostatic solvation energy. The MMGBSA. py (Hou et al., 2011; Miller et al., 2012) command which is implemented in Amber18 software was used to calculate energy values.
Molecular Dynamics (MD) Simulation
The dynamic behavior of Bro with WT, B.1.1.7, B.1.351, P.1, and B.1.427 lineage complexes was studied by MD simulation performed on Amber18 (Case et al., 2018) using the FF14SB protein force field. The solvation of the system was done by using the TIP3P water model with a margin of 10 Å (Jorgensen et al., 1983; Maier et al., 2015). Following that, the system was neutralized by the addition of counter ions. Thenceforward, two consecutive minimization steps were performed before the long MD simulation. In the first, 3,000 iterations (1,000 steepest descent and following 2,000 conjugate gradient) were submitted. In the second minimization step, 4,000 iterations (2,000 steepest descent and following 2,000 conjugate gradient) were performed. In the first minimization step, atom coordinates for the whole system were restrained to their initial coordinates with a force constant of 10 kcal/mol Å−2. In the second minimization step, the whole system was freely minimized to relieve atomic clashes/contacts in the entire system. The system was then heated at 300 K through 100 ps of MD with a time step of 2 fs per step while the whole system was restrained again with a force constant of 10 kcal mol-1 Å-2. That was followed by density evaluation of the system through 100 ps of MD with a time step of 2 fs per step. Before the long MD simulation, an equilibration step (a short MD for 200 ps) was carried out to equilibrate the system. Finally, a routine MD simulation for 100 ns was applied during which the temperature was kept at 300 K. The LEaP module of AMBER was used for protein-protein complex preparation, counter ion addition, solvation, and preparation of topology files. The PMEMD. CUDA GPU implementation was used for MD simulation production and the CPPTRAJ package of Amber18 was used for trajectory processing (Salomon-Ferrer et al., 2013).
Results and Discussion
The mutated variant of SARS-CoV-2 has spread worldwide and raised concerns about the vaccine’s effectiveness. The B.1.1.7 lineage is more difficult to neutralize than the parental virus, thus weakening the neutralization process by some members of a major class of public antibodies through light-chain interacting with the residue 501. This means that the virus can change the immune domains in a variety of ways to escape human immunity while still infecting and causing disease (Supasa et al., 2021).
The B.1.1.7 lineage was said to have started in the United Kingdom (Challen et al., 2021). The deletions at 69–70, 144, and substitutions K417N, K417T, E484K, N501Y, A570D, D614G, P681H, T716I, S982A, D1118H, and many others in the spike protein are predicted to affect SARS-CoV-2’s ability to transmit and infect. The N501Y is located on the ACE2 interacting surface (Supasa et al., 2021). This could have the consequence that the currently developed vaccine will not be effective against this B.1.1.7 lineage (Harrington et al., 2021). The variants BR and SA have also been reported to be more contagious. Many researchers around the world are thus using computational methods and techniques to identify new nature-based agents to solve this problem and reduce the fear of a pandemic.
Based on the information above, we considered natural resource-based management options for COVID-19. Bromelain, the active compound contained in pineapple fruit and stem, was discovered to have the ability to disintegrate spike protein in recent research (Sagar et al., 2020), and also alter the viral protein (Akhter et al., 2020). The interaction between Bro and RBDs is discussed in this article.
Fruit Bromelain 3D Structure
GenBank provides information about proteins that have been published in previous research. The platform provides information, among others, about catalytic activity residues. Our study utilizes this platform to obtain information on protein sequences from the fruit bromelain (Ananas comosus (L.) Merr.), which we then use to do homology modelling. The sequence of Bro extracted from pineapple fruit obtained from GenBank is shown in Figure 1.
FIGURE 1. The amino acid sequence of fruit bromelain of pineapple (Ananas comosus) retrieved from GenBank with accession number QIM61761.1.
The sequence has a cathepsin pro-peptide inhibitor domain (inhibitor 129 superfamily) at the N-terminus and a papain family cysteine protease (peptidase C1 superfamily) at the other end. The 3D structure of the fruit bromelain prepared using homology modelling on the SWISS-MODEL web-server is displayed in Figure 2A. A good model will typically have 90% of its residues in the Ramachandran plot’s allowable regions (Laskowski et al., 1993), and this information is in accordance with the model generated by SWISS-MODEL structure assessment (Figure 2B).
FIGURE 2. (A) The 3D structural model of bromelain fruit prepared using homology modelling (B) Ramachandran plot showing the residues in the most favoured regions of 90.55%.
Multiple Alignment Analysis
The spike glycoprotein comprises of the S1 subunit (14–685 residues), which is responsible for receptor binding, and the S2 subunit (686–1,273 residues), which is responsible for membrane fusion. The RBD is located at 319–541 residues (Huang et al., 2020). In this report, we used the S1 subunit of spike protein, where there is a receptor-binding domain (RBD), which is responsible for the attachment of the virus to hACE2. The RBD amino acid sequences of all variants analyzed in this study were aligned using chimera release 1.15 (Figure 3). The locations of all mutations were indicated. The WT has no mutations, the P.1 lineage (BR) has mutations at K417T, E484K, and N501Y, the B.1.351 lineage (SA) has mutations at K417N, E484K, and N501Y, the B.1.1.7 lineage (UK) has only one mutation at N501Y (Davies et al., 2021; Khan et al., 2021; Madhi et al., 2021; Wang et al., 2021), and the B.1.427 lineage (US) carries L452R mutation (Deng et al., 2021). Khan et al. (2021) further suggested that the N501Y-E484K variants hypothetically alter the binding affinity of RBD to hACE2, establish new interprotein contacts, and change the internal structural dynamics, resulting in increased binding and infectivity. The L452R mutation circumvents human leukocyte antigen (HLA)-restricted cellular immunity and enhances binding affinity for the viral receptor hACE2 (Motozono et al., 2021). The SA and BR mutations are significantly more lethal than the United Kingdom variant. However, it is reported that each of these mutants changes the binding affinity, establishes new inter-protein contacts, and modifies the internal structural dynamics, thereby increasing binding and eventually infectivity (Khan et al., 2021).
FIGURE 3. Multiple sequence alignment of the amino acid sequences of wild type (WT) SARS-CoV-2 spike receptor-binding domain (RBDs) and the four variants BR, SA, United Kingdom, and US.
Molecular Docking Analysis
Protein-protein docking is a molecular modelling technique that aims to predict the mutual orientation and position of two molecules forming a complex using computer algorithms and techniques. Simulation of molecular docking was performed on the Cluspro web-server (https://cluspro.bu.edu/login.php) to determine the possibility of interaction between RBDs and Bro. Docking is a method of interacting between molecules to determine the lowest energy created when a stable complex forms as a result of docking events. ClusPro ranks docking models based on the size of the conformation cluster to which they belong. ClusPro also provides two types of docking energies: the lowest energy among the conformations within a cluster of conformations and the core energy of the cluster (Xue et al., 2011).
The ClusPro score reflects the attempt to find the native site with the lowest free binding energy. Table 1 shows the size (number of members) of each cluster, the weighted energy score of the cluster center (i.e., the structure with the most neighboring structures in the cluster), and the energy score of the cluster’s lowest energy structure (Kozakov et al., 2017). As a result of mutations in RBD, the ClusPro score increased from −817.2 to −844.3 kJ/mol in Bro-UK variant complex and decreased to −744.3, −744.5, and −636.0 kJ/mol, in the Bro-BR, SA, and US variant complexes, respectively. It appears that Bro has a higher preference for the RBD United Kingdom variant compared to other RBDs. However, according to Kozakov et al. (2017), the ClusPro score should not be regarded as a measure of binding affinity. The native fold was typically the cluster with the largest number of low-energy structures (Comeau et al., 2004). The prodigy results for the ΔG calculation indicated that Bro binds to the BR and SA variants slightly stronger. The dissociation constant (KD) represents the equilibrium constant that exists when molecules bonded together in a complex dissociate. The smaller the dissociation constant value, the more tightly bonded a molecular complex is, or the greater the affinity between the molecules in the complex. Inferring from the binding affinity and dissociation constant values, Bro binds slightly more strongly to RBD variants than to WT.
TABLE 1. Weighted scores, binding affinity (ΔG), and dissociation constant (Kd) of the interaction of each RBD with bromelain.
The molecular docking complexes generated by ClusPro were used to analyze the amino acid residue interaction between Bro and RBDs. The interaction was generated by the EMBL-EBI tool PDBsum and LigPlot+ (Figure 4–8). Cannalire et al. (2020) reported that the RBD of SARS-CoV Tyr442, Leu472, Asn479, Thr487, and Tyr491 to bind to ACE2. On the other hand, several studies have shown that the key amino acids from the RBD of SARS-CoV-2, Leu455, Phe486, Gln493, and Asn501 form stronger interactions with the host receptor, together with Lys417, thus causing SARS-CoV-2’s affinity for ACE2 to be 20 folds compared to SARS-CoV (Lan et al., 2020; Yan et al., 2020). According to Koley et al. (2021), the interacting residues of RBD are Gln493, Tyr449, and Gly446. Khan et al. (2021), on the other hand, reported that hACE2-RBD forms hydrogen bonds between Glu30-Lys417, Glu35-Gln493, Glu38-Tyr449, Glu38-Gly496, Tyr41-Thr500, Tyr41-Thr500, Gln42-Gln498, Asn330-Thr500, Lys353-Gly502, Lys353-Gly496, and Lys353-Gln498, and also a salt bridge between amino acids Glu30 and Lys417.
FIGURE 4. Docking representation of the WT and bromelain complex. (A) the binding interface of the complex, (B) the binding interaction between the amino acids, and (C) interaction representation including hydrogen, salt bridges, and nonbonded interactions. Chain A is RBD WT and chain B is bromelain.
FIGURE 5. Docking representation of the BR variant and bromelain complex. (A) the binding interface of the complex, (B) the binding interaction between the amino acids, and (C) interaction representation including hydrogen, salt bridges, and nonbonded interactions. Chain A is RBD BR and chain B is bromelain.
FIGURE 6. Docking representation of the SA variant and bromelain complex. (A) the binding interface of the complex, (B) the binding interaction between the amino acids, and (C) interaction representation including hydrogen, salt bridges, and nonbonded interactions. Chain A is RBD SA and chain B is bromelain.
FIGURE 7. Docking representation of the United Kingdom variant and bromelain complex. (A) the binding interface of the complex, (B) the binding interaction between the amino acids, and (C) interaction representation including hydrogen, salt bridges, and nonbonded interactions. Chain A is RBD United Kingdom and chain B is bromelain.
FIGURE 8. Docking representation of the United States variant and bromelain complex. (A) the binding interface of the complex, (B) the binding interaction between the amino acids, and (C) interaction representation including hydrogen, salt bridges, and nonbonded interactions. Chain A is RBD United States and chain B is bromelain.
Table 2 shows the interacting residues between Bro and RBDs, where Bro occupies the site where hACE2 binds to RBD. Bro binds with Tyr449, Glu484, Gln493, Gly496, and Gln498 in RBD-WT, with Tyr449, Gly446, Glu484, Gln493, Gly496, Gln498, and Tyr501 in RBD-UK, with Tyr449, Gly446, Lys484, Gln493, Gly496, Gln498, and Tyr501 in RBD-BR, with Tyr449, Gly446, Lys484, Gln493, Gly496, Gln498, and Tyr501 in RBD-SA, and with Tyr449, Gly446, Arg452, Glu484, and Gly496 in RBD-US. This suggests that bromelain binds more to BRD variants than WT, except with the US variant. The substitution sites E484K (Glu484 → Lys484) and N501Y (Asn501 → Tyr501) in RBD variants also have bonds and nonbonded contacts with Bro. Bro has several nonbonded contacts with Arg452 in the US variant in the substitution site (L452R). Bro does not have interaction in the substitution site K417T (Lys417 → Thr417).
TABLE 2. The position of the interacting residues of bromelain with RBDs pocket and mutation sites (highlighted in bold). Key amino acid residues that play a role in binding RBD to hACE2 are marked with an italic font.
ACE2 and Bro can be considered competitive because they do not show a complete allosteric binding. It may be connected right next to the binding site of RBD and prevent the virus from binding with ACE2. Figure 9 illustrates the competitive binding model between bromelain and the ACE2 receptor at the RBD binding site. Several previous studies have shown that Bro is capable of removing the spike protein from a variety of viruses (Kennedy, 1974; Kharitonenkov et al., 1978). Bro has also been shown to be capable of inhibiting SARS-CoV-2 infection in VeroE6 cells via a disulfide bond-mediated mechanism (Sagar et al., 2020). Due to the paucity of research on the inhibition of SARS-CoV-2 RBD by plant-based proteins/enzymes, comparative data is difficult to obtain. However, considerable research has been conducted on the inhibitory effect of peptides on RBD (Barh et al., 2020; Salman et al., 2020; Schütz et al., 2020). By comparing to these studies, it appears promising that Bro could be used as an anti-SARS-CoV-2 agent.
FIGURE 9. The competitive binding model between bromelain (green) and the ACE2 receptor (yellow) on the binding site of RBD (light blue) of the spike glycoprotein.
The ability of Bro to recognize and interact with RBD is evidenced by the number of interface residues and hydrogen bonds as presented in Table 3. The hydrogen bond has an important role in the stabilization of macromolecular interactions by helping to stabilize the three-dimensional structures of the protein (Gurusaran et al., 2016). The number of interface residues in the RBD variants is lower than in WT. The number of H-bonds is higher in Bro and RBD variants, although the number of bonds does not differ much among RBDs. Table 4 shows that Bro and RBDs exhibit numerous H-bonds and nonbonded interactions, indicating that Bro interacts well with RBDs. Nonbonded interactions are typically divided into two types: electrostatic and van der Waals interactions in force field representation (Duan et al., 2020). The salt bridges were contributed by Glu484 (E484) in RBD WT, United Kingdom, and US variants, where the mutation occurs, with Lys162 of Bro. Salt bridges are frequently associated with structural driving forces that increase the stability of the interaction (Pylaeva et al., 2018). The interface area indicates that the two groups of complexes of Bro and RBDs overlap in terms of affinities and buried surface areas, with no clear boundary separating them (Chen et al., 2013).
TABLE 4. H-bonds between RBDs and Bromelain. Amino acids that are important in the interaction between RBD and hACE2 are marked with italics font, while bold font indicates the mutated amino acids.
Molecular Dynamics Simulation Analysis
A molecular dynamics simulation was conducted to analyze the binding stability of Bro and RBD complexes, where multiple descriptors were analyzed to understand the flexible and stable nature of the complexes. The last frame interaction between Bro and RBDs in MD simulation is shown in Table 5. Between the two molecules, there was an interacting shift in the amino acid. However, interactions with amino acids continue to exist at the active sites of the RBDs. In RBD WT, Bro still interacted with Tyr449, Glu484, Gln493, Gly496, although the interacting amino acids of Bro were slightly different. The interaction of Bro with the RBD UK variant was quite stable with amino acids Tyr449, Gly446, Glu484, Gln493, Gly496, Gln498, and Tyr501, although the amino acids in the interacting Bro underwent slight changes. There was a slight shift in the interaction between Bro and the RBD BR variant, where the interaction only occurred at Tyr449, Lys484, Gln498, Thr500, and Tyr501. The Bro no longer interacted with Gly446, Gln493, and Gly496. Bro had another interaction, namely with Thr500. In the RBD SA variant, Bro interacted with Tyr449, Gly446, Lys484, Gly496, Gln498, and Tyr501, and no longer with Gln493. Meanwhile, in the final frame, Bro interacts more with the RBD US variant, specifically with residues Tyr449, Gly446, Glu484, Arg452, Gln493, Gly496, and Gly498, which did not initially interact with Gln493 and Gln498. Given this fact, the last frames of the Bro-RBD complex were reliable.
TABLE 5. Interaction between bromelain and RBDs of the last frame MD simulation on amino acids that are important in the interaction between RBD and hACE2 (marked with italics font). The bold font indicates mutated amino acids.
Figure 10 illustrates that bromelain-RBD complexes, RBD, and bromelain were relatively stable across the simulation times where lower flexibility was observed. The lower deviation in RMSD was not found in the complexes of the BR variant where high RMSD was observed after 20ns simulation time. This complex maintained a similar profile across the simulation trajectories, but it showed a comparatively less stable profile due to a higher trend in RMSD. The Bro-SA variant complex was relatively stable and exhibited lower RMSD. The complex followed a similar pattern, but after 40 ns, it began to rise slightly due to flexibility. However, the complex remained stable after that and maintained a similar pattern until the end of the simulations. The complex of the Bro-UK variant had lower RMSD than other variants, which defines a more rigid profile of the complex. However, the complex Bro-US variant had a comparatively large difference in RMSD. The RMSD value of Bro-WT was lower than other complexes, which indicates the stable profile of the complex.
FIGURE 10. A molecular dynamics simulation. (A) Root mean square deviation of bromelain, RBDs, and the complex of bromelain and RBDs, and (B) Radius of gyration of the complex of bromelain and RBDs.
Moreover, the radius of gyration of the simulated complexes was analyzed to understand the mobile nature of the complexes and the complexes’ flexibility (Figure 10). The higher Rg correlates with the higher level of labile nature of the complexes, whereas the lower Rg profile indicates the contracted nature of the complexes. The Bro-WT, United States, United Kingdom, and SA variant complexes had relatively stable and lower Rg compared to the Bro-BR variant complex. They maintained a steady-state throughout the whole simulation time.
The Bro-BR variant complex had flexibility after 20 ns time, and higher than the other variants, which indicates the labile nature of this complex. The bromelain-WT RBD complex stayed stable with a 2 Å RMSD value up to 10 ns, then the RMSD of complex slowly rose between 10–15 ns and reached 2.5 Å, and the RMSD value, which continued to increase slowly, was 3 Å at 25 ns and stabilization was seen at this value until 50 ns. Between 50–70 ns, the RMSD value fluctuated at around 3 Å after 70 ns, then at 92 ns, the complex was uninformed and maintained at the level of 2.5 Å. At 92 ns, the RMSD rose sharply to 4 Å with some ripple, and the RMSD reached 4.5 Å, but from 98 ns, the RMSD stabilized at 3.5 Å. The Rg graph was observed to support this result. The average Rg value for the complex during the simulation was calculated as 27.25 Å.
The Bro-US complex stayed stable with a 3 Å RMSD value up to 27 ns. Subsequently, the RMSD of the complex slowly rose between 27–85 ns and reached 4 Å, and slowly decreased to about 3 Å end of simulations. The Rg value of the complex was stable between 27 and 28 Å throughout 100 ns with an average value of 27.75 Å.
The Bro-UK complex stayed stable until the 8 ns of the long MD, while the RMSD value gradually increased until 3.5 Å among the 8–30 ns interval. Following that, the complex’s RMSD remained around 2.5 Å until the 100 ns long MD. The Rg value of the complex showed an upward trend from 26.75 to 28.25 Å during the first 40ns or 4k frames. Afterward, in the early 40 ns, the Rg value decreased and stayed steady at around 27.5 Å until the end of 100 ns The average Rg value for the simulation was determined as 27.66 Å.
The Bro-SA complex was stable until 40 ns. After that, the RMSD value rose to 3 Å and even with some fluctuations, the RMSD value was stable at around 3 Å until the 100 ns The Rg graph of the complex showed an upward trend from 27.00 to 28.5 Å throughout the first 65 ns From 65 to 85 ns, some waves were observed in the Rg value, but after 85 ns until 100 ns, the value was stable, where the Rg value was observed at around 28.00 Å. The average value was observed as 27.76 Å.
The Bro-BR complex was stable until the 20 ns at below 2 Å, and then increased to about 4 Å up to 30 ns, and was stable with a small deviation between 35 and 100 ns The complex Rog is similar to the RMSD profile, and it reached 29 Å value until 30 ns and was stable at about 28.5 Å after 40 ns Lastly, the average Rg value of the complex was calculated as 28.16 Å.
Moreover, the root mean square fluctuations (RMSFs) of the simulated complexes were analyzed to understand the flexibility across the amino acid residues of the complexes Figure 11). The RMSF profile in their amino acid residues had lower RMSF than 2.5 Å, except for the BR variant on Asn221, SA variant on Asn107, Ala220, Asn221, Ser222, non on United Kingdom variant, United States variant on Asp31, Pro32, Asn107, Ser222, and WT variant on Asp31, Ala213, Ala219, Ala22. Moreover, the RBD had a similar lower RMSF profile for maximum residues, except for BR variant on Ile332, Thr333, Ser530, Thr531, Asn532, SA on variant Ile332, United Kingdom variant on Ile332, Thr333, Asn532, United States variant on Ile332, Thr333, Asn532, and WT on Ile332, Asn370, Asn532, which defines the stable and lower flexibility of the complex. However, there exceptions. In all RBD variants, regions 366–372, 445–448, 475–486, and 517–521 are relatively more flexible than the rest of the protein with RMSF values between 1.5–2.5Å.
Moreover, binding free energy (BFE) calculations were also performed for the 100 ns long MDs with 10,000 frames, with the 50 frames interval of each complex. Calculations were done as explained in the Methods part. Respectively, the complexes binding free-energy (∆G) were calculated as −125.89 kcal/mol (WT), −74.24 (BR variant), −80.80 (SA variant), −98.47 (United Kingdom variant) and −106.6 (US variant) (Table 6). Binding analyses revealed differences in the interactions of Bro and WT, BR variant, SA variant, United Kingdom variant, and US variant, with Bro’s binding affinity being stronger on WT than its variants. The complex Bro-BR variant has the least interaction.
The interaction between bromelain and various RBD variants of SARS-CoV2 was reported for the first time in this study. The detailed in silico analysis suggests that bromelain effectively binds to RBD’s active site, especially with the residues Gln493, Tyr449, Gly446, Gly496, and Gln498. As a result, the interaction between RBD and hACE2 was hypothetically interfered with significantly, thus suggesting the potential use of bromelain to prevent viral entry into the host cells. The current study could lay the groundwork for bromelain to be used as an inhibitor against SARS-CoV-2, but more in vitro and in vivo testing is needed before making any final conclusion.
Data Availability Statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Concept and design: TT, F. and IC. Acquisition of data: TT, IC, AY, TY, SM, and DK. Analyses and interpretation of data: TT, RI, IC, TE, TY, WS, and SM. Drafting the manuscript: TT, F, IC, TY, and TE. Critical revision of the manuscript for important intellectual content: TT, F, RI, TE, WS, SM, TH, AA, SA, MR, RJ, MK, and IC, Technical or material support: AY, IC, TA, AA, SA, MR, RJ and TY.
The authors would like to show gratitude to the Ministry of Education, Culture, Research, and Technology of the Republic of Indonesia for providing a research fund for the fiscal year 2021 through the Basic Research Scheme (Contract number: 145/E4.1/AK.04.PT/2021; Sub-contract number: 1997/UN12.13/LT/2021).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
All MD simulations presented here were carried out using resources provided by the Institute of Pharmacy, Martin-Luther University of Halle-Wittenberg, 06120 Halle/Saale, Germany.
Akhter, J., Pillai, K., Badar, S., Mekkawy, A., Valle, S. J., and Morris, D. L. (2020). In vitro study of BromAc on SARS-CoV-2 spike and envelope protein shows synergy and disintegration at modest concentrations. bioRxiv. doi:10.1101/2020.09.07.286906
Barh, D., Tiwari, S., Silva Andrade, B., Giovanetti, M., Almeida Costa, E., Kumavath, R., et al. (2020). Potential chimeric peptides to block the SARS-CoV-2 spike receptor-binding domain. F1000Res 9, 576. doi:10.12688/f1000research.24074.1
Cannalire, R., Stefanelli, I., Cerchia, C., Beccari, A. R., Pelliccia, S., and Summa, V. (2020). SARS-CoV-2 entry inhibitors: Small molecules and peptides targeting virus or host cells. Int. J. Mol. Sci. 21 (16), 5707. doi:10.3390/ijms21165707
Castell, J. V., Friedrich, G., Kuhn, C. S., and Poppe, G. E. (1997). Intestinal absorption of undegraded proteins in men: Presence of bromelain in plasma after oral intake. Am. J. Physiol. 273 (1 Pt 1), G139–G146. doi:10.1152/ajpgi.1997.273.1.g139
Chakraborty, A. J., Mitra, S., Tallei, T. E., Tareq, A. M., Nainu, F., Cicia, D., et al. (2021). Bromelain a Potential Bioactive Compound: A Comprehensive Overview from a Pharmacological Perspective. Life (Basel) 11 (4), 317. doi:10.3390/life11040317
Challen, R., Brooks-Pollock, E., Read, J. M., Dyson, L., Tsaneva-Atanasova, K., and Danon, L. (2021). Risk of mortality in patients infected with SARS-CoV-2 variant of concern 202012/1: Matched cohort study. BMJ 372, n579. doi:10.1136/bmj.n579
Chandler, D. S., and Mynott, T. L. (1998). Bromelain protects piglets from diarrhoea caused by oral challenge with K88 positive enterotoxigenic Escherichia coli. Gut 43 (2), 196–202. doi:10.1136/gut.43.2.196
Chen, J., Sawyer, N., and Regan, L. (2013). Protein-protein interactions: General trends in the relationship between binding affinity and interfacial buried surface area. Protein Sci. 22 (4), 510–515. doi:10.1002/pro.2230
Davies, N. G., Abbott, S., Barnard, R. C., Jarvis, C. I., Kucharski, A. J., Munday, J. D., et al. (2021). Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England. Science 372 (6538), eabg3055. doi:10.1126/science.abg3055
Deng, X., Garcia-Knight, M. A., Khalid, M. M., Servellita, V., Wang, C., Morris, M. K., et al. (2021). Transmission, infectivity, and neutralization of a spike L452R SARS-CoV-2 variant. Cell 184, 3426–3437. doi:10.1016/j.cell.2021.04.025
Duan, G., Ji, C., and Zhang, J. Z. H. (2020). Developing an effective polarizable bond method for small molecules with application to optimized molecular docking. RSC Adv. 10 (26), 15530–15540. doi:10.1039/d0ra01483d
Gautam, S. S., Mishra, S. K., Dash, V., Goyal, A. K., and Rath, G. (2010). Comparative study of extraction, purification and estimation of bromelain from stem and fruit of pineapple plant. Thai J. Pharm. Sci. 34 (2), 67–76.
Gheblawi, M., Wang, K., Viveiros, A., Nguyen, Q., Zhong, J. C., Turner, A. J., et al. (2020). Angiotensin-converting enzyme 2: SARS-CoV-2 receptor and regulator of the renin-angiotensin system: Celebrating the 20th Anniversary of the Discovery of ACE2. Circ. Res. 126 (10), 1456–1474. doi:10.1161/CIRCRESAHA.120.317015
Gurusaran, M., Sivaranjan, P., Dinesh Kumar, K. S., Radha, P., Thulaa Tharshan, K. P. S., Satheesh, S. N., et al. (2016). Hydrogen Bonds Computing Server (HBCS): An online web server to compute hydrogen-bond interactions and their precision. J. Appl. Cryst. 49 (2), 642–645. doi:10.1107/S1600576716002041
Harrington, D., Kele, B., Pereira, S., Couto-Parada, X., Riddell, A., Forbes, S., et al. (2021). Confirmed Reinfection with Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Variant VOC-202012/01. Clin. Infect. Dis., ciab014. doi:10.1093/cid/ciab014
Hou, T., Wang, J., Li, Y., and Wang, W. (2011). Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. J. Chem. Inf. Model. 51 (1), 69–82. doi:10.1021/ci100275a
Huang, Y., Yang, C., Xu, X. F., Xu, W., and Liu, S. W. (2020). Structural and functional properties of SARS-CoV-2 spike protein: potential antivirus drug development for COVID-19. Acta Pharmacol. Sin 41 (9), 1141–1149. doi:10.1038/s41401-020-0485-4
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79 (2), 926–935. doi:10.1063/1.445869
Khairan, K., Idroes, R., Tallei, T. E., Nasim, M. J., and Jacob, C. (2021). Bioactive Compounds from Medicinal Plants and their Possible Effect as Therapeutic Agents against COVID-19: A Review. Cnf 17 (1), 621–633. doi:10.2174/1573401317999210112201439
Khan, A., Zia, T., Suleman, M., Khan, T., Ali, S. S., Abbasi, A. A., et al. (2021). Higher infectivity of the SARS‐CoV‐2 new variants is associated with K417N/T, E484K, and N501Y mutants: An insight from structural data. J. Cel Physiol. 2021, 1–14. doi:10.1002/jcp.30367
Kharitonenkov, I. G., Siniakov, M. S., Grigoriev, V. B., Arefiev, I. M., Eskov, A. P., and Klimontovich, A. V. (1978). The length of the influenza virus spikes measured by photon correlation spectroscopy. FEBS Lett. 96 (1), 120–124. doi:10.1016/0014-5793(78)81075-8
Koley, T., Madaan, S., Chowdhury, S. R., Kumar, M., Kaur, P., Singh, T. P., et al. (2021). Structural analysis of COVID-19 spike protein in recognizing the ACE2 receptor of different mammalian species and its susceptibility to viral infection. 3 Biotech. 11 (2), 109. doi:10.1007/s13205-020-02599-2
Lan, J., Ge, J., Yu, J., Shan, S., Zhou, H., Fan, S., et al. (2020). Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor. Nature 581 (7807), 215–220. doi:10.1038/S41586-020-2180-5
Laskowski, R. A., MacArthur, M. W., Moss, D. S., and Thornton, J. M. (1993). PROCHECK: a program to check the stereochemical Quality of protein structures. J. Appl. Cryst. 26 (2), 283–291. doi:10.1107/s0021889892009944
Madhi, S. A., Baillie, V., Cutland, C. L., Voysey, M., Koen, A. L., Fairlie, L., et al. (2021). Efficacy of the ChAdOx1 nCoV-19 Covid-19 vaccine against the B.1.351 variant. N. Engl. J. Med. 384, 1885–1898. doi:10.1056/nejmoa2102214
Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theor. Comput 11 (8), 3696–3713. doi:10.1021/acs.jctc.5b00255
Miller, B. R., McGee, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012). MMPBSA.py: An efficient program for end-state free energy calculations. J. Chem. Theor. Comput 8 (9), 3314–3321. doi:10.1021/ct300418h
Motozono, C., Toyoda, M., Zahradnik, J., Ikeda, T., Seng Tan, T., Ngare, I., et al. (2021). An emerging SARS-CoV-2 mutant evading cellular immunity and increasing viral infectivity. BioRxiv. doi:10.1101/2021.04.02.438288
Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M., Meng, E. C., et al. (2004). UCSF Chimera--a visualization system for exploratory research and analysis. J. Comput. Chem. 25 (13), 1605–1612. doi:10.1002/jcc.20084
Rathnavelu, V., Alitheen, N. B., Sohila, S., Kanagesan, S., and Ramesh, R. (2016). Potential role of bromelain in clinical and therapeutic applications. Biomed. Rep. 5 (3), 283–288. doi:10.3892/br.2016.720
Sagar, S., Rathinavel, A. K., Lutz, W. E., Struble, L. R., Khurana, S., Schnaubelt, A. T., et al. (2020). Bromelain inhibits SARS-CoV-2 infection in VeroE6 cells. bioRxiv. doi:10.1101/2020.09.16.297366
Salman, S., Shah, F. H., Chaudhry, M., Tariq, M., Akbar, M. Y., and Adnan, M. (2020). In silico analysis of protein/peptide-based inhalers against SARS-CoV-2. Future Virol. 15 (9), 557–564. doi:10.2217/fvl-2020-0119
Salomon-Ferrer, R., Götz, A. W., Poole, D., le Grand, S., and Walker, R. C. (2013). Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh ewald. J. Chem. Theor. Comput 9 (9), 3878–3888. doi:10.1021/ct400314y
Schütz, D., Ruiz-Blanco, Y. B., Münch, J., Kirchhoff, F., Sanchez-Garcia, E., and Müller, J. A. (2020). Peptide and peptide-based inhibitors of SARS-CoV-2 entry. Adv. Drug Deliv. Rev. 167, 47–65. doi:10.1016/j.addr.2020.11.007
Singh, A. K., Singh, A., Shaikh, A., Singh, R., and Misra, A. (2020). Chloroquine and hydroxychloroquine in the treatment of COVID-19 with or without diabetes: A systematic search and a narrative review with a special reference to India and other developing countries. Diabetes Metab. Syndr. 14 (3), 241–246. doi:10.1016/j.dsx.2020.03.011
Supasa, P., Zhou, D., Dejnirattisai, W., Liu, C., Mentzer, A. J., Ginn, H. M., et al. (2021). Reduced neutralization of SARS-CoV-2 B.1.1.7 variant by convalescent and vaccine sera. Cell 184 (8), 2201–e7. e7. doi:10.1016/j.cell.2021.02.033
Tallei, T. E., Tumilaar, S. G., Niode, N. J., Fatimawali, B. J., Kepel, B. J., Idroes, R., et al. (2020). Potential of Plant Bioactive Compounds as SARS-CoV-2 Main Protease (Mpro) and Spike (S) Glycoprotein Inhibitors: A Molecular Docking Study. Scientifica 2020, 1–18. doi:10.1155/2020/6307457
Wang, P., Casner, R. G., Nair, M. S., Wang, M., Yu, J., Cerutti, G., et al. (2021). Increased resistance of SARS-CoV-2 variant P.1 to antibody neutralization. Cell Host & Microbe 29, 747–751. doi:10.1016/j.chom.2021.04.007
Waterhouse, A., Bertoni, M., Bienert, S., Studer, G., Tauriello, G., Gumienny, R., et al. (2018). SWISS-MODEL: Homology modelling of protein structures and complexes. Nucleic Acids Res. 46 (W1), W296–W303. doi:10.1093/nar/gky427
Wu, R., Wang, L., Kuo, H. D., Shannar, A., Peter, R., Chou, P. J., et al. (2020). An Update on Current Therapeutic Drugs Treating COVID-19. Curr. Pharmacol. Rep. 6 (3), 1–15. doi:10.1007/s40495-020-00216-7
Xie, X., Liu, Y., Liu, J., Zhang, X., Zou, J., Fontes-Garfias, C. R., et al. (2021). Neutralization of SARS-CoV-2 spike 69/70 deletion, E484K and N501Y variants by BNT162b2 vaccine-elicited sera. Nat. Med. 27 (4), 620–621. doi:10.1038/s41591-021-01270-4
Xue, L. C., Jordan, R. A., El-Manzalawy, Y., Dobbs, D., and Honavar, V. (2011). Ranking Docked Models of Protein-Protein Complexes Using Predicted Partner-specific Protein-Protein Interfaces: A Preliminary Study. ACM BCB 2011, 441–445. doi:10.1145/2147805.2147866
Xue, L. C., Rodrigues, J. P., Kastritis, P. L., Bonvin, A. M., and Vangone, A. (2016). PRODIGY: A web server for predicting the binding affinity of protein-protein complexes. Bioinformatics 32 (23), 3676–3678. doi:10.1093/bioinformatics/btw514
Keywords: bromelain, receptor-binding domain, SARS-CoV-2, mutation, variants
Citation: Tallei TE, Fatimawali , Yelnetty A, Idroes R, Kusumawaty D, Emran TB, Yesiloglu TZ, Sippl W, Mahmud S, Alqahtani T, Alqahtani AM, Asiri S, Rahmatullah M, Jahan R, Khan MA and Celik I (2021) An Analysis Based on Molecular Docking and Molecular Dynamics Simulation Study of Bromelain as Anti-SARS-CoV-2 Variants. Front. Pharmacol. 12:717757. doi: 10.3389/fphar.2021.717757
Received: 31 May 2021; Accepted: 09 August 2021;
Published: 20 August 2021.
Edited by:Luca Rastrelli, University of Salerno, Italy
Copyright © 2021 Tallei, Fatimawali, Yelnetty, Idroes, Kusumawaty, Emran, Yesiloglu, Sippl, Mahmud, Alqahtani, Alqahtani, Asiri, Rahmatullah, Jahan, Khan and Celik. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.