In silico identification of Tretinoin as a SARS-CoV-2 envelope (E) protein ion channel inhibitor.

Viroporins are oligomeric, pore forming, viral proteins that play critical roles in the life cycle of pathogenic viruses. Viroporins like HIV-1 Vpu, Alphavirus 6 K, Influenza M2, HCV p7, and Picornavirus 2B, form discrete aqueous passageways which mediate ion and small molecule transport in infected cells. The alterations in host membrane structures induced by viroporins is essential for key steps in the virus life cycle like entry, replication and egress. Any disruption in viroporin functionality severely compromises viral pathogenesis. The envelope (E) protein encoded by coronaviruses is a viroporin with ion channel activity and has been shown to be crucial for the assembly and pathophysiology of coronaviruses. We used a combination of virtual database screening, molecular docking, all-atom molecular dynamics simulation and MM-PBSA analysis to test four FDA approved drugs - Tretinoin, Mefenamic Acid, Ondansetron and Artemether - as potential inhibitors of ion channels formed by SARS-CoV-2 E protein. Interaction and binding energy analysis showed that electrostatic interactions and polar solvation energy were the major driving forces for binding of the drugs, with Tretinoin being the most promising inhibitor. Tretinoin bound within the lumen of the channel formed by E protein, which is lined by hydrophobic residues like Phe, Val and Ala, indicating its potential for blocking the channel and inhibiting the viroporin functionality of E. In control simulations, tretinoin demonstrated a lower binding energy with a known target as compared to SARS-CoV-2 E protein. This work thus highlights the possibility of exploring Tretinoin as a potential SARS-CoV-2 E protein ion channel blocker and virus assembly inhibitor, which could be an important therapeutic strategy in the treatment for coronaviruses.

Introduction

Although viruses from different families have distinct life cycles; there exist striking mechanistic similarities in host interaction pathways. A few such commonalities include host membrane fusion by enveloped viruses, membrane disruption by non-enveloped viruses orchestrated by amphipathic viral peptides, and viroporin mediated membrane alteration/remodelling for facilitating virus replication and egress. Both enveloped and non-enveloped viruses utilize virally encoded membrane proteins or “viroporins” at different stages of their life cycles; and membrane remodelling by viroporins appears to be a crucial step during enveloped virus assembly and budding [1]. Examples of well-studied viroporins include M2 protein of Influenza A Virus, 6 K proteiene">n of n class="Species">Chikungunya Virus (CHIKV), Vpu of HIV-1 and p7 of HCV [1,2]. Small molecule inhibitors block multiple viroporins from different virus families; in spite of these components lacking in sequential or structural similarities [2]. This shows a degree of mechanistic equivalency between viroporins, indicating that they may have similar roles to execute in the viral life cycle [2,3].

The “E” or “Envelope” protein of coronaviruses is a viroporiene">n produced duriene">ng n class="Disease">viral infection, although only a small percentage is incorporated into assembled virions [4]. The majority of the protein remains associated with ER, golgi and the ERGIC pathway, which are the sites for coronavirus assembly and release [5]. E is a small membrane protein of 76–109 amino acids, with a short N-terminal region of 7–12 amino acids, a central transmembrane domain (TMD) and a cytoplasmic tail, which is largely hydrophilic [4,6]. The major function of E in coronavirus biology appears to be to support assembly and release of particles [4,7]. It has been shown by several groups that association between the cytoplasmic tails of E and M (Membrane protein) of coronaviruses is essential for particle formation [8]. While M interacts with multiple structural proteins, the involvement of E is crucial for membrane bending and scission during virus budding, leading to the “pinching off” of the virus from ERGIC [4,7,8]. Removal of E from the structural cassette of SARS-CoV has a devastating effect on virus assembly and replication, resulting in the production of immature or incompetent progeny [9].

The transmembrane domain of E is essential for its viroporin activity [3,6]. A synthetic version of the TMD assembles iene">nto peene">ntamers that form cation-selective membrane pores. The ion channels geene">nerated by E of SARS-CoV in ERGIC membranes transports Ca2+, which was shown to activate the NLRP3 inflammasome activation, and IL-1β overproduction [10]. This phenomenon was linked to the worsening of lung damage observed in SARS-CoV infected mice, and thus to the immunopathology associated with viral infection [10]. The ion-channel activity of E is crucial for the virus life cycle. Point mutations N15A and V25F introduced in SARS-CoV E causes disruption of its ion-channel activity [4,6,7], which, in turn, results in the accumulation of compensatory mutations in the protein that restore its membrane activity. However, the mechanistic link between the ion channel activity of E, and its membrane bending functionality that facilitates virus budding, is not clearly understood.

