Introduction

SARS-CoV-2 is an enveloped; single stranded positive sense RNA virus which affects the respiratory tract of warm-blooded animals, including humans [1]. SARS-CoV-2 is associated with the common cold, pneumonia and severe acute respiratory syndrome (SARS). Coronavirus was first isolated in 1937 from an infectious bronchitis virus in chickens, which decimated the poultry population. These viruses cause 15–30% of all common colds (https://www.medicalnewstoday.com/articles/256521). Human coronavirus was discovered in the year 1960 in the nose of a patient suffering from common cold, and OC43 and 229E were the two human coronaviruses responsible for common cold [2]. A novel coronavirus, formally named as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has caused coronavirus infection 2019 (COVID-19) around the world. According to World Health Organization (WHO), SARS-CoV-2 has caused a beginning outbreak in Wuhan, city of China [3••]. China was initially unaware of the direct transmission of the 2019 novel coronavirus (2019-nCoV). However, patients without symptoms were recognized as the major source of infection, and human-to-human transmission was also confirmed [4]. Four structural proteins have been identified in SARS-CoV-2, i.e. envelope, membrane, nucleocapsid and spike proteins [5]. The envelope and membrane proteins play an important role in the assembly of coronavirus. In contrast, nucleocapsid protein is needed for viral RNA synthesis [6]. In SARS-Cov-2, virus binding and its entry are governed by critical spike glycoprotein having subunits S1 with 685 amino acids and S2 with 588 amino acids [5, 7]. S1 subunit has receptor binding domain (RBD) that facilitates the entry of the virus into the host cells via ACE2 (angiotensin-converting enzyme 2) receptor [8]. External sub-domain of RBD has highest variation of amino acids, which accounts for direct interaction between the host ACE2 receptor and the spike protein of the virus [5, 8]. S2 subunit remains conserved among the SARS-CoV-2 viruses and shows 99% similarity with the bat SARS-CoVs [5]. The receptor binding domain (RBD) of SARS CoV-2 spike protein binds more strongly to human ACE2 receptors than SARS-CoV and can act as main target to prevent the viral attachment to ACE2-expressing cells [9]. The SARS-CoV-2 pandemic has taken a huge toll across the world, while therapeutics are still being discovered or developed. The primary target for most of the viral diseases including SARS-CoV-2 is the upper or lower respiratory tracts. Apart from vaccines, there are many natural phytocompounds that have a protective effect on the respiratory system. To build up the drug discovery pipeline for the treatment of SARS-CoV-2 that majorly affects the respiratory system, scientists across the world are taking new approaches such as repurposing the existing antiviral and anti-malarial drugs and in silico approach to explore and identify new molecules or phytocompounds that could be further tested under in vitro and in vivo systems. The in silico approach could help in building the repository and pipeline for clinical trials and reduce the time for discovering new drug candidates.

One such important class of compounds is methylxanthines present in coffee, tea and cacao [10], The most famous beverage and consumed by two-thirds of the world’s population. Caffeine/thiene, theobromine and theophylline are the most important methylxanthines.

Drinking beverages such as tea and coffee has been considered as health-promoting benefits since old-fashioned times. Caffeine-related methylxanthines are products of caffeine metabolism in vivo (paraxanthine, theobromine and theophylline). Theophylline is in clinical use for the treatment of asthma and is a well-known bronchodilator and reduces the symptoms of asthma, wheezing, coughing and breathlessness. Theophylline has also been reported as anti-inflammatory. The inhibition of phosphodiesterase (PDE) 3 is thought to cause bronchodilation, but the anti-inflammatory activity may be due to inhibition of PDE4 and activation of histone deacetylase-2, which switches off the activated inflammatory genes [11]. Theophylline reverses corticosteroid resistance via this process, which may be essential in extreme asthma and COPD (chronic inflammatory lung disease), where histone deacetylase-2 activity is reduced (Barnes, 2013). Caffeine (1,3,7-trimethylxanthine) is a methyl-xanthine alkaloid that is gaining popularity due to its wide range of physiological and pharmacological effects [12, 13]. Caffeine has been shown to have antiviral properties against HSV-1 (1,2). Yamazaki and Tagaya (1980) [14] published tentative findings of caffeine’s antiviral activity in viruses such as poliovirus, influenza virus, HSV-1 and vaccinia virus. Olson and Consigli (1979) [15] reported that caffeine has antiviral properties against NDV (Newcastle disease virus) by inhibiting the synthesis of viral RNA and proteins in infected cells. According to current consensus [12, 14,15,16,17], caffeine tends to prevent the development of NDV, human immunodeficiency virus, polyomavirus, HSV-1 and vaccinia virus. Caffeine and other methylxanthines have recently been shown to inhibit the replication of infectious HIV-1 strains, with the integration stage of the HIV-1 life cycle being the caffeine’s target [18]. Keeping in mind the importance of methylxanthines in inhibiting certain classes of viruses and health-promoting benefits, the current study was designed to test the interaction of methylxanthines (caffeine/thiene, methylxanthine, theobromine, theophylline, xanthine) with three different target proteins (spike, main protease and RNA protein) of SARS-CoV-2 using in silico tools (Fig. 1).

Fig. 1
figure 1

Hypothetical binding mechanisms of methylxanthines to SARS-CoV-2. Methylxanthines binds spike protein (6LZG), nucleocapsid (6M3M) and main protease (6LU7). Based on our in silico studies, it is hypothesized that methylxanthine competes with spike protein of SARS-CoV-2 for binding to ACE2 receptor. The proposed interactions between methylxanthines and ACE2 receptor need to be experimentally verified in the future. Normal entry of SARS-CoV-2 through ACE2 receptor is shown in the boxed part of the diagram

Full size image

Methods

Bioinformatic Tools

Open Babel GUI [19], UCSF Chimera 1.8.1 [20]. Pubchem (www.pubchem.com), RCSB PDB (http://www.rscb.org/pdb), PDBsum (www.ebi.ac.uk/pdbsum), SwissADME and Autodock vina were used in the present investigation.

Ligand Preparation

Five major methylxanthines such as caffeine/thine, methylxanthine, theobromine, theophylline and xanthine (which are either found in coffee or the degradation products of caffeine) were selected in the current study. The 3-dimensional structures of all the methylxanthines, chloroquine (standard antimalarial drug) and lopinavir (standard antiretroviral drug) were retrieved from the pubchem (www.pubchem.com) in.sdf format. Open Babel was used to convert the phytocompounds.sdf file into PDB format [21]. Table 1 shows the molecular structure, compound CID and molecular weight of selected phytocompounds chloroquine (standard antimalarial drug) and lopinavir (standard antiretroviral drug).

Table 1 Molecular structure, compound CID and molecular weight of selected compounds and selected drug as standard drug
Full size table

Retrieval of SARS CoV-2 Proteins

The crystal structures of three SARS CoV-2 proteins were downloaded in pdb format from the protein databank, namely nucleocapsid protein N-terminal RNA binding domain (PDB ID: 6M3M) [22], main protease (6LU7) [23] and spike receptor-binding domain (PDB ID: 6LZG) [24]. For virtual screening, xanthine derivatives such as caffeine/thiene, methylxanthine, theobromine, theophylline and xanthine were used.

Target Protein Preparation

The high-resolution structure of SARS-CoV2 spike receptor-binding domain (PDB ID: 6LZG as shown in Fig. 2A,B), main protease (PDB ID: 6LU7 shown in Fig. 2C,D), nucleocapsid protein N-terminal RNA binding domain (PDB ID: 6M3M shown in Fig. 2E,F) and the Auto-Dock tool 1.5.4 was used for the preparation of the protein structures. The selected protein structures were processed for creating disulfide bonds, assigning proper bond orders and the addition of missing hydrogens in the raw structure. Hydrogen bond optimization was assigned using non-hydrogen atoms of protein.

Fig. 2
figure 2

Structures of selected SARS-CoV-2 S proteins: spike protein 6LZG (A: surface view representation and B: secondary structure representation); main protease 6LU7 (C: surface view representation and D: secondary structure representation); nucleocapsid protein N-terminal RNA binding domain (E: surface view representation and F: secondary structure representation). Helices are shown in red, -sheets in yellow and loops in green color (A–F)

Full size image

ADMET and Toxicity Prediction of Phytocompounds

Absorption, Distribution, Metabolism, Excretion and Toxicity (ADMET) screening was done to determine the absorption, toxicity and drug-likeness properties of selected methylxanthines (ligands). The 3D structures of ligands (caffeine, methylxanthine, theobromine, theophylline, thiene and xanthine) were saved in.smiles format and were uploaded on SWISSADME (Molecular Modeling Group of the SIB (Swiss Institute of Bioinformatics), Lausanne, Switzerland), admetSAR (Laboratory of Molecular Modeling and Design, Shanghai, China) and PROTOX webservers (Charite University of Medicine, Institute for Physiology, Structural Bioinformatics Group, Berlin, Germany) for ADMET screening. SWISSADME is a web tool used for the prediction of ADME and pharmacokinetic properties of a molecule. The predicted outcome includes lipophilicity, water solubility, physicochemical properties, pharmacokinetics, drug resemblance, medicinal chemistry and the Brain or Intestinal Estimated Permeation Method (blood–brain barrier and PGP prediction) [25]. admetSAR can predict about fifty ADMET properties and provides ADMET profiles for the query molecules. The following are the toxicity indicators: classes: I Compounds with LD50 values less than 50 mg kg−1, (ii) Compounds with LD50 values greater than 50 mg kg−1 but less than 500 mg kg−1 (Cheng et al., 2012; Yang et al., 2019). PROTOX is a server that predicts the LD50 value and toxicity class of a question molecule in rodents. The following are the toxicity classes: (Class 1: fatal if swallowed (LD50 5) and Class 2: fatal if swallowed (55,000) [26, 27••, 28••].

Molecular Docking

All the xanthine derivatives were used to perform molecular docking studies against three different proteins of SARS-CoV-2 to find the potential hits for further drug discovery experiments. In the current study, AutoDock vina was used for performing the docking studies [29, 30]. AutoDock vina uses a Lamarckian Genetic Algorithm (LGA) and is based on a semi empirical free energy force field. In all the proteins, the grid was manually adjusted within the key active binding amino acid residues. The high scoring compounds were chosen based on the binding energy of the ligands with the receptor. The lower ligand binding energy with the receptor is an indication of high affinity for the target receptor.

The.pdb complex of protein and ligands was analyzed by discovery studio (www.ebi.ac.uk/pdbsum) to classify the protein–ligand interactions. Chimera performed detailed visualization and comparison of docked sites of target proteins and ligands.

Molecular Dynamics Simulation

MD simulation was done to study the stability of protein ligand complexes over the period of 100 ns. For the authenticity of the data, MD simulation was done with two different software such as GROMACS 2018.3 (simulation in duplicate) and Desmond program version 2.0 academic version (simulation in triplicate). Among all the selected compounds, theophylline showed the best interactions with three selected targets (6LZG1, 6LU7 and 6M3M) of SARS-CoV-2.

MD Simulation by Using GROMACS 2018.3 Software

Molecular dynamics study of the theophylline complexes with 6LZG1, 6LU7 and 6M3M was done using GROMACS 2018.3 in duplicate [31] software which was installed in ubuntu 18.04 LTS. Docked structures of the protein–ligand complexes (theophylline with 6LZG1, 6LU7 and 6M3M) were used in the simulation study. Proteins were processed, and the topology file was prepared by using pdb2gmx and GROMOS54a7_atb.ff (imported from the automated topology builder website) force field, while the ligand topology file was prepared by using the Automated Topology Builder (ATB) version 3.0. The solvent addition was done in a cubic box by using a box distance of 1.0 nm from closest atom in the protein. To neutralize the device, the Cl− ions calculated from the genion module for each protein were used. The energy was minimized by using the steepest descent algorithm with 50,000 steps and a cumulative force of 5 kJ mol−1, as well as the Verlet cutoff scheme with particle mesh Ewald (PME) columbic interactions. During the equilibration process, position restraints were used. Following that, NVT equilibration was performed at 300 K with 100 ps steps, and NPT equilibration was performed with Parrinello-Rahman (pressure coupling), 1 bar reference pressure and 100 ps steps. The LINCS algorithm was used to constrain the length of all bonds. For long-range electrostatics, the particle mesh Ewald (PME) algorithm was used. The protein–ligand complex’s MD was run for 100 ns (in duplicate). Following efficient completion of a 100-ns molecular dynamic simulation, the root means square deviation (RMSD) of backbone residues, the number of hydrogen bonds and the root mean square fluctuations (RMSF), radius of gyration (RG) and solvent accessible surface area (SASA) were calculated.

Estimation of Free Energy of Binding

The free energy of binding calculations was done using the standalone program, g mmpbsa, following the protocol described by Kumari et al. (2014) [32]. Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) is a widely used method to calculate the free binding energy and to predict the stability of the complex. Modified molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) is an open-source software and employed to calculate the free energy of binding between two defined groups. Recently, MM-PBSA algorithm has been utilized as a scoring function in computational drug design [33, 34]. The free energy of binding is based on the theory:

$${G}_{\mathrm{binding}}={G}_{\mathrm{complex}}-({G}_{\mathrm{protein}}+{G}_{\mathrm{ligand}})$$

where Gcomplex is the total free energy of the protein-inhibitor complex, and Gprotein and Gligand are total free energies of the separated form of protein and inhibitor in solvent, respectively [33]. The average binding energy calculations were done by a python script provided in g mmpbsa stand-alone program.

MD Simulation by Using Desmond Program Version 2.0 (Academic Version)

To check the accuracy of docking observations and the stability of protein–ligand interactions, dynamic simulations (MD) simulations were performed in triplicate using the Desmond software version 2.0 (academic version) [28••, 30, 34,35,36,37,38]. All the three complexes were solvated in TIP3P water model with cubic periodic box containing Simple Point Charge (SPC) (101,010) with Optimized Potentials for Liquid Simulations (OPLS) all-atom force field 2005. MD simulation was performed to check the accuracy of docking observations and the stability of protein–ligand interactions. To run the MD simulations, the OPLS 2005 force field parameters were combined with periodic boundary conditions in the NPT ensemble method, which included a relaxation time of 1 ps at a constant temperature of 300 K, a constant length, a smooth particle mesh Ewald (PME) (with a 10–9 tolerance limit) and a cutoff distance of 9.0 Å. 100-ns production time had been chosen for the analysis of protein structure evaluation by every 1 ns. To decipher the stability, an average structure was chosen from the MD simulation corresponding to production phase. In addition, root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bond, radius of gyration (Rg), histogram for torsional bonds were performed for the investigation of structural changes with dynamic role of the protein–ligand complexes for 100 ns. After successful completion of the MD simulation for 100 ns, MM-GBSA calculation was done in triplicates.

Results

ADMET Prediction and Toxicity Analysis of Selected Phytocompounds

The ADMET properties predicted by SwissADME of selected phytocompounds and drugs Chloroquine and Lopinavir are summarized in Table 2. A Log Po/w value of 5 implies strong aqueous solubility, which means that a sufficient amount of drug will enter and stay within the body after oral administration. The Log Po/w value of all the selected phytocompounds is 5 (Table 2). The capacity of a compound to permeate into the cells is indicated by TPSA (total polar surface area). For good permeation of the compound through the cell membrane, a TPSA value of 94.40 Å is needed, and a value of 90 Å is required to pass through the blood–brain barrier. Except for xanthine, all other phytocompounds (caffeine/thiene, methylxanthine, theobromine and theophylline) showed a TPSA value of 90 Å, indicating good permeability through blood–brain barrier. Lipinski’s rule of five aids in determining a compound’s drug likeness; an orally active drug does not break more than one of Lipinski’s rules. None of the selected ligands violated any of Lipinski’s rule of five (Table 2). Table 3 summarizes the toxicity expected by PROTOXII and admetSAR. Table 3 shows that all the selected phytocompounds are non-carcinogenic and non-cytotoxic in nature and are thus safe to use. Caffeine/thiene, methylxanthine, theobromine, theophylline and xanthine had lower LD50 values than all other phytocompounds as measured using Protox II.

Table 2 ADMET properties of selected phytocompounds and standard drugs (Chloroquine and Lopinavir) predicted by SwissADME
Full size table
Table 3 Toxicity prediction of selected phytocompounds and standard drugs (Chloroquine and Lopinavir) predicted by admetSAR and PROTOX‑II software
Full size table

Docking Analysis

The selected phytocompounds showed relatively good binding affinity with all the selected target proteins of SARS CoV-2 (spike receptor-binding domain PDB ID: 6LZG, nucleocapsid protein N-terminal RNA binding domain PDB ID: 6M3M and main protease PDB ID: 6LU7) of SARS-CoV-2. Theophylline showed the best interaction binding affinity with 6LZG (− 5.7 kcal mol−1), followed by methylxanthine (-5.3 kcal mol−1), xanthine (− 5.2 kcal mol−1), theobromine (− 4.9 kcal mol−1) and caffeine (− 4.8 kcal mol−1) as shown in Figs. 3, 4 and Table 4. Interestingly, nucleocapsid protein N-terminal RNA binding domain (6M3M) also showed best interaction with theophylline (− 5.8 kcal mol−1), followed by theobromine (− 5.5 kcal mol−1), xanthine (− 4.9 kcal mol−1), methylxanthine (− 4.7 kcal mol−1) and caffeine (− 4.6 kcal mol−1) as shown in Figs. 4, 5, 6 and Table 4. Similarly, the main protease (6LU7) also showed best binding interaction with theophylline (kcal mol−1 − 6.5), followed by theobromine (− 5.8 kcal mol−1), caffeine (− 5.7 kcal mol−1), methylxanthine and xanthine (− 5.5 kcal mol−1) as shown in Figs. 7, 8 and Table 4. Binding affinity of Lopinavir was much stronger (6LZG: − 7.1 kcal mol−1, 6M3M: − 6.4 kcal mol−1, 6LU7: − 8.2 kcal mol−1) as compared to selected phytocompounds (Table 4) and chloroquine also represented in Table 4. In contrast to Lopinavir, Chloroquine showed weak interaction with all the three targets (6LZG: − 5.1 kcal mol−1, 6M3M: − 4.9 kcal mol−1, 6LU7: − 5.7 kcal mol−1) but comparable to the selected compounds (Table 4). The results of this study suggested that theophylline showed the best interaction with all the three selected target protein of SARS-CoV-2. Complexes of theophylline with 6LZG, 6M3M and 6LU7 were further analyzed by MD simulation.

Fig. 3
figure 3

2D interactions of xanthine derivatives with spike receptor-binding domain (PDB ID: 6LZG): (A) caffeine, (B) methylxanthine, (C) theobromine, (D) theophyllin, (E) xanthine and (F) chloroquine. Nature of interactions is shown by different colors as indicated

Full size image
Fig. 4
figure 4

Electrostatic surface view representation of receptor-binding domain (PDB ID: 6LZG) with bound xanthine derivatives: (A) caffeine, (B) chloroquine, (C) methylxanthine, (D) theobromine, (E) theophylline and (F) xanthine. The color coding is by electrostatic potential; negative, positive and neutral regions are shown in blue, red and white respectively (image generated using Pymol)

Full size image
Table 4 Interacting amino acid residues and E- total of xanthine derivatives with three different target proteins (6LZG, 6M3M and 6LU7) of SARS CoV-2
Full size table
Fig. 5
figure 5

2D interaction of xanthine derivatives with nucleocapsid protein N-terminal RNA binding domain (PDB ID: 6M3M): (A) caffeine, (B) methylxanthine, (C) theobromine, (D) theophylline, (E) xanthine and (F) chloroquine. Nature of interactions is shown by different colors as indicated

Full size image
Fig. 6
figure 6

Electrostatic surface view representation of nucleocapsid protein N-terminal RNA binding domain (PDB ID: 6M3M) with bound xanthine derivatives: (A) caffeine, (B) chloroquine, (C) methylxanthine, (D) theobromine, (E) theophylline and (F) xanthine. The color coding is by electrostatic potential; negative, positive and neutral regions are shown in blue; red and white respectively (image generated using Pymol)

Full size image
Fig. 7
figure 7

2D interactions of xanthine derivatives with SARS-CoV-2 main protease (PDB ID: 6LU7): (A) caffeine, (B) methylxanthine, (C) theobromine, (D) theophylline, (E) xanthine and (F) chloroquine. Nature of interactions is shown by different colors as indicated

Full size image
Fig. 8
figure 8

Electrostatic surface view representation of SARS-CoV-2 main protease (PDB ID: 6LU7) with bound xanthine derivatives: (A) caffeine, (B) chloroquine, (C) methylxanthine, (D) theobromine, (E) theophylline and (F) xanthine. The color coding is by electrostatic potential; negative, positive and neutral regions are shown in blue, red and white respectively (image generated using Pymol)

Full size image

Each amino acid residue is numbered followed by three letter codes. Amino acid residues in boldface letters are the phosphorylation sites on serine and threonine, predicted by Davidson et al. (2020).

MD Simulation by Using GROMACS 2018.3 Software

MD simulation of theophylline in complex with 6LZG, 6LU7, 6M3M was done for validation of molecules docking study and stability of theophylline in complex with SARAS-CoV-2 proteins (6LZG, 6LU7 and 6M3M). MD simulation was done for 100 ns.

RMSD Analysis of Protein–Ligand Complexes

RMSD of theophylline in complex with 6LZG was stabilized at 10ns/0.4 nm in both the runs. RMSD of theophylline in complex with 6LZG remains stable up to 80ns and further stabilized at 90–100ns at 0.7–0.8 nm in first run. No major fluctuation was observed in the second run (Fig. 9A), and minor fluctuations are however still within the acceptable limits. The RMSD of theophylline in complex with 6LU7 and 6M3M was stable throughout the entire simulation (Figure 9B,C) during both the runs. The estimated RMSD for the complexes shows that the theophylline complex with 6LU7 and 6M3M is more stable as compared to the spike protein (6LZG).

Fig. 9
figure 9

RMSD graph of theophylline with SARS-CoV-2 proteins for 100 ns: (A) theophylline complex with 6LZG (B) theophylline complex with 6LU7 and (C) theophylline complex with 6M3M. The duplicate runs are indicated by different line colors as indicated in A, B and C panels

Full size image

RMSF Analysis of Protein–Ligand Complexes

RMSF was done to study the flexibility of protein-ligand complexes and fluctuation in interactive amino acid residues in secondary structure of target proteins (6LZG, 6LU7 and 6M3M) SARS-CoV-2.

The ligand RMSF plot for theophylline fit over the entire 6LZG protein and showed very minor residual fluctuation in alpha helical and beta strands. The fluctuation was observed in ARG 127 and TYR 143 during both the runs (Fig. 10A). In case of 6LU7 complex with theophylline, no residual fluctuation was observed in alpha helical and beta strands in both runs (Fig. 10B). Similarly, in case of 6M3M with theophylline, we observed a perfect fit within the binding pocket of target protein, and residual fluctuation was observed at ALA 91 in both the runs of simulation (Fig. 10C).

Fig. 10
figure 10

RMSF graph of best protein–ligand complexes of theophylline with SARS-CoV-2 proteins for 100 ns: (A) theophylline complex with 6LZG, (B) theophylline complex with 6LU7 and (C) theophylline complex with 6M3M. The duplicate runs are indicated by different line colors as indicated in A, B and C panels

Full size image

Radius of gyration is in the range of 1.8–1.95 nm for 6LZG complex with theophylline, 2.2–2.5 nm for 6LU7 in complex with theophylline and 1.35–1.4 nm for 6M3M in complex with theophylline during both the runs as shown in Figure 11A, B, C. The radius of gyration plots establishes the compactness of the theophylline complexes with the SARS-CoV-2 proteins (6LZG, 6LU7 and 6M3M) and confirms their stable interactions. Solvent accessible range of theophylline complex with 6LZG is between 106 and 125nm, theophylline complex with 6LU7 is between 145 and 160 nm and theophylline complex with 6M3M is 73–84 nm during both the runs as shown in supplementary figures (Fig. S1A, B, C). The changes in the solvent accessible surface area did not change much about a mean position through the course of the simulation of theophylline with the SARS-CoV-2 proteins (6LZG, 6LU7 and 6M3M). This suggest that the interaction of the ligands with the selected proteins did not overtly change the conformation of the protein about the mean position. However, the hydrogen bonding (Fig. S2A, B, C) observed for the three complexes during the course of the simulation showed that theophylline had the most polar interactions with the 6LU7 protein, which must have ultimately contributed to it having the highest estimated free energy of binding. Essentially, the data obtained from the RMSD, RMSF, ROG, SASA and hydrogen bonding estimates are in good agreement with the observed free energy of binding for the three complexes of theophylline with the SARS-CoV-2 target proteins (6LZG, 6LU7 and 6M3M). The complexes are stable, and the changes over the course of the simulation are limited, but significantly different from each other. The free energies are different, and the highest free energy of binding was observed for the theophylline complex with 6LU7. The breakdown of the contributions of all the energy terms to the binding energy is presented in Table 5.

Fig. 11
figure 11

Radius of gyration graph of best protein–ligand complexes of theophylline with SARS-CoV-2 proteins for 100 ns: (A) theophylline complex with 6LZG, (B) theophylline complex with 6LU7 and (C) theophylline complex with 6M3M. The duplicate runs are indicated by different line colors as indicated in A, B and C panels

Full size image
Table 5 Breakdown of the energy contributions to the estimated free energy of binding
Full size table

MD Simulation by Using Desmond Program Version 2.0 (Academic Version)

MD simulation of theophylline in complex with 6LZG, 6LU7 and 6M3M was done in triplicate (100ns) to study the stability of protein-ligand complexes using Desmond program version 2.0 to verify the binding stability.

RMSD of Protein–Ligand Complexes

The RMSD graph of theophylline in complex 6LZG showed the stability of backbone and Cα at 40 ns at 4–4.5 Å and remains stable up to 100-ns time scale. The RMSD of ligand seems to be stable from the beginning of simulation (0.2–0.5 Å) and remains stable up to 100 ns as shown in Fig. 12A. RMSD of theophylline in complex with 6LU7 showed stability in the backbone and Cα at 30 ns at 3.2–3.6 Å. Ligand RMSD was stable from the beginning of simulation as shown in Fig. 12B. In case of theophylline complex with 6LZG, RMSD of backbone (3–3.5 Å) and Cα (2–2.5 Å) was stable from the initial phase of simulation and remains stable up to 100 ns. RMSD of the ligand was stable from 0 to 100 ns at 0.5 Å (Fig. 12C).

Fig. 12
figure 12

RMSD graph of best protein–ligand complexes of theophylline with SARS-CoV-2 proteins for 100 ns: (A) theophylline complex with 6LZG, (B) theophylline complex with 6LU7 and (C) theophylline complex with 6M3M. The different color lines for Cα, backbone, side chains, heavy atoms, ligand fit on protein and ligand fit on ligand are as indicated at the top of panel A, B and C

Full size image

RMSF of Protein–Ligand Complexes

RMSF explained the flexibility of protein structures and fluctuation of interactive residues in secondary structure elements of target proteins. The RMSF plot for theophylline in complex with 6LZG showed low residual fluctuation in α-helical and β-strands (Fig. 13A) and showed hydrophobic interactions with Val 73, Ile 75; hydrogen bonds with Gln 71, Pro 75, Asn 76 and water bridges with Gly 70, Gln 71, Gln 84, Tyr 113. (Fig. 14A). RMSF plot for theophylline fits over the entire 6LU7 protein, and no residual fluctuation was observed in α-helical and β-strands (Fig. 13B). Theophylline in complex with 6LU7 showed hydrophobic interactions with Met 49, Met 165, Leu 167, Ala 173, Phe 185; hydrogen bond interactions with His 41, Cys 44, Asn 142, Gly 143, Cys 145, His 164, Arg 188, Gln 189, Thr 190 and water bridges with Thr 26, Arg 40, His 41, Asn 119, His 164, Glu 166, Val 185, Val 186, Asp 187, Arg 188, Gln 189, Thr 190, Gln 192, Ala 193 (Fig. 14B). Similarly, RMSF of theophylline fit over the entire 6M3M, and no residual fluctuation was observed in α-helical and β-strands (Fig. 14C). Theophylline in complex with 6M3M showed the hydrophobic interactions with Phe12, Leu 122, Leu 125, Phe 126, Lys 128, Tyr 143, Phe 156, Tyr 159, Phe 160, Leu 162, Leu 188; hydrogen bonding with Ser 41, Tyr 91, Arg 124, Lys 128, Ser 129, Glu 129, Glu 141, Ser 147, Asn 151, Glu 154, Gly 155, Phe 156, Lys 156, Asn 157, Lys 158, Tyr 159, Phe 160, Leu 187, His 189 and water bridges with Asn 40, Lys 87, Tyr 91, Arg 124, Leu 125, Arg 127, Lys 128, Ile138, Ser 139, Glu 141, Tyr 143, Gln 144, Ala 445, Thr 148, Pro 149, Asn 157, Lys 158, Tyr 159, Phe 160, Leu 162, Gln 163, Ser 164, Glu 186, Leu 187, His 189 (Fig. 14C). The breakdown of the contributions of all the energy terms to the binding energy is presented in Table 6.

Fig. 13
figure 13

RMSF graph of best protein–ligand complexes of theophylline with SARS-CoV-2 proteins for 100 ns: (A) theophylline complex with 6LZG, (B) theophylline complex with 6LU7 and (C) theophylline complex with 6M3M. The different color lines for Cα, backbone, side chains, heavy atoms and B factor are as indicated at the top of panel A, B and C

Full size image
Fig. 14
figure 14

Histogram of ligand contacts with amino acid residues of target proteins: (A) theophylline complex with 6LZG, (B) theophylline complex with 6LU7 and (C) theophylline complex with 6M3M. The different colors indicating H-bonds, hydrophobic, ionic and water bridges are as indicated at the bottom of panel C

Full size image
Table 6 Breakdown of the energy contributions to the estimated free energy of binding
Full size table

Discussion

Molecular docking and computational techniques are widely used to predict drug likeness, toxicity and protein–ligand interactions for speeding up drug discovery. Molecular docking and in silico approach could help in building the pipeline for wet lab experiments and clinical trials and advance the drug discovery challenges. The current study was focused on major xanthine derivative phytocompounds present in coffee and tea to provide relief and reduce the symptoms of Covid-19 infection and related viruses. The research idea stems from the fact that theophylline (a known methylxanthine) acts as a bronchodilator by relaxing the muscles in the lungs and widening the airways (bronchi). Nasal and lungs are the primary target expressing ACE2 receptor and spike protein of SARS-CoV-2, and related viruses bind ACE2 receptor and initiate pathogenesis.

Molecular docking approach is one technique that helps in identification of interacting molecules. Lower score in molecular docking is considered as higher stability in ligand receptor binding [39]. Using molecular docking approach, we found that caffeine, methylxanthine, theobromine, theophylline, thiene and xanthine are effective in binding affinity or pharmacokinetic properties with three different proteins (6LZG, 6LU7 and 6M3M) of SARS-CoV-2. Among all the selected phytocompounds, theophylline showed most promising interaction with all the three different targets of SARS-CoV-2 with docking scores of − 5.7, − 6.5 and − 5.8 kcal mol−1 with 6LZG, 6LU7 and 6M3M respectively. Similar to our study, Celik et al. (2020) [40••] reported the mechanism of hydroxychloroquine on SARS-CoV-2 proteins with docking scores − 3.074 (6VXX) and − 6.227 (6LZG). Peele et al. (2020) [41] screened 68 US-FDA approved drugs and phytocompounds for binding with SARS-CoV-2 main protease (6LU7), and the best docking score was found with lopinavir, amodiaquine and theaflavin digallate. Bhowmik et al. (2020) reported the antiviral activity of rutin, doxycycline, caffeic acid, ferulic acid, simeprevir and grazoprevir with envelop (E), membrane (M), nucleocapsid (N) proteins of SARS-CoV-2 proteins. They found that rutin and doxycycline showed promising interaction with envelop protein (E) with binding scores of − 9.3 kcal mol−1 and − 9.1 kcal mol−1) respectively. On the other hand, caffeic acid and ferulic acid showed best binding with membrane protein (M)with binding scores of − 8.4 kcal mol−1 and − 8.3 kcal mol−1 respectively. Simeprevir and grazoprevir showed promising binding with nucleocapsid (N) protein with binding scores of − 8.7 kcal mol−1 and − 8.5 kcal mol−1 respectively. Utsunomiya et al. (2008) [42] reported the antiviral activity of hot water extract and instant coffee solution against coffee virus and herpes type 1 virus. Clarck et al. (1998) [43] reported the antiviral activity of crude theaflavin extracted from black tea against coronavirus with EC50 value of 34.7 µg ml−1. Caffeine inhibits the SARS coronavirus main peptidase (6LU7) [44••] and Paine like protease (6W9C) [45]. Theophylline is used as a bronchodilator, and at lower concentration it also showed anti-inflammatory effect in asthma and COPD [9]. Murayama et al. (2008) [46] reported that HSV1, poliovirus and influenza are sensitive to caffeine and have potential to inhibit both DNA and RNA viruses. Theophylline and pyrimidine derivatives are possible inhibitors of RNA bound N-terminal domain [47]. Weber et al. (2020) [48] reported tea infusions inhibited both adenovirus infection and viral protein adenain in HepG2 cells. A report from Clarck et al. (1998) [43] showed the anti-rotaviral activity of black tea theaflavin, whereas the crude tea extract was active against coronavirus. The black tea theaflavin also inhibited the replication of HIV [49]. Polyphenolic compounds of black tea inhibited HIV-1 entry into the target cell by blocking the HIV entry through its envelop (glycoprotein-mediated membrane fusion). Elzupis [50••] studied the inhibitory effect of caffeine and caffeine containing pharmaceuticals known as 3CPs (theophylline, linagliptin, dyphylline, bromotheophylline, pentoxifylline and istradefylline) against SARS-CoV-2 (3-chymotrypsin-like protease 3CLpro). This in silico study showed for the first time that inexpensive drugs available in large quantities could potentially act as inhibitors against 3CLpro. It was also suggested that linagliptin and caffeine could be more specifically used for the treatment of SARS-CoV-2 after in vitro and in vivo studies. LeBlanc and Colpitts [51] proposed that green tea catechin, epigallocatechin gallate (EGCG), prevents murine and human coronavirus infection. Jang et al. [52] showed that tea polyphenols EGCG and theaflavin inhibit the activity of SARS-CoV-2 3CL-protease in vitro. EGCG and theaflavin (active ingredients of green tea and black tea) showed inhibitory activity against the SARS-CoV-2 3CL-protease with half inhibitory concentration (IC50) was 7.58 μg/ml for EGCG and 8.44 μg/ml for theaflavin. No cytotoxicity was observed for either EGCG or theaflavin at the concentrations tested up to 40 μg/ml in HEK293T cells. Bhardwaj et al. [53••] showed that Oolonghomobis flavan, a bioactive molecule found in tea, acts as an inhibitor for the Mpro of SARS-CoV-2. In spite of the tremendous health benefits and potential drug for treating various viral infections, Gottwalt and Tadi [54] have proposed that methylxanthines have a narrow therapeutic range. Mild effects such as nausea, vomiting, increased gastric acid secretion, polyuria, insomnia, palpitations, headaches and tremors could occur even below concentrations of 20 mcg/ml. Severe side effects such as intractable vomiting, arrhythmias, irregular heartbeat (slow or fast), cardiac arrest, allergic skin reactions or seizures could occur with serum concentrations that exceed 20 mcg/ml. Furge and Guengerich [55] have reported Cytochrome P450 may deactivate the caffeine and could reduce its efficacy.  

Conclusion

Caffeine/thiene, methylxanthine, theobromine, theophylline and xanthine were screened on the basis of binding energy against SARS-CoV-2 spike protein S1 and receptor-binding domain (S1RBD). All the selected phytocompounds showed promising binding affinity with SARS-CoV-2 spike protein and S1 receptor-binding domain (S1RBD) and could potentially prevent the binding of virus to ACE2 receptor and stop the infection. ADMET prediction revealed that caffeine/thiene, methylxanthine, theobromine, theophylline and xanthine could be developed as potential drugs for SARS-CoV-2 treatment. Based on our in silico studies, the future challenge is to test all the methylxanthine derivatives using in vitro and in vivo assays for developing broad spectrum therapeutics against different viruses including SARS-CoV-2.