Montelukast is a dual-purpose inhibitor of SARS-CoV-2 infection and virus-induced IL-6 expression identified by structure-based drug repurposing.

Drug-repurposing has been instrumental to identify drugs preventing SARS-CoV-2 replication or attenuating the disease course of COVID-19. Here, we identify through structure-based drug-repurposing a dual-purpose inhibitor of SARS-CoV-2 infection and of IL-6 production by immune cells. We created a computational structure model of the receptor binding domain (RBD) of the SARS-CoV-2 spike 1 protein, and used this model for insilico screening against a library of 6171 molecularly defined binding-sites from drug molecules. Molecular dynamics simulation of candidate molecules with high RBD binding-scores in docking analysis predicted montelukast, an antagonist of the cysteinyl-leukotriene-receptor, to disturb the RBD structure, and infection experiments demonstrated inhibition of SARS-CoV-2 infection, although montelukast binding was outside the ACE2-binding site. Molecular dynamics simulation of SARS-CoV-2 variant RBDs correctly predicted interference of montelukast with infection by the beta but not the more infectious alpha variant. With distinct binding sites for RBD and the leukotriene receptor, montelukast also prevented SARS-CoV-2-induced IL-6 release from immune cells. The inhibition of SARS-CoV-2 infection through a molecule binding distal to the ACE-binding site of the RBD points towards an allosteric mechanism that is not conserved in the more infectious alpha and delta SARS-CoV-2 variants.

Introduction

The Severe Acute Respiratory Syndrome coronavirus 2 (SARS-CoV-2) emerged in late 2019 and developed into a pandemic [1]. SARS-CoV-2 infection can cause disease (COVID-19) that is unpredictable in its course and severe COVID-19 is characterized by failure of multiple organs [2]. Infection of SARS-CoV-2 is mediated by binding of the receptor-binding domain (RBD) of the viral spike protein to its cellular receptor, Angiotensin-converting enzyme 2 (ACE2) [3]. The spike-protein consists of the S1 subunit that entails the RBD which interacts with ACE2 at the binding site, i.e. LYS-391, TYR-425, GLN-448, PHE-460, PHE-472, PHE-474 of the RBD [4]. Activation of the serine protease TMPRSS2 then induces changes in the S2 subunit of the spike protein leading to conformational changes of the spike protein to a pre-fusion state that leads to infection of the cell [3]. Since many cells in the human body express ACE2 and TMPRSS2, SARS-CoV-2 can infect multiple organs. Upon infection, SARS-CoV-2 induces numerous changes in the cell that support its replication [5]. However, critical organ and tissue-damage is not caused by viral replication but rather by an overzealous immune response against SARS-CoV-2 infection that is characterized by widespread immune cell activation and induction of a cytokine storm [6], [7], [8], [9]. Although COVID-19 vaccines have proven to be highly efficient in preventing SARS-CoV-2 infection and also reduce severity of COVID-19 [10], [11], there is still an urgent need for additional therapeutic options against SARS-CoV-2 infection and COVID-19. Therapeutic candidates which are based on monoclonal antibodies suffer high costs with limited outcomes [12]. The continuous development of SARS-CoV-2 variants may pose a threat to the efficacy of COVID-19 vaccines, as has been shown for the SARS-CoV-2 beta variant [13].

Drug repurposing has attracted a lot of attention as a possibility to employ existing drugs with a known safety profile to fight SARS-CoV-2 infection and COVID-19 [14]. The best example is the drug Remdesivir, that was originally developed as antiviral agent against Ebola virus infection, but also proved to be efficient in limiting replication of SARS-CoV-2, when given during the early stages of infection [15]. Drug repurposing approaches can be used on structural analysis searching for similar binding sites, by targeted searches based on transcriptome analyses of SARS-CoV-2-infected cells or by experimental testing of candidate molecules, which have all been employed to repurpose drug candidate molecules for use against SARS-CoV-2 [16], [17], [18]. While all steps of the SARS-CoV-2 life cycle may serve as target for a drug repurposing search, the initial step of SARS-CoV-2 spike binding to its cellular receptor ACE2 bears the advantage to prevent infection of cells. Drug repurposing approaches have predicted several drugs that may interfere with binding of SARS-CoV-2 spike protein to ACE2 or inhibit the protease activity of TMPRSS2. Drugs like Zanamivir, Indinavir, Saquinavir or Remdesivir were predicted to have an impact on SARS-CoV2 infection by in-silicostudies, and other in-silico candidates like Baricitinib or Ivermectin entered clinical trials [19], [20], [21].

Here, we identify the drug montelukast to bind to the RBD of wildtype (WT) and variant SARS-CoV-2 and to interfere with infection of SARS-CoV-2. Furthermore, montelukast also reduced the SARS-CoV-2-induced release of IL-6 from immune cells, thus attributing a dual function of montelukast in reducing SARS-CoV-2 infection and limiting release of the pro-inflammatory mediator IL-6, that has been shown to be involved in severe COVID-19 [22].

Materials and methods

Sequence retrieval, structural modeling, and validation

The amino acid sequence of the SARS-CoV-2 spike (S) surface glycoprotein was retrieved from NCBI protein database (accession no: QHD43416.1). The viral S protein consists of two parts, the S1 protein that contains the receptor binding domain (RBD) engaging the ACE2 (residues: 13–685 residues), and the S2 protein that is involved in cell fusion and binding to TMPRSS2 [3]. We took only RBD protein sequence for further analysis. The 3D structure of RBD of S protein was modelled based on the recently reported experimental structure of spike protein [23] (PDB id:6VSB) using SWISS-MODEL [24]. The SWISS-MODEL gave 11 models we took the 3D model based on the Global Model Quality Estimation (GQME) and Qualitative Model Energy Analysis (QMEAN) values. GQME values range from 0 to 1, higher value represents reliability of the predicted 3D structure, while for QMEAN value < 4.0 shows reliability. The S1 subunit of the computational RBD model was refined using ModRefiner [25], which optimizes atomic details like turning side chains into favorable positions, and Chiron tools [26] was used to remove steric clashes. Molprobity server [27] was used to create Ramachandran plot and z-score [28] of torsion angels of each amino acids and to predict the stereochemical quality of the modeled 3D RBD structure. ProSA-web algorithm [29] was used to compute the Z-score of the RBD structure model to check for potential errors compared to experimentally resolved structures derived from X-ray or NMR analyses. We computed RMSD between the RBD protein structural model and the experimentally resolved partial structure of the RBD from Yan R et al. [4] (PDB id: 6 m17, chain E) using Pymol software. The surface plot of the 3D RBD model was rendered by v2.7 Pymol software.