The E protein of coronaviruses contaiene">ns several motifs, apart from the n class="Disease">TMD, that are necessary for its functionality. A β-coil-β-motif, containing a highly conserved proline residue, is necessary for the targeting of the protein to the golgi complex. Mutation of this residue retargets E to the plasma membrane, which inhibits the assembly and release of CoV particles [11]. A PDZ-binding motif (PBM) located in the terminal four amino acids of E is involved in several, crucial host factor interactions that appear to contribute to the pathogenicity of the virus [12,13]. Mutations in the PBM also rapidly revert in cell culture. Some of the reported host partners of the E-protein, which interact through the PBM, include the B-cell lymphoma extra-large protein (Bcl-xL), which has anti-apoptotic properties; the protein associated with Caenorhabditis elegans lin-7 protein 1 (PALS1); syntenin, sodium/potassium (Na+/K+) ATPase α-1 subunit, and stomatin [4,13] Some of these interactions have been shown to contribute significantly to the pathogenicity of CoVs. For example, the interaction of E protein with PALSI has been shown to disrupt tight junctions in the lungs, allowing virus particles to cross the alveolar barrier [14]. Likewise, interaction of E with syntenin results in overexpression of inflammatory cytokines that is thought to contribute to the tissue damage caused by coronaviruses [4,13].

Given the importance of E in virus propagation and new virus assembly and budding, this protein can be considered as a crucial drug target for antiviral generation. Previous studies with viroporins from different virus families have reported amantadine as a broadly functional iene">nhibitor of ion channel activity [2,6,15,16]. Also, a receene">nt study has reported Gliclazide and Memantine as potential inhibitors of E protein channel activity [17].

In this work, we targeted the E protein of the novel coronavirus n class="Species">SARS-CoV-2 (nCoV19) for drug repurposing studies using computational techniques. Utilizing a combination of in silico molecular docking, 200ns all-atom molecular dynamics simulations, H-bonding and binding energy analysis (MM-PBSA), we tested four specific FDA-approved drugs, Tretinoin, Mefenamic acid, Ondansetron and Artemether were evaluated for their ability to bind SARS-2 E protein. Out of the four drugs Tretinoin was selected as the best candidate due to its ability to form extensive H-bonding interactions and high binding energy value (− 412.8 kJ/mol). It is hoped that blocking of the ion-channels formed by the viroporin E can have a substantial detrimental effect on SARS-CoV-2 assembly and propagation.

Materials and methods

In silico 3D structure generation and homo-oligomerization of SARS-CoV-2 E protein

The structure for SARS-CoV-2 E proteiene">n monomer was geene">nerated usiene">ng the I-TASSER (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) web server [18,19]. The E proteiene">n sequeene">nce utilized for structure prediction was from the Wuhan-1 isolate (Geene">nBank ID: QHD43418.1). The template used for structure prediction was the n class="Gene">NMR structure of the E protein from SARS CoV (PDB ID: 5X29)[20]. The model with the best C score (−0.75) was chosen for further processing. The predicted model was then refined using GalaxyRefine server [21]. The final selected 3D structure had a clash score of 1.6 with 94.5% of residues in the favored region. The secondary structure content of the predicted E protein monomer had a 66.7% alpha-helical content. Since the E protein from SARS-CoV has been shown to form pentamers [6], the predicted 3D structure of the SARS-CoV-2 E protein monomer was subjected to homo-oligomerization using the GalaxyHomomer server [22] with the SARS CoV E protein pentamer as the template. The resultant pentameric arrangement was energy minimized prior to being used for docking and simulation studies.

Virtual screening and molecular docking

All ligand structures were obtained from the DrugBank database (https://www.drugbank.ca/) [23]. For ligand screening, the DrugScreen server (http://cao.labshare.cn/drugscreen/) was utilized. Ligand search was conducted on a database of 1806 FDA-approved drug molecules, which generated a list of 50 top ranking molecules as output. Further screening was carried out based on the chemical nature, therapeutic target, side-effects and the status of FDA approval of molecules, which gave a total of 15 candidates [Supplemeene">ntary Table T1]. Out of these four compounds were further tested.

The pentameric arrangement of SARS-CoV-2 E was docked with pre-selected compounds usiene">ng the dockiene">ng server CB-Dock (http://cao.labshare.cene">n/cb-dock/)[24]. CB-Dock performs bliene">nd dockiene">ng usiene">ng automated cavity-detection[24], followed by geene">neration of docked poses usiene">ng AutoDock Viene">na [25]. Selected ligand candidates were subjected to geometry optimization usiene">ng the Uene">niversal Force Field [26,27]; and eene">nergy miene">nimization with the steepest desceene">nt algorithm [28,29]usiene">ng the Auto Optimize tool of Avogadro [30] prior to dockiene">ng. n class="Chemical">Hydrogens and charges were added using the DockPrep tool of UCSF Chimera [31]. Binding affinity of docked ligands were evaluated using the PRODIGY web server [32]. Protein-ligand interactions of docked complexes were analyzed using Protein Ligand Interaction Profiler (PLIP) server (https://projects.biotec.tu-dresden.de/plip-web/plip/index) [33]. Visualization of protein-ligand interactions was carried out using UCSF Chimera [31] and PyMOL [34,35].

MD simulation

In order to understand the binding dynamics and stability of the protein-ligand complexes, multiple all-atom molecular dynamics simulations were carried out using GROMACS (v5.1.1) and the Gromos force field (43a1). Topology and parameters for each ligand was calculated usiene">ng the PRODRG server [36]. A system consistiene">ng of E proteiene">n assembly, either with or without docked ligands, was embedded withiene">n a pre-equilibrated 392 POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membrane. System solvation was carried out using the Simple Point Charge (SPC) water model [37] followed by ionization with chloride (Cl−) ions for system charge neutralization. 50,000 steps of energy minimization was performed using steepest descent algorithm [28,29]with a tolerance of 1000 kJ mol−1 nm−1. This was followed by 500 ps of NVT (isothermal-isochoric ensemble) and 5ns of NPT (isothermal-isobaric ensemble) with position-restraints. The system temperature was maintained at 300 K using Noose-Hoover thermostat [38] and pressure was kept at 1 bar with Parrinello-Rahman pressure-coupling [39] with a compressibility of 4.5 × 10−5 bar−1 and 2ps time constant. Simulation timestep was kept at 2fs.Bond lengths were constrained using the LINCS (Linear Constraint Solver) algorithm [40]. Electrostatics were treated with PME (Particle mesh Ewald) [41]. A 1.2 nm cut-off was used for both coulombic and van der Waals's interactions. Post equilibration, unrestrained production run was carried out for 200 ns. Trajectory analysis was carried out using Gromacs analysis tools and VMD [42].

Additionally, two control simulations were also carried out with the known targets of Tretinoin (n class="Chemical">DB00755) and Mefenamic acid (DB00784) which are Human Retinoic Acid Receptor Gamma and Cyclooxygenase-2 respectively. Briefly, the crystal structures of the complexes (2LBD and 5IKR) were subjected to 200 ns all-atom MD simulation under explicit solvent conditions at 300 K and 1 bar atmospheric pressure. Triplicate MD runs with random initial distributions of the atomic velocities were performed amounting to total simulation time of 600 ns for each system.

MM-PBSA analysis

MM-PBSA stands for Molecular Mechanics Poisson-Boltzmann Surface Area. It is a method of choice for calculatiene">ng the Gibb's free eene">nergy of biene">ndiene">ng iene">n proteiene">n-ligand and proteiene">n-proteiene">n simulation studies [[43], [44], [45]]. It is considered an iene">ntermediate iene">n terms of accuracy and computational cost betweeene">n empirical scoriene">ng and alchemical perturbation techene">niques [46]. g_mmpbsa was utilized for quantifying the Gibb's free energy of binding of the protein-ligand complexes from the generated trajectories. It uses the following equations for its calculations.ΔG ΔG ΔE ΔGwhere.

ΔG = total free energy of protein-ligand complex

ΔG , ΔG = free energies of protein and ligand in solvent respectively

ΔG  = free energy for each entity

ΔE  = vacuum potential energy (total energy of bonded and nonbonded interactions)

ΔG  = solvation free energy including polar (ΔG ) and nonpolar (ΔG ) energies

TΔS = entropic contribution of free energy in vacuum. T and ΔS is the temperature and entropy respectively.

Since the g_mmpbsa tool does not calculate the eene">ntropic term (S), the term TΔS MM was not iene">ncluded iene">n this study. Thus, the ΔG biene">nd was calculated as follows, where ΔGSASA is equivalent to ΔGnonpolar since the calculation of nonpolar solvation energy was based on the SASA model in this study, ΔG

Comprehensive information regarding the concept and the protocol for use in Gromacs is described [46]. In addition to determining overall free energy of binding between the protein and ligand, residue level contribution towards the overall binding energy was also determined for in-depth insight on the mechanism of binding. The following formula was used,where Ai bound and Ai free represent energies of ith atom from xth residue with and without ligand whereas n is total number of atoms which makes up the residue.

Results

Novel coronavirus envelope protein is a membrane spanning viroporin

The envelope or E protein of Coronaviruses is a siene">ngle transmembrane domaiene">n contaiene">niene">ng proteiene">n with hydrophilic termiene">ni. Sequeene">nce based analysis of n class="Species">SARS-CoV-2 E protein using hydrophobicity plots [Fig. 1 A] and multiple transmembrane domain prediction servers [Fig. 1B] suggested the presence of a conserved, prominent, central transmembrane domain between residues 12–39. Sequence composition of the predicted TMD (12–39 aa) of the SARS-CoV-2 E protein indicated the presence of multiple hydrophobic residues like Alanine, Valine, Phenylalanine and Leucine [Fig. 1C]. The N and C termini were hydrophilic, as expected [Fig. 1A]. Both the hydrophobicity profile and the sequence of E were found to be highly conserved between SARS-CoV and SARS-CoV-2 [Fig. 1A and D]. The SARS CoV E protein structure was earlier solved by solution NMR [6]. In the absence of available 3D structures for the SARS-CoV-2 E protein, in silico structural prediction was carried out using I-TASSER [Fig. 2 A], with the NMR structure of SARS CoV E as the template. The predicted 3D structure of SARS-CoV-2 E protein was dominated by the central helical transmembrane domain with N and C-terminus oriented opposite to each other [Fig. 2A]. Ramachandran plot analysis revealed that 94.5% of the residues are in the favored region with a total of two outliers [Fig. 2C]. The pentameric assembly [Fig. 2B and D] showed the most likely transmembrane channel [47]lined with residues Val14, 25, 29, Leu18, 21, 34, Ala22, Phe26, Thr30,Ile33 and Leu37 [Fig. 2B and D]contributed by each monomer. Channel lining residues were identified using PoreWalker server [48]. The narrowest region of the pentameric channel, which had a diameter of 6.5 Å, was lined by Phe26contributed by each monomer. This region of the channel formed a constricted hydrophobic pocket [Fig. 2D].

Fig. 1

(A) Kyte-Doolittle hydrophobicity plot of SARS CoV and SARS-CoV-2 E protein sequences. (B) Table showing the predicted transmembrane domain of SARS-CoV-2 E protein by three different transmembrane domain prediction servers [[49], [50], [51]] (C) Helical wheel plot of the transmembrane domain of SARS-CoV-2 E protein showing the amino acid composition, (D) Clustal Omega alignment of the full-length SARS CoV SARS-CoV-2 E protein sequences.

Fig. 2

(A) Predicted structure of the SARS-CoV-2 E protein monomer in ribbon representation (purple: alpha-helix; black: coil), (B) Top view of the pentameric arrangement of SARS-CoV-2 E protein, (C) Ramachandran plot of the predicted SARS-CoV-2 E protein monomer, (D) Line representation of the pentameric SARS-CoV-2 E protein showing the predicted ion channel. Inset shows the positioning of the Phe26 residues (red) at the narrowest point of the ion channel.

(A) Kyte-Doolittle hydrophobicity plot of SARS CoV and n class="Species">SARS-CoV-2 E protein sequences. (B) Table showing the predicted transmembrane domain of SARS-CoV-2 E protein by three different transmembrane domain prediction servers [[49], [50], [51]] (C) Helical wheel plot of the transmembrane domain of SARS-CoV-2 E protein showing the amino acid composition, (D) Clustal Omega alignment of the full-length SARS CoV SARS-CoV-2 E protein sequences.

(A) Predicted structure of the SARS-CoV-2 E proteiene">n monomer iene">n ribbon represeene">ntation (purple: alpha-helix; black: coil), (B) Top view of the peene">ntameric arrangemeene">nt of n class="Species">SARS-CoV-2 E protein, (C) Ramachandran plot of the predicted SARS-CoV-2 E protein monomer, (D) Line representation of the pentameric SARS-CoV-2 E protein showing the predicted ion channel. Inset shows the positioning of the Phe26 residues (red) at the narrowest point of the ion channel.

Docking studies highlight important residues involved in ligand binding

Blind molecular docking was carried out with the 15 selected drugs from the DrugBank database using the CB-DOCK web server. Based on the number of binding modes within the channel lumen- Tretinoin (n class="Chemical">DB00755), Mefenamic acid (DB00784), Ondansetron (DB00904) and Artemether (DB06697) were selected [Table 1 , Supplementary Fig. S1] for further analysis. Binding affinity calculation of docked complexes [Fig. 3 E] showed that Tretinoin possessed the highest binding affinity of −v8.6 kcal/mol followed by Ondansetron (−8.5 kcal/mol), Mefenamic Acid (−7.9 kcal/mol) and Artemether (−7.0 kcal/mol). Three of the four ligands -Tretinoin, Mefenamic acid and Ondansetron -bound to the narrowest region of the channel, lined by Phe26 residues [Fig. 2D]. Artemether, on the other hand, bound towards the N terminal region [Fig. 2D]. In all cases, the ligands interacted via hydrophobic interactions with the protein residues from all chains lining the channel [Fig. 3A–E]. Tretinoin formed a total of 11 interactions involving the residues Leu18 (Chain C, D), Leu21 (Chain C), Ala22 (Chain C), Val25 (Chain C) and Phe26 (Chains A-E). Ondansetron was observed to form 7 hydrophobic interactions with Ala22 (Chain B, C, D), Val25 (Chain B) and Phe26 (Chain B, C). Mefenamic acid exhibited 6 interactions with residues Leu18 (Chain B, C), Leu21 (Chain D) and Ala22 (Chain B, C) [Fig. 3A–E]. Lastly, Artemether showed interactions with Phe4 (Chains A-E) only [Fig. 3A–E].

Table 1

Details of selected drugs for the study from DrugBank.

S No.	DrugBank ID	Approval Status	MW (Average)	Approved Usage
1	DB00755	Approved, Investigational, Nutraceutical	300.4	Acne treatment, photodamaged skin, keratinization disorders and Acute Promyelocytic Leukemia.
2	DB00784	Approved	241.3	Anti-inflammatory agent
3	DB00904	Approved	293.4	Prevention of nausea and vomiting associated with emetogenic cancer chemotherapy
4	DB06697	Approved	298.4	Anti-malarial agent

Fig. 3

(A–D) Relative positioning of ligands (A: Tretinoin-DB00755; B: Mefenamic acid-DB00784; C: Ondansetron-DB00904 and D: Aremether-DB06697) with respect to SARS-CoV-2 E. Protein is represented as ribbons (Monomer 1–5: yellow, magenta, cyan, green and beige) and ligands are represented as ball and stick. (E) PRODIGY server predicted binding affinity values in ΔG (Kcal/mol), for each tested compound and their interacting partner.

(A–D) Relative positioning of ligands (A: Tretinoin-n class="Chemical">DB00755; B: Mefenamic acid-DB00784; C: Ondansetron-DB00904 and D: Aremether-DB06697) with respect to SARS-CoV-2 E. Protein is represented as ribbons (Monomer 1–5: yellow, magenta, cyan, green and beige) and ligands are represented as ball and stick. (E) PRODIGY server predicted binding affinity values in ΔG (Kcal/mol), for each tested compound and their interacting partner.

Molecular dynamics studies

All-atom molecular dynamics simulation was carried out with the SARS-CoV-2 E proteiene">n peene">ntamer and E proteiene">n-ligand complexes embedded withiene">n a pre-equilibrated POPC membrane for a total of 200 ns. RMSD (Root Mean Square Deviation) analysis of the C-alpha backbones of the E protein monomers showed that the systems were stable during the simulation time period [Fig. 4 A-D], Supplementary Fig. S2]. All the four ligands also showed a relatively stable RMSD trace except Artemether which showed relatively large fluctuations during simulation [Fig. 4E-H]. Radius of gyration values of SARS CoV-2 E showed a steady state in presence of all four ligands which indicates that the protein maintained a relatively compact state during the simulation time period [Fig. 5 A-D].

Fig. 4

(A–D) Shows the RMSD plots of the C-alpha backbone of SARS CoV-2 E protein in presence of complexed ligands. (A) SARS CoV-2 E-DB00755, (B) SARS CoV-2 E-DB00784, (C) SARS CoV-2 E-DB00904, and (D) SARS CoV-2 E-DB06697 respectively (triplicate runs). (E–H) shows the RMSD trace of ligands in complex with SARS CoV-2 E protein where, (E) DB00755 (Tretinoin, green), (F) DB00784 (Mefenamic acid, purple), (G) DB00904 (Ondansetron, grey) and (H) DB06697 (Artemether, blue) respectively (). For (A–D) color coding is run1 (blue), run2 (red), run3 (yellow).

Fig. 5

(A–D) Shows the radius of gyration plots of SARS CoV-2 E protein in presence of complexed ligands. (A) SARS CoV-2 E-DB00755, (B) SARS CoV-2 E-DB00784, (C) SARS CoV-2 E-DB00904, and (D) SARS CoV-2 E-DB06697 respectively (triplicate runs). Color coding is run1 (blue), run2 (red), run3 (yellow).

(A–D) Shows the RMSD plots of the C-alpha backbone of n class="Species">SARS CoV-2 E protein in presence of complexed ligands. (A) SARS CoV-2 E-DB00755, (B) SARS CoV-2 E-DB00784, (C) SARS CoV-2 E-DB00904, and (D) SARS CoV-2 E-DB06697 respectively (triplicate runs). (E–H) shows the RMSD trace of ligands in complex with SARS CoV-2 E protein where, (E) DB00755 (Tretinoin, green), (F) DB00784 (Mefenamic acid, purple), (G) DB00904 (Ondansetron, grey) and (H) DB06697 (Artemether, blue) respectively (). For (A–D) color coding is run1 (blue), run2 (red), run3 (yellow).

(A–D) Shows the radius of gyration plots of SARS CoV-2 E proteiene">n iene">n preseene">nce of complexed ligands. (A) n class="Species">SARS CoV-2 E-DB00755, (B) SARS CoV-2 E-DB00784, (C) SARS CoV-2 E-DB00904, and (D) SARS CoV-2 E-DB06697 respectively (triplicate runs). Color coding is run1 (blue), run2 (red), run3 (yellow).

Out of four compounds, only Tretinoin and n class="Chemical">Mefenamic acid formed stable hydrogen bonds [Table 2 ] with E protein residues within the channel lumen, apart from hydrophobic interactions. Tretinoin formed multiple high occupancy H-bonds (≥10%) with Phe20, Leu21 Ala22 and Phe23, with the strongest interaction being with Phe23 that exhibited an occupancy value of 58.9%. Mefenamic acid on the other hand, formed H-bonds with Ala22, Val24 and Val25, with the strongest interaction with Val24having an occupancy value of 44.9%. Quantification of the total number of H-bonds formed during the simulation time period also showed a steady decrease in the number and frequency of H-bonds in the following order (Tretinoin > Mefenamic Acid > Artemether > Ondansetron) [Fig. 6 A–D]. In terms of hydrophobic interactions, Tretinoin and Mefenamic acid formed multiple hydrophobic interactions with Val17, Leu18, Leu21, Ala22, Phe23, Val25 and Phe26. However, unlike in case of H-bonds, multiple hydrophobic interactions were observed for Ondansetron and Artemether with residues like Phe4, Glu8, Ile13, Leu21, Ala22 and Phe26.

Table 2

Hydrogen bond occupancy values for SARS-CoV-2 E − drug complexes.

DrugBank ID	Protein-ligand Hydrogen bond	Percentage Occupancy (≥10%)
DB00755	Phe23	58.9
	Phe20	35.6
	Leu21	25.0
	Ala22	13.9
DB00784	Val24	44.9
	Val25	39.5
	Ala22	10.6

Fig. 6

Graphs showing number of hydrogen bonds formed during the simulation period for the ligands - (A) Tretinoin (DB00755), (B) Mefenamic acid (DB00784), (C) Ondansetron (DB00904) and (D) Artemether (DB06697).

Graphs showing number of hydrogen bonds formed duriene">ng the simulation period for the ligands - (A) Tretinoin (DB00755), (B) Mefenamic acid (DB00784), (C) Ondansetron (DB00904) and (D) Artemether (DB06697).

To further evaluate the mechanistics of biene">ndiene">ng of the ligands to the SARS-CoV-2 E protein pentamer during simulation, MM-PBSA analysis was carried out.

It was observed that for Tretinoin and n class="Chemical">Mefenamic acid, electrostatic energy and polar solvation energies were the major contributors towards the binding energy; whereas, for binding of Ondansetron and Artemether to the channel, van der Waals's energy played the major role [Table 3 , Supplementary Table 3]. Out of the four compounds tested, Tretinoin showed the maximum binding energy value of −412.8 kJ/mol followed by Mefenamic acid (−381.4 kJ/mol), Ondansetron (−227.1 kJ/mol) and lastly Artemether (−173.5 kJ/mol) respectively. At the residue level, the major contributors for binding to Tretinoin and Mefenamic acid were Met1, Leu19, Phe20, Leu21, Ala22, Phe23, Val24, Val25, Phe26, Arg38, Lys53, Arg61, Lys63 and Arg69. For Ondansetron and Artemether, Phe4, Ile13 and Leu18 were majorly involved in binding [Fig. 7 A-D].

Table 3

Binding free energy (MM-PBSA) analysis of SARS-CoV-2 E-drug complexes.

DrugBank ID	Binding Energy (kJ/mol)	Van der waal's Energy (kJ/mol)	Electrostatic Energy (kJ/mol)	Polar Solvation Energy (kJ/mol)	SASA Energy (kJ/mol)
DB00755	−412.8 ± 12.8	−185.2 ± 13.4	−502.7 ± 14.6	293.0 ± 7.3	−17.9 ± 0.9
DB00784	−381.4 ± 11.9	−191.7 ± 11.1	−424.9 ± 13.2	249.7 ± 6.8	−14.4 ± 0.6
DB00904	−227.1 ± 9.9	−232.5 ± 9.5	−6.8 ± 2.5	28.6 ± 2.2	−16.4 ± 0.7
DB06697	−173.5 ± 9.5	−187.2 ± 8.7	−2.4 ± 0.8	31.5 ± 3.0	−15.4 ± 0.8

Fig. 7

Binding energy contribution per residue of the binding sites corresponding to the ligands - (A) Tretinoin (DB00755), (B) Mefenamic acid (DB00784), (C) Ondansetron (DB00904) and (D) Artemether (DB06697).

Discussion

The current pandemic caused by the novel coronavirus n class="Species">SARS-CoV-2 has resulted in 4.5 million infections and 308,000 deaths so far. To combat this unprecedented public health emergency, developing methods of prevention and cure is crucial. Besides designing a safe and efficacious vaccine, effective treatment options are absolutely essential for already infected patients or identified asymptomatic carriers. Thus, understanding common host-virus interactions, and developing inhibitors targeted towards these processes, appears to be a viable strategy. Structure-based drug design and drug repurposing efforts against protein components of SARS-CoV-2 have primarily focused on the main protease (3CLpro), the Spike protein (S) and the viral replicase [[52], [53], [54], [55], [56], [57], [58], [59]], because of their prominent roles in the life cycle of the virus. However, existing data on other coronaviruses indicate that the small E protein has an outsize role in the life cycle of these viruses [9]. The functional similarity of E with viroporins from other virus families, and the previous success in identifying inhibitors for these components [2,6,15,16], indicates that similar strategies may been utilized to identify or repurpose inhibitors against E of SARS-CoV-2.

Recent studies have highlighted the potential of targeting SARS-CoV-2 E proteiene">n as a therapeutic strategy. Compounds like n class="Chemical">Rutin, Doxycycline [60], Nimbolin A [61], Belachinal, Macaflavanone E, Vibsanol B [62] have been proposed as potential inhibitors of the SARS-2 E protein. However, their experimental validation remains. We initiated a thorough search for potential inhibitors of the SARS-CoV-2 E protein by screening 1806 FDA-approved drugs using the DrugScreen web server. The top 50 hits were further screened based on specific attributes of the drugs, resulting in a final list of 15 candidates [Supplementary Table T1]. Four of these candidates (Tretinoin, DB00755; Mefenamic Acid, DB00784; Ondansetron, DB00904 and Artemether, DB06697) were tested for their suitability as inhibitors of SARS-CoV-2 E protein by docking and simulation studies. Ligand binding analysis of docked complexes showed a predominance of hydrophobic interactions in all cases [Fig. 3A–E]. Except Artemether, all other ligands occupied a region in the pentameric channel formed by the E protein. This channel includes protein residues 18–26 and is characterized by the predominance of hydrophobic amino acids like Leucine, Phenylalanine and Valine. The region forms a constriction in the channel diameter, in the form of a hydrophobic cavity. Interestingly, in two studies [[61], [62]] residues Val25 and Phe26 of the SARS-2 E protein were found to play an important role in binding to potential inhibitors. Our study validates the importance of these residues in forming hydrophobic interactions and intermolecular H-bonds with potential inhibitors. This indicates that the hydrophobic stretch encompassing residues 18–26 might be the major ligand-binding site in SARS-CoV-2 E protein for a variety of small molecules. The blocking of the SARS-CoV-2 E ion channels by small molecules has the potential to strongly inhibit the viroporin activity of E, and consequently nullify its contribution to virus assembly. A previous NMR study on SARS-CoV E protein had demonstrated binding of an inhibitor Hexamethylene Ameloride (HMA) mediated primarily by Asn15 [6]. Based on recent literature and our current work, the mode of inhibitor binding appears to be different for the SARS-CoV-2 E protein channel, in spite of strong sequence similarities between the E proteins. This variation might be attributed to the differences in the molecular architecture of the tested compounds.

In our study, Tretinoin formed the strongest hydrophobic iene">nteractions withiene">n this cavity as en class="Chemical">valuated by binding affinities [Fig. 3E]. Control simulations of Tretinoin and Mefenamic acid in complex with their respective known targets (Human Retinoic Acid Receptor Gamma (PDB ID: 2LBD) and Human Cyclooxygenase-2 (PDB ID: 5IKR) respectively) showed comparatively lower binding energy values [Supplementary Table 2], which indicates that there is a strong possibility for these compounds to bind to SARS-2 E protein preferentially, and with higher affinity. The complexes were then subjected to 200 ns of explicit all-atom MD simulations in presence of a lipid membrane. The presence of the membrane was necessary to provide a native-like environment for the pentameric ion channel. Evaluation of protein drug interactions during the simulation time period revealed formation of high occupancy hydrogen bonds in case of Tretinoin and Mefenamic acid apart from hydrophobic interactions [Table 2]. However, Ondansetron and Artemether established H-bonds sporadically [Fig. 5]. This indicated that Tretinoin and Mefenamic acid are able to interact more extensively with the channel. MM-PBSA analysis also showed that Tretinoin exhibited the strongest binding energy value among the four drugs tested. It is hoped that the role of proposed E-protein channel inhibitors on SARS-CoV-2 assembly can be biochemically validated in future.

Conclusion

In this work, we have utilized in silico-based techniques of molecular docking and MD simulations to identify potential inhibitors against the Envelope (E) protein of SARS-CoV-2, which is a relatively less explored but very critical component of the virus. Based on docking and simulation studies, Tretinoin emerged as the best candidate with the ability to inhibit the viroporin activity of SARS-CoV-2 E and affect virus assembly; and should be further tested in vitro and in vivo. The other three FDA approved drugs could also be evaluated against SARS-CoV-2 for their ability to inhibit the virus propagation by disrupting the ion channel functionality of E.

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.