Binding site similarity analysis and drug identification

To perform for the binding site similarity analysis of RBD of S protein 3D model we utilized the SPRITE program in Drug ReposER (Drug REPOsitioning Exploration Resource) web server [30]. We uploaded the RBD of S protein 3D model to the SPRITE program which matches the 3D model with a database of known drug binding sites. We took the binding sites which showed lowest (≤0.7 Å) RMSD with the S1 subunit of RBD protein model binding sites for further analysis. The druggability and volume of the binding sites was assessed using the PockDrug-Server [31]. The drugs chlorambucil, tigecycline and montelukast were extracted and used for docking and molecular dynamics simulation analysis.

In-silico mutagenesis to generate models of SARS-CoV-2 variant S1-RBD

We downloaded mutations related to alpha, beta and delta variants from (https://covariants.org/variants). Then we generated alpha (ASN->TYR-475, ASN->ASP-544), beta (ARG->ILE-220,ASP->ALA-54,ASP->GLY-189, GLU->LYS-458,ASN->TYR-475, and LYS->ASN-391) and delta (ARG->GLY-131,LEU->ARG-426,THR->LYS-452 and ASP->GLY-588) variants of RBD 3D protein structural models using mutagenesis wizard in PyMOL software v2.7 ().

Molecular docking

Molecular docking analysis was performed using Auto Dock tools [32] and Auto Dock Vina v1.1.2 software [33], the three drugs chlorambucil, tigecycline and montelukast found via binding site similarity analysis were docked on to RBD of S protein 3D model. First, we downloaded the 3D structures of drug molecules in .SDF format from PubChem database (). We utilized openbabelv2.4 to build ionization states of chlorambucil (pH 5.8), tigecycline (pH 7.3) and montelukast (pH 5.7) [34], [35], [36]. These drug molecules in .MOL2 format were loaded into the Auto Dock Tools software and their torsions were adjusted based on the total number of rotatable bonds and saved in the .pdbqt format. Next, we added polar hydrogen partial charges and Auto Dock atom types were assigned to the RBD model. The Auto Grid 4 module was used to build three grid boxes (three binding sites in RBD model) 119.432 × 102.140 × 101.264 Å, 66.570 × 48.380 × 126.724 Å and 126 × 116 × 100 Å which were constructed to comprise the entire binding sites of wildtype, alpha, beta and delta variants of the S1-RBD structure. The in-silico docking analysis was performed between the experimentally resolved protein structure of the cysteinyl leukotriene receptor (PDB ID: 6RZ4) and montelukast using the grid box (33.988 × 30.777 × 23.884), which encompass the protein binding site. Finally, binding affinities were computed using Auto Dock Vina v1.1.2 software [33]. We used grid spacing 0.375 Å and exhaustiveness was set at 20. Docked poses were visualized using Pymol v2.7 software (). The discovery studio visualizer was used to visualize the molecular interactions between the drug and protein structure.

Molecular dynamics simulation

To perform molecular dynamics (MD) simulation we extracted the both wildtype and variant RBD 3D models and drug molecules from docked poses. Simulations were carried out using GROMACS 2020 package [37] and the CHARMM27 force field was adopted for the RBD model in GROMACS. Force fields for drugs were prepared using the SwissParam web server [38]. In total we performed three simulations, which included two S1-RBD-drug complexes and one apo form of S1-RBD without a drug-molecule that served as a negative control during simulations. We used pdb2gmx tool to add hydrogens to the RBD protein structure. The S1-RBD (apoform) and the S-1 RBD-drug complexes (bonded forms) were engaged in a water box using the explicit TIP3 water model at a buffering distance of 1.2 nm. To neutralize the net charges in the system NA+ and CL+ counterions were added to the system to attain equilibration. The S1-RBD drug complexes and apoform of RBD were subjected to minimization procedure with the Steepest Descent method for 2000 steps. All simulations were run in periodic boundary conditions with NVT and NPT ensemble. The Berendsen coupling algorithm was used for this process, and the temperature was kept at a constant 310 K with pressure at 1 bar. The Particle-Mesh Ewald (PME) is used to compute the electrostatic interactions with an interpolation of order 4 and grid spacing of 0.12 nm, and all bonds were constrained using the LINCS algorithm. Finally, the time step for all simulations was set at 2 fs, and 100 ns of molecular dynamic simulations were performed. An analysis of molecular dynamic trajectories was performed with built-in GROMACS tools. The program RMS was used to analyze the stability of apo-form of RBD and drug complexes. RMSF was used to analyze the fluctuation of each amino acid residue in the RBD. Radius of gyration analysis was performed using the gyrate program. We calculated the binding free energy at every 20 ns interval for protein-drug complexes using the molecular mechanics energies combined with Poisson-Boltzmann (MM-PBSA) using the g_mmpbsa tool [39].

Drug compounds

Montelukast (5 mg, PZN 09731805, 1A Pharma, Germany) was dissolved in DMSO (Sigma-Aldrich, Germany) at a concentration of 20 mM, sterile filtered (0.22 µm) and stored in aliquots at −20 °C until usage. Tigecycline (Y0001961, Sigma-Aldrich, Germany) was dissolved in DMSO and at a concentration of 20 mM, sterile filtered (0.22 µm) and stored in aliquots at −20 °C until usage. For the neutralization assay and the cytokine array, montelukast sodium CRS (European Pharmacopoeia Reference Standard, Y0001434, batch 2.1, EDQM, France) was dissolved in DMSO at a concentration of 20 mM, sterile filtered and stored in aliquots at −20 °C until usage. Chlorambucil was dissolved (PZN 01263654, aspen pharma, Germany) in 0.5% DMSO at a concentration of 400 µg/ml, sterile filtered and stored in aliquots at −20 °C until usage.

SARS-CoV-2 amplification

SARS-CoV-2-Munich-TUM-1 (EPI-ISL-582134) derived from patient material in Germany was propagated and the viral titer was determined by plaque assay, as described previously [40]. Briefly, SARS-CoV-2 was isolated from a patient who was infected during the COVID-19 outbreak in Germany in January 2020 with a virus imported from Wuhan via a single contact in Shanghai [41]. Different SARS-CoV-2 variants were also used in the study, alpha variant was isolated from a patient in Bundeswehr and kindly provided by Institute of Microbiology, TUM. The Beta SARS-CoV-2 variant was isolated as per the instructions provided by Bayerisches Landesamt für Gesundheit und Lebensmittelsicherheit. All SARS-CoV-2 isolates were propagated in Vero E6 cells in DMEM medium with 5% FCS, 200 mM L-glutamine, 1% MEM-non-essential amino acids, 1% sodium pyruvate and 1% penicillin/streptomycin (all Gibco, Gaithersburg, MD, USA). For virus titration, Vero E6 cells were seeded in 12-well plates with 5 × 105 cells/well in supplemented DMEM medium and infected after overnight incubation with serial dilutions of SARS-CoV-2 for one hour at 37 °C. Thereafter, cells were incubated with 1 ml per well of 5% carboxymethylcellulose (Sigma) in Minimum Essential Media (Gibco) at 37 °C until plaques were visible. After fixation with 10% paraformaldehyde (ChemCruz, US) for 30 min at room temperature and washing with PBS, 1% crystal violet diluted in 20% methanol (Sigma) was added for 15 min at room temperature. After PBS washing, the plate was dried and plaque forming units (PFU) were determined by counting the average number of plaques per dilution and considering the dilution factor.

SARS-CoV-2 infection inhibition assay

Viral neutralization followed by in-cell ELISA was performed as described previously with a multiplicity of infection (MOI) of 0.03 PFU/cell. Briefly, 1.6 × 104 Vero E6 cells per well of a 96-well plate were seeded one day before infection with SARS-CoV-2 in DMEM supplemented as described above for SARS-CoV-2 propagation. Serial concentrations of montelukast or chlorambucil were diluted in supplemented DMEM and incubated with SARS-CoV-2 for 1 h at 37 °C and added thereafter onto the cells at an MOI of 0.03 PFU/ml. After inoculation for 1 h at 37 °C, supernatants were removed and replaced by supplemented DMEM medium. 23.5 h after infection of Vero E6 cells with SARS-CoV-2, a CellTiter-Blue Cell Viability Assay (Promega, Fitchburg, WI, USA) was performed. CellTiter-Blue reagent was added to cell culture medium in a 1:6-dilution and incubated for 25 min at 37 °C. Fluorescence was measured on a Tecan infinite F200 plate reader at 550 nm (excitation) and 590 nm (emission). Cells were washed once with PBS thereafter and the in-cell ELISA was continued with fixation 24 h after infection by adding 4% paraformaldehyde for 15 min. After washing with PBS, cells were permeabilized for 15 min at room temperature using 0.5% saponin (Roth) in PBS and blocked with 0.1% saponin and 10% goat serum (Sigma) in PBS for 1 h at room temperature. Then, primary antibody SARS-CoV-2 Nucleocapsid antibody T62 (40143-T62, Sino Biological, Beijing, China) in 1% FCS in PBS was added in a 1:1500-dilution at room temperature for 2 h. After 3 times washing (PBS with 0.05% Tween-20 (Roth)), secondary antibody goat anti-rabbit IgG-HRP (12–348, EMD Millipore, USA) in 1% FCS in PBS in a 1:4000-dilution was added for 1 h at room temperature. After 5 washing steps (PBS with 0.05% Tween-20), in-cell ELISA was developed with TMB for 10 min, stopped with 2 N H2SO4 and optical density was measured on a Tecan infinite F200 pro plate reader at 450 nm (correction with 560 nm). Percent of neutralization was calculated by setting uninfected control cells to 100% and infected control cells to 0% after background subtraction. The results were analyzed with the software PRISM version 8 (GraphPad Software, USA), and the drug concentration, where 50% of inhibition (IC50) was achieved, was calculated by non-linear regression.

Cytokine measurements in immune cells after incubation with SARS-CoV-2

Human PBMCs were isolated from peripheral blood of healthy volunteers. PBMCs were separated by Biocoll density gradient centrifugation (Biochrom, Germany) at 800g for 20 min without break. PBMCs were collected and washed twice with PBS (350g, 5 min) and resuspended in RPMI supplemented with RPMI, 10% FCS, 1% penicillin/streptomycin, 2 mM L-glutamine, 1% sodium pyruvate, 1% non-essential amino acids, 10 mM HEPES, 20 µg/ml gentamicin. After overnight resting, 1x106 PBMCs per replicate were incubated for 24 h with SARS-CoV-2 at an MOI of 0.1 PFU/cell in presence or absence of 200 µM montelukast. Cell culture supernatants were harvested and subjected to centrifugation before analysis through a cytokine array. A Cytokine Human Membrane Antibody Array (ab133997, Abcam, Cambridge, UK) was used according to manufacturer’s instructions using pooled supernatants from biological replicates.

Results

Homology modeling of RBD spike (S) protein of SARS-CoV-2

To identify novel molecules that prevent SARS-CoV-2 infection, we chose an integrated computational analysis pipeline to make use of an existing molecular data set on drug-target interactions (Supplementary Fig. S1A). SARS-CoV-2 infects cells by binding to its receptor ACE through its spike protein [3] The wild type SARS-CoV-2 spike protein consists of 1273 amino acids, and is composed of two subunits: the S1 subunit contains the receptor binding domain (RBD) from aa 13–685 (673 amino acid residues) (Fig. 1A, Supplementary Fig. S1B) and the S2 subunit facilitates fusion to the cell membrane. While the protein structure of the RBD and of the spike protein were resolved experimentally by Cryo-EM [4], [23], we found that some residues were missing from this analysis, i.e. the RBD structure covers only aa 316 to 518 and the full length spike structure covers only aa 27 to 672 with several amino acid clusters further missing (Fig. 1B). We therefore generated a full length structure of RBD using a homology modelling approach to predict the 3-dimensional (3D) protein structure of RBD using the SWISS-MODEL software [24]. We used the already known spike protein structure as template to model the full-length RBD protein structure. This approach allowed us to generate a computational 3D model of spike 1 RBD (S1-RBD) covering aa 27 to aa 672 (Fig. 1C, Supplementary Movie 1). The model revealed particular regions, i.e. the ACE2-binding site lies in a loop region with two adjacent beta sheets, further beta sheets and an alpha helix region 1 in proximity to the ACE binding site, two beta sheet regions and one loop-alpha helix beta sheet region are located at the hinge region, one beta sheet, alpha helix and loop region next to the S1-RBD structure.

Fig. 1

Structure analysis of the receptor binding domain (RBD) of SARS-CoV-2. (A) full sequence of the SARS-CoV-2 (1.177) spike protein with the S1 protein sequence (S1-RBD) highlighted. (B) Sequence features highlighting the missing structural data in published protein structures of S1-RBD and spike protein. (C) structural model of the S1-RBD by computational analysis of the aa 13–685 amino acids; ACE2 binding residues of the S1-RBD are highlighted. (D) validation of the computational S1-RBD structural model by ProSA (z-score: −8.37) using X-ray- and NMR-based experimental structures. (F) Ramachandran plot shows that 96.8% of the S1-RBD structural model is in a favored region (cyan lines), 3.2% allowed region (purple blue lines) and 0% unfavorable region (white space outside cyan and dark blue lines); Rama Z-score (-0.04) of residue distributions in the Ramachandran plot. (G) The overlap between the S1-RBD computational model (cyan color) and experimental structure of RBD (orange color) with an RMSD 0.75 Å… (Å). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We performed a quality check for our computational S1-RBD model that revealed a GQME value of 0.65 and a QMEAN score of −2.79, indicating that the computational RBD structure was reliable. To further improve the computational S1-RBD model with respect to hydrogen bonds, backbone topology and sidechain positioning we performed energy minimization using modrefiner [25]. As in-silico 3D models tend to have abnormal steric clashes, that may produce unfavorable bond lengths, it is important to further optimize the structural model of S1-RBD to remove steric clashes [26]. After this we subjected the improved protein structure of S1-RBD to a verification process to assess the its quality. We used the ProSA algorithm [29], which compares the predicted protein structure with experimentally resolved protein structures through x-ray and NMR techniques and assigns a score z-score from 10, −20 to quantify the overall model quality, which evaluates the deviation of the total energy by the predicted structure. The ProSA analysis of the computational S1-RBD structure revealed a z-score −8.35 (Fig. 1D), which indicates that it fits very well with experimentally resolved protein structures.

To evaluate the quality of the computational S1-RBD structure in more detail, we performed an analysis at the level of single residues using Ramachandran plot analysis, which assesses favored and allowed regions of protein residues in the context of steric hindrance in the protein backbone by calculating phi and psi torsion angles. These angles define the rotations of the backbone bonds between the N terminal-Cα (Phi) and Cα- C (Psi) terminus of protein structures and shows the quantitative distribution of these angles [42]. The Ramachandran plot analysis revealed a high number of residues (96.8%) in the computational RBD structure model to be in an energetically favorable region (Fig. 1E), which suggested that these residues had minimum steric hinderance. Only a small percentage of residues (3.2%) in the computational S1-RBD structure model were in allowed region, which indicates low level steric hinderance and most importantly, there were no residues found to be in unfavorable regions (Fig. 1E). Recently, the z-score calculated from the Ramachandran plot analysis (Rama-Z score) was proposed as a validation metric to locate abnormal distributions of the protein structure conformation [28], with a Rama-Z score > 3 indicating protein structure backbone with unlikely geometry and a Rama-Z score < 2 indicating a favorable backbone geometry. The Rama Z-score calculated for the computational S1-RBD structure model showed revealed a very low Rama-Z score of −0.04 (<2) (Fig. 1E), further confirming that the geometry of computational S1-RBD structure was reliable. Moreover, the overlay of the computational S1-RBD structure with the experimentally resolved RBD structure [4] revealed a perfect overlap with very low root mean square deviation (RMSD) of 0.75 a (<2 Å) (Fig. 1F). Thus, the computational S1-RBD structure model generated here appeared to reflect the structural confirmation of the RBD protein, and allowed us to use this model to perform in-silico screens to discover molecules that prevent to ACE2 and thereby prevent infection.

Identification of candidate molecules from the ReposER database for RBD-binding by repositioning and molecular docking analyses

Drug repurposing is a powerful approach to identify new uses for approved drugs that are outside the scope of the original medical indication [14], [43]. We used our computational S1-RBD structure model to perform an in-silico screen using the known binding sites of drug molecules contained in the DrugReposER database [30] and employing the SPRITE program to determine binding site similarities between the originally identified binding site and the S1-RBD binding site. This Calculation of binding site similarities of the 6171 structurally defined binding sites of drug molecules in the DrugReposER database to the computational S1-RBD structure model identified 199 candidate drug molecules and their binding sites in the RBD (Fig. 2A and Supplementary Table S1). Structural similarity at the level of a particular binding site was quantified by root mean square deviation (RMSD) with low RMSD values denoting high similarity. This analysis showed RMSD values that ranged between 0.52 and 1.41 Å for the binding to RBD from binding sites of drug molecules contained in the DrugReposeER database (Fig. 2B). Among the candidate molecules with RMSD ≤ 0.7 Å indicating strong binding were adenosine, montelukast, chlorambucil, benzylpenicillin, niacin, tobramycin, bortezomib, mycophenolic acid, tigecycline that all bound at different sites to the RBD (Fig. 2C). We decided to further investigate three candidate drug molecules: montelukast, an inhibitor of the leukotriene receptor CysLT1 [44]; chlorambucil, a chemotherapeutic agent against hematological tumors [45]; and tigecycline an anti-microbial agent [46], because of their low RMSD values (montelukast and chlorambucil) and their distinct binding site (tigecycline binding to the ACE2-binding region of the RBD). We used isocaloric titration and NMR-based structural analysis to confirm binding of the candidate compounds to the RBD protein. However, due to the low solubility of the compounds in water it was not possible to obtain reliable results.

Fig. 2

Drug repositioning analysis. (A) drug repositioning analysis for S1-RBD based on binding site similarity. (B) Histogram of numbers of molecules with similar binding sites in computational S1-RBD structure with protein binding sites. (C) Drug molecules extracted from the database with <0.7 Å. (D) Calculated binding affinities (Ki) using molecular docking analysis of drug molecules and the RBD binding regions. (E-G) Illustration of predicted chlorambucil, tigecycline and montelukast binding to the computational S1-RBD structure model and insert showing the respective binding sites in the S1-RBD at higher resolution.

We next determined the exact molecular binding sites within the computational S1-RBD structure model and the binding affinities of three candidate drug molecules using an in-silico molecular docking analysis (Fig. 2D). Chlorambucil (C14H19CI2NO2) (Supplementary Fig. S2A) was predicted to bind at ARG-8, GLY-9, GLY-242, GLU-272, THR-273 (Fig. 2E); tigecycline (C29H39N5O8) (Supplementary Fig. S2B) was predicted to bind at SER-447, SER-451, PRO-453, ALA-449, GLY-450 (Fig. 2F), within the ACE2 binding site region of the RBD; montelukast (C35H36CINO3S) (Supplementary Fig. S2C) was predicted to bind at ARG-331, ASN-368, VAL-369, THR-404, SER-488, a region close to ACE2-binding site of RBD (Fig. 2G). We utilized the PockDrug server, which estimates the ligand proximity and surface amino acids and predicts druggability probability score ranges from 0 to 1 and > 0.5 considered to be druggable and it also predicts the volume of the binding site. We employed this server to assess the druggability score and volume of the drug binding sites used in our study. This analysis revealed highest druggability score (P = 0.98) for montelukast binding site (Supplementary Table S2A) in comparison to chlorambucil (P = 0.76) (Supplementary Table S2B) and PockDrug did not predict the druggability score for tigecycline binding site in RBD. The binding affinity scores predicted by in-silico docking analyses were −4.8 kcal/mol for chlorambucil, −6.0 kcal/mol for montelukast and −7.5 kcal/mol for tigecycline. We subjected these docked complexes to ligand interaction analysis this revealed that ARG-8, GLY-9 and GLY-242 forms van der Waals interactions with chlorambucil (Supplementary Fig. 3A). The SER-451 forms hydrogen bonds and ALA-449 and PRO-453 forms van der Waals interactions with tigecycline (Supplementary Fig. 3B). The ARG-331, THR-404, SER-488 forms van der Waals interactions with montelukast (Supplementary Fig. 3C). Together, these in-silico results suggested firm binding of the three candidate molecules to the S1-RBD structure model.

Molecular dynamics simulation of the impact of candidate molecules on S1-RBD structure

In a next step, we addressed whether binding of candidate drug molecules at the different sites within the RBD would alter the S1-RBD structure. For this, we performed molecular dynamics simulation analysis, which was carried out for a simulated time frame of 100 ns to predict stability and conformational changes of the computational S1-RBD structure following binding the candidate molecules using gromacs tools (radius of gyration, RMSD, RMSF).

First, we analyzed whether binding of the candidate molecules caused changes in compactness of the S1-RBD structure by quantifying the radius of gyration (Rg). The native (apo) RBD structure showed an initial Rg of 3.8 nm, and remained stable throughout the rest of the simulation time. The S1-RBD structure in complex with chlorambucil or tigecycline showed similar or even lower Rg values, i.e., 4 and 4.2 nm, respectively (Fig. 3A). The S1-RBD in complex with montelukast, however, showed a rise to an Rg of 4.5 nm and declined after 60 ns still this complex showed highest gyration values in comparison to the other complexes (Fig. 3A), suggesting that montelukast-binding had an influence on the S1-RBD structure.

Fig. 3

Molecular dynamics simulation of the changes in S1-RBD structure after binding of candidate molecules. (A) Radius of gyration (Rg) of S1-RBD model in complex with the candidate molecules over 100 ns. (B) Root mean square deviation (RMSD) of the S1-RBD model after binding of candidate molecules over 100 ns. (C-E) Root mean square fluctuation (RMSF) profile depicting the effect of chlorambucil(C), tigecycline(D) and montelukast(E) on the S1-RBD model at level of single aa residues. (F) Boxplots represent the RMSF distribution of montelukast and ACE2 binding sites. Significance was calculated using the non-parametric Mann-Whitney test. The boxes indicate the 25th percentile, median and 75th percentile.

The molecular dynamics simulation analysis further allowed us to determine the root-mean-square deviation (RMSD) that predicts the distance between atoms in the RBD backbone structure, for S1-RBD in complex with the candidate molecules compared to its native (apo) state. The RMSD of the S1-RBD - chlorambucil complex rose to 1.5 nm at 20 ns, but remained stable thereafter with similar patterns as the native S1-RBD (Fig. 3B), suggesting a low impact on the overall S1-RBD structure. The RMSD of the tigecycline- S1-RBD complex showed a rise to 1.2 nm showed slight fluctuations throughout the simulation (Fig. 3B), also indicating a low impact on the S1-RBD structure. The RMSD of the S1-RBD - montelukast complex, however, exhibited an initial steep rise to 2.3 nm until 25 ns, after which the complex remained stable (Fig. 3B), suggesting a more profound disturbance of the S1-RBD structure. Furthermore, we investigated ligand stability along MD trajectory in each complex we found low RMSD of drug molecules ranging from 0.15 to 0.3 (Supplementary Fig. 4A-C)

We concluded the assessment of the changes of S1-RBD structure after binding of the candidate drug molecules by defining the average RMSD over time, i.e. the root-mean-square fluctuation (RMSF), for the entire range of aa residues in the S1-RBD. This analysis showed that chlorambucil did not induce fluctuations of the residues at its binding-site to RBD nor of the residues forming the ACE2 binding-site within the RBD. Rather, chlorambucil binding appeared to stabilize the S1-RBD structure (Fig. 3C). The RMSF calculated for the S1-RBD in complex with tigecycline also induced only minor changes in the residues forming the ACE2 binding-site (Fig. 3D). However, higher RMSF values were observed for the binding sites of montelukast to S1-RBD, but the RMSF values for the residues in the ACE2 binding-site were not increased (Fig. 3E-F).

Furthermore, we also assessed the free energies of binding of complexes of chlorambucil, tigecycline and montelukast using molecular mechanics – Poisson Boltzmann surface area (MM-PBSA). It has been reported that MM-PBSA delivers truthful estimates of free energies of protein-drug complexes [47]. The MM-PBSA method computed following energy terms van der Waals, Solvent Accessible Surface Area (SASA), electrostatic, polar solvation and final binding energy terms. We calculated these terms from 20 ns to 100 ns with 20 ns interval. We observed that except polar solvation energy all other terms contributed to interaction between the drugs and RBD. The chlorambucil – RBD complex revealed final binding energy score −98.521 KJ/mol. Tigecycline-RBD complex displayed −136.427 KJ/mol final binding energy. The montelukast-RBD complex exhibited final binding energy −163.589 KJ/mol. (Supplementary Fig. 5A-C).

Together, the molecular dynamics simulation analysis suggested that montelukast induced overall changes in the S1-RBD structure with strong binding during simulation, induced distinct changes at its binding site to RBD but had no disturbing impact on the ACE2-binding-site, whereas only minor effects on the S1-RBD structure were predicted after binding of chlorambucil or tigecycline.

Prevention of SARS-CoV-2 infection by montelukast but not chlorambucil and tigecycline

To address the functional relevance of the three candidate drugs against SARS-CoV-2 infection, we performed in vitro infection inhibition experiments with wildtype SARS-CoV-2 (1.177). To this end, Vero cells expressing ACE2 as receptor for SARS-CoV-2 infection were incubated with infectious SARS-CoV-2 and different concentrations of the candidate drug molecules ranging from 1 µM to 1 mM. There was no direct cell toxicity after exposure to the candidate molecules (Fig. 4A). After 24 h, intracellular SARS-CoV-2 RNA was quantified using an in-cell ELISA to determine the inhibitory concentration of the candidate drug molecules, where 50% of infection of cells was blocked (IC50). Chlorambucil did not have a significant inhibitory effect on SARS-CoV-2 infection (Fig. 4B). Notwithstanding the predicted binding of tigecycline to the ACE2 binding-site of the RBD, the infection inhibition assay, which employed infectious SARS-CoV-2, failed to provide any evidence for tigecycline-mediated inhibition of infection. Strikingly, montelukast prevented SARS-CoV-2 infection albeit only at a high IC50 values above 137.3 µM (Fig. 4B), although the molecular dynamics analysis did not predict structural perturbation of the ACE2 binding-site. Taken together, montelukast but not chlorambucil or tigecycline inhibited SARS-CoV-2 infection and indicated that binding of montelukast to the RBD at a site distinct from the ACE2-binding site perturbed infection, which points towards an allosteric mechanism in infection inhibition.

Fig. 4

Montelukast inhibits SARS-CoV-2 infection in vitro. (A) Viability of the Vero cells cell-culture showing that neutralization is not induced by increased cell-death. (n = 3). (B) Infection-inhibition assay for candidate molecules employing SARS-CoV-2 for infection of Vero cells measured by In-cell-ELISA relative to uninfected cells. Data from three independent experiments (n = 3 biological replicates per experiment) or one representative experiment (n = 3 biological replicates).

Docking and molecular dynamics simulation analysis predict distinct changes for the S1-RBD structure from SARS-CoV-2 variants after binding of montelukast

Since SARS-CoV-2 has continuously undergone mutations over the last 18 months [13], we wondered whether molecular dynamics simulation would predict changes in the S1-RBD structure of variant SARS-CoV-2 after montelukast binding, that might be relevant for prevention of infection. We therefore created in-silico mutated S1-RBD models from the alpha (B1.1.7) and the beta (B1.351) variants of SARS-CoV-2 for subsequent docking and molecular dynamics simulation. We replaced the ASN to TYR at position 475 and ASN to ASP at position 544 to create the computational S1-RBD structure of the alpha variant of SARS-CoV-2 (Fig. 5A). For montelukast binding to the S1-RBD of the alpha SARS-CoV-2 variant we calculated a binding affinity of −7.1 Kcal/mol (Fig. 5A), compared to the binding affinity of −6.0 Kcal/mol to the native S1-RBD. This suggested continuous strong binding of montelukast to the S1-RBD of variant SARS-CoV-2, which was not surprising as the mutations in the RBD of the alpha SARS-CoV-2 variant were outside of the montelukast binding site [13]. We further created a computational S1-RBD model for the beta-variant of SARS-CoV-2 by changing ARG to ILE at position 220, ASP to ALA at position 54, ASP to GLY at position 189, ASN to TYR at position 475, LYS to ASN at position 351 and GLU to LYS at position 458 (Fig. 5B). Molecular docking analysis predicted montelukast binding affinity to this S1-RBD of −6.5 Kcal/mol (Fig. 5B). Furthermore, the montelukast interacts with ARG-331, THR-404 and SER-488 via van der Waals interaction in alpha variant (Supplementary Fig. 6A). Drug interacts with the beta variant through the hydrogen bond with ARG-331 and THR-404 and SER-488 forms van der Waals interactions (Supplementary Fig. 6B). Together these docking analyses suggested that binding of montelukast to the S1-RBD of the SARS-CoV-2 variants was not substantially changed.

Fig. 5

Molecular dynamics simulation analysis of RBD from alpha and beta SARS-CoV-2 variants in complex with montelukast. (A, B) Docking analysis of montelukast binding to the S1-RBD of alpha (B1.1.7) and beta SARS-CoV-2 (B.1.351) variants. (C, D) Molecular dynamics simulation analysis for changes in S1-RBD structure of SARS-CoV-2 variants after montelukast binding by predicting the radius of gyration for 100 ns. (E, F) Root mean square deviation (RMSD) of S-1RBD from SARS-CoV-2 variants in complex with montelukast for 100 ns. (G, H) Root mean square fluctuation (RMSF) for S1- RBD from variant SARS-CoV-2 in complex with montelukast over 100 ns. Boxplots of the distribution of RMSF with significances for the differences between the montelukast and ACE2 binding sites and variant sites. Significance was calculates using the non-parametric Mann-Whitney test. The boxes indicate the 25th percentile, median and 75th percentile.

We continued to analyze by molecular dynamics simulation whether S1-RBDs from variant SARS-CoV-2 would behave similarly after montelukast binding. Determination of gyration (Rg) revealed that there was no substantial lasting change in the compactness of the S1-RBD structure from the alpha-variant SARS-CoV-2 in complex with montelukast (Fig. 5C). The compactness of the S1-RBD structure of the beta-variant SARS-CoV-2, however, was changed when in complex with montelukast (Fig. 5D), which suggested that the S1-RBD structure from the alpha variant of SARS-CoV-2 was less susceptible to changes inflicted by binding of montelukast compared to the original SARS-CoV-2 or the beta variant. The RMSD analysis further confirmed this notion, because no substantial change in the RMSD of the S1-RBD from the alpha variant SARS-CoV-2 in complex with montelukast. In contrast, molecular dynamics simulation analysis predicted a change in the RMSD of the S1-RBD from the beta variant SARS-CoV-2 upon montelukast binding (Fig. 5E, F). Interestingly, ligand RMSD analysis showed relatively high value (0.4–0.5 nm) of montelukast in comparison to alpha variant (Supplementary Fig. 7A). Whereas, its value is 0.2–0.25 nm confirming that the ligand is more stable when it is bound to beta variant (Supplementary Fig. 7B). Analyzing the RMSF predicted that binding of montelukast rather stabilized the fluctuations of the residues of the drug-binding sites and the ACE2-binding site of the alpha-variant S1-RBD (Fig. 5G), whereas it increased fluctuations in the beta-variant S1-RBD (Fig. 5H). We also computed binding free energies of montelukast alpha and beta complexes. This revealed montelukast binds to beta with higher binding energy (-130.246 KJ/mol) when compared to alpha −120.071 KJ/mol (Supplementary Fig. 8A-B). Of note, we also generated an in-silico mutated model of the delta variant (B1.617.2) RBD of SARS-CoV-2 and predicted changes in the S1-RBD structure upon montelukast binding docking this revealed SER-488 and THR-404 forms van der Waals interactions with drug molecule (Supplementary Fig. S9A-B). Molecular dynamics simulation, which did not reveal significant changes (Supplementary Fig. S9C-F), and free binding energy analysis showed that montelukast bound to delta variant with −113.761 KJ/mol comparable to the findings for the alpha variant of SARS-CoV-2 (Supplementary Fig. S9G). Together, this predicted disturbance of the beta-variant S1-RBD structure after montelukast binding that was not conserved in alpha and delta SARS-CoV-2 variants.

Montelukast inhibits infection with the beta but not the alpha variant of SARS-CoV-2

To investigate whether the predicted changes in the S1-RBD structure of the SARS-CoV-2 variants after montelukast binding correlated with the capacity to prevent infection, we performed infection inhibition experiments using alpha (B1.1.7) and beta (B.1.351) variant SARS-CoV-2. As detected for infection with wild SARS-CoV-2 infection, there was no direct toxicity of the candidate drug molecules in combination with SARS-CoV-2 infection on Vero cells (Fig. 6A, B). As expected, neither incubation with chlorambucil nor tigecycline reduced infection of alpha and beta variant SARS-CoV-2 of Vero cells (Fig. 6C, D). Strikingly, montelukast did not inhibit infection with the alpha-variant SARS-CoV-2 anymore (IC50 > 700 µM) (Fig. 6C, D). In contrast, montelukast showed similar infection inhibition efficacy (IC50 of 164.6 µM) for the beta variant of SARS-CoV-2 (Fig. 6C, D) as compared to the wildtype SARS-CoV-2 (IC50 > 130 µM). Together these results supported the notion that changes in the overall RBD structure predicted by high RMSD values in molecular dynamics simulation were also predictive for inhibition of infection with SARS-CoV-2 variants.

Fig. 6

Montelukast but not chlorambucil and tigecycline inhibits infection of the beta variant SARS-CoV-2. (A, B) Viability of the Vero cells cell-culture showing that neutralization as seen in A-C is not induced by increased cell-death. (C, D) Neutralization of SARS-CoV-2 infection (Alpha and Beta) of Vero cells by montelukast or chlorambucil measured by In-cell-ELISA relative to uninfected cells. Data from three independent experiments (n = 3 biological replicates per experiment) or one representative experiment (n = 3 biological replicates).

Montelukast reduces SARS-CoV-2-induced IL-6 production by peripheral blood immune cells

Since the original target of montelukast is the cysteinyl-leukotriene-receptor, where montelukast functions as antagonist and prevents leukotrienes from triggering cell activation and production of inflammatory mediators such as IL-6, we reasoned that montelukast might also reduce SARS-CoV-2 induced cytokine expression by immune cells. In-silico docking analysis of montelukast and cysteinyl leukotriene receptor (PDB ID: 6RZ4) revealed that it is bound with -8.8 Kcal/mol binding affinity and showed the montelukast binding sites in the receptor(Supplementary Fig. 10A). Clearly, montelukast binding sites to the RBD of SARS-CoV-2 and the cysteinyl leukotriene receptor were distinct and at the opposite ends of the molecule (Fig. 7A). This led us to investigate whether montelukast inhibited SARS-CoV-2-induced release of cytokines from human peripheral blood mononuclear cells (PBMCs) using a multiplex cytokine array. SARS-CoV-2 is known to trigger innate immune cell activation [48]. Incubation of human PBMCs with SARS-CoV-2 led to increased expression of IL-6 (Fig. 7B, Supplementary Fig. S10B-E). There was constitutive IL-8 expression presumably from PBMC contact with tissue culture material that was further increased by incubation with SARS-CoV-2 (Fig. 7B). SARS-CoV-2 induced release of IL-6 and IL-8 was partially inhibited by montelukast (Fig. 7B). This indicated a role of montelukast in controlling IL-6 release through its antagonist function of the leukotriene receptor, although we cannot formally exclude that a reduction of infection was responsible.

Fig. 7

Prevention of SARS-CoV-2 induced IL-6 release from human PBMCs by montelukast. (A) Schematic illustration of montelukast with residues relevant for binding of the RBD domain (cyan) and the cysteinyl leukotriene receptor (green). (B) Bead-array based quantification of cytokine release from peripheral blood mononuclear cells (PBMCs) from healthy volunteers exposed in vitro to SARS-CoV-2 for 24 hrs in presence or absence of montelukast. Data from three independent experiments are shown (n = 3 biological replicates per experiment) or one representative experiment (n = 3 biological replicates). Statistical significance was calculated using two-way-anova with p* = 0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion

Here, we used a structure-based drug-repurposing approach to screen in a large database of binding sites for drug molecules to identify candidate molecules to interfere with SARS-CoV-2 infection. Following homology modeling of the SARS-CoV-2 RBD to generated a complete computational S1-RBD structure model, that was further optimized by atomic-level energy minimization and correction of steric clashes, we used a structure-based drug repurposing approach using a publicly available molecular interaction database to identify candidate drug molecules with predicted high RBD binding affinity. Molecular dynamics simulation analysis predicted only one candidate molecule, i.e., montelukast– a cysteinyl leukotriene receptor antagonist, to cause structural changes of the S1-RBD molecule of the original SARS-CoV-2 (B1.117) and the beta variant of SARS-CoV-2 (B1.351) but not the alpha variant SARS-CoV-2 (B1.1.7). Notwithstanding its binding site to S1-RBD outside of the ACE-binding region, montelukast inhibited in vitro infection with wildtype and the beta but not the alpha variant of SARS-CoV-2. In addition to its antiviral activity, montelukast also reduced the expression of the cytokine IL-6 by human peripheral blood-derived immune cells exposed to SARS-CoV-2.

Large data platforms have been generated that cover the entire range of possible anti-viral drug targets for SARS-CoV-2 [18], [49], [50] and a detailed molecular landscape of the virus-host interaction has been established [5], [49], [51], [52]. Most drug repurposing approaches targeted molecules relevant for viral replication, such as RNA-dependent RNA polymerase. A large number of candidate molecules with potential anti-viral activity has been predicted, but a recent publication has casted doubt on the value of many of these candidate molecules, since drug-induced phospholipidosis rather than targeted antiviral activity may be their non-specific mode of action [53], which may explain why most of these drugs with predicted anti-viral activity failed to show clinical efficacy. Targeting the early steps of SARS-CoV-2 infection also is a promising strategy, i.e., binding of the viral RBD to ACE2 on target cells and the subsequent proteolytic cleavage of the S2 domain of the viral spike protein by the cellular protease TMPRSS2 [3], [54]. Screening of a small molecule library has recently been successful in identification of inhibitors of TMPRSS2 [55], but still no inhibitors acting on the spike 1 part of the RBD have been identified.

Using a structure-based drug-repurposing approaches to find molecules with anti-viral activity against SARS-CoV-2 we identified a molecule with predicted binding to the ACE2-binding site of S1-RBD, which nevertheless failed to show an effect on inhibition of infection, further confirming that structure-based drug repurposing often fails to yield functional anti-viral molecules. However, our approach identified one molecule, montelukast, that prevented SARS-CoV-2 infection, although it was predicted to bind to the S1-RBD outside of the ACE2-binding site at ARG-331, ASN-368, VAL-369, THR-404, SER-488. Montelukast has already been investigated in clinical trials to improve the COVID-19 disease course and was reported to have an effect on the viral proteinase thereby reducing viral replication [56]. Molecular dynamics simulation suggested that montelukast-binding did not disturb the structure of the ACE2-binding site but rather affected the overall S1-RBD structure. Rather than inhibiting the interaction of RBD with its receptor ACE2, these results indicated that montelukast binding to the S1-RBD had an allosteric effect and disturbed the subsequent steps of infection. For SARS-CoV-2 infection to occur, the trimeric spike protein needs to enter into a prefusion state that is mediated through the serine protease function of TMPRSS2 [3], [23], [57], [58], [59], [60]. The increased flexibility of the S1-RBD structure after montelukast binding predicted by molecular dynamics simulation might interfere with TMPRSS2 function or spike protein folding to from the prefusion state and thereby result in reduction of SARS-CoV-2 infection. Since montelukast did not show protease inhibition activity [61], its binding to the S1-RBD may likely interfere with formation of the prefusion state of the spike protein. Interestingly, the inhibitory effect of montelukast on SARS-CoV-2 infection was predicted to be lost in SARS-CoV-2 variants with higher infectivity, i.e., the alpha and the delta variants. Since the effect of montelukast on infection was conserved in the beta variant of SARS-CoV-2, the montelukast binding site appears to be a “sweetspot” to interfere with virus infection. In analogy, antibodies generated against the spike protein after infection and binding to the montelukast binding site may also have infection-inhibition activity. Interestingly, a SARS-CoV-2 infection neutralizing monoclonal antibody was identified [62], that bound at the montelukast-binding region of the S1-RBD. Our finding of a defined binding site in the S1-RBD that interferes with viral infection apparently through an allosteric mechanism will inform further studies to search for molecules with infection inhibition activity across SARS-CoV-2 variants.

CRediT authorship contribution statement

Max Luedemann: Data curation, Invesitagtion, Software, Resources, Formal analysis, Visualization, Writing – original draft, Writing – review & editing. Daniela Stadler: Validation, Invesitagtion, Resources, Formal analysis, Visualization, Writing – review & editing. Cho-Chin Cheng: Validation, Visualization, Formal analysis. Ulrike Protzer: Supervision, Funding acquisition. Writing- review. Percy A. Knolle: Supervision, Funding acquisition, Formal analysis, Writing- draft, Writing – review & editing. Sainitin Donakonda: Conceptualization, Data curation, Methodology, Project administration, Software, Supervision, Writing – original draft, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.