Microsecond molecular dynamics simulations revealed the inhibitory potency of amiloride analogs against SARS-CoV-2 E viroporin

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) encodes small envelope protein (E) that plays a major role in viral assembly, release, pathogenesis, and host inflammation. Previous studies demonstrated that pyrazine ring containing amiloride analogs inhibit this protein in different types of coronavirus including SARS-CoV-1 small envelope protein E (SARS-CoV-1 E). SARS-CoV-1 E has 93.42% sequence identity with SARS-CoV-2 E and shared a conserved domain NS3/small envelope protein (NS3_envE). Amiloride analog hexamethylene amiloride (HMA) can inhibit SARS-CoV-1 E. Therefore, we performed molecular docking and dynamics simulations to explore whether amiloride analogs are effective in inhibiting SARS-CoV-2 E. To do so, SARS-CoV-1 E and SARS-CoV-2 E proteins were taken as receptors while HMA and 3-amino-5-(azepan-1-yl)-N-(diaminomethylidene)-6-pyrimidin-5-ylpyrazine-2-carboxamide (3A5NP2C) were selected as ligands. Molecular docking simulation showed higher binding affinity scores of HMA and 3A5NP2C for SARS-CoV-2 E than SARS-CoV-1 E. Moreover, HMA and 3A5NP2C engaged more amino acids in SARS-CoV-2 E. Molecular dynamics simulation for 1 μs (1,000 ns) revealed that these ligands could alter the native structure of the proteins and their flexibility. Our study suggests that suitable amiloride analogs might yield a prospective drug against coronavirus disease 2019.

SARS-CoV-1 encodes VPs from small envelope protein E, ORF3a and ORF8a genes. Protein E and ORF3a interact with cellular proteins via PDZ-binding motif and exerts IC activity. These functions are essential for optimal viral replication. However, among these 3 VP-forming genes, protein E was indispensable for viral virulence [12]. Protein E forms calcium ion (Ca 2+ ) channels in the endoplasmic reticulum golgi apparatus intermediate compartment (ER-GIC)/Golgi membranes. These ICs alter calcium homeostasis in host cells and activate the cytosolic innate immune signaling receptor NLR family pyrin domain containing 3 inflammasome [13,14]. Absence of protein E attenuates the viral infectious activity by reducing nuclear factor kB mediated inflammation [15]. Therefore, protein E can be a plausible therapeutic target for COVID -19. Previous studies demonstrated that VP-like prototypic M2 proton selective channel of IAV can be an ideal target for antiviral development [16]. All coronavirus E proteins are assumed to form cation-selective ion channels which participate in viral replication and virus's life cycle [17,18]. The actions of VPs in SARS-CoV-2 may also be figured out from other CoVs e.g., Middle East respiratory syndrome coronavirus and human coronavirus 229E (HCoV-229E) [5,6,19]. Deletion of E gene from SARS-CoV-1 significantly decreased the viral pathogenesis [15]. Moreover, the usage of the channel blocking compounds such as amiloride analogs had substantially alleviated the viral replications of HCoV-229E, murine hepatitis virus (MHV), and SARS-CoV-1 [20,21]. Here, the amiloride analogs interacted with the E proteins and inhibited the functions of VP [22,23]. These results have raised a possibility to develop a broad spectrum antiviral drug against coronaviruses [20][21][22][23]. Aside from coronaviruses, Amiloride derivatives (particularly hexamethylene amiloride [HMA]) were found to be efficient inhibitors of ICs in hepatitis C virus, influenza virus, and HIV-1 [7][8][9]. Since E gene of SARS-CoV-1 and SARS-CoV-2 are highly identical (Table 1, Supplementary Data 1), implementation of different in silico tools can unveil the anti-SARS-CoV-2 E mechanisms of amiloride analogs [24,25].

Characterization of SARS-CoV-2 E VP
For characterization of the SARS-CoV-2 E, the E protein (accession: YP_009724392.1) went under Protein Basic Local Alignment Search Tool (BLASTp) in National Center for Biotechnology Information (NCBI) database. To find homologous sequences in other CoVs, SARS-CoV-2 was excluded during this BLASTp. Top 4 hits in BLAST were taken for characterization. Afterward, E protein of SARS-CoV-1 GD01 (accession: AAP51230.1) was also included. The protein sequences were uploaded in ProtParam, Pfam, and THMM to analyze their physicochemical properties, domains, and transmembrane helices [30][31][32]. Their phylogenetic characterization was conducted by CLC Drug Discovery Workbench 3.0 using default parameters (https://digitalinsights.qiagen.com).
Three-dimensional structure generations of E proteins E proteins of SARS-CoV-1 and SARS-CoV-2 were generated via Robetta (https://robetta.bakerlab.org) [33]. The structures were refined with three-dimensional refine and Galaxy Refine [34,35]. The qualities of the generated structures were assessed using SWISS structure assessment [36].

Pharmacokinetic property exploration of the ligands and preparation for molecular docking
The molecules were selected after going through literature [27,28]. The cannonical smiles of 5-(N,N-hexamethylene) amiloride or HMA and 3A5NP2C were uploaded in SWISS-ADME for exploring the pharmacokinetic properties (absorption, distribution, metabolism, excretion) [37]. The canonical SMILES (Simplified molecular-input line-entry system) were collected from PubChem database and transferred into protein data bank file via Online SMILES Translator and Structure File Generator (https://cactus.nci.nih. gov/translate/) [38].

Molecular docking and dynamic simulations
The selected ligands went under blind molecular docking simulation against SARS-CoV-1 and SARS-CoV-2 E proteins via AutoDock tools 1.5.6 [39]. The ligand-receptor interactions between drug-receptor complexes were visualized by BIOVIA Discovery Studio (http://www.discover.3ds.com) and PyMOL Molecular Graphics System, Version 2.3.3 Schrödinger, LLC [40] and UCSF Chimera [41]. The dynamics and stability of SARS-CoV-1 and SARS-CoV-2 E proteins with amiloride analogs bound complexes were compared by carrying out MD simulations. For MD simulations, the drug-bound complexes were embedded in a lipid bilayer membrane composed of 36% phosphatidylcholine, 21% phosphatidylethanolamine, 21% cholesterol, 6% phosphatidylserine. This lipid bilayer mimics the lipid composition of golgi complex [42]. The SARS-CoV-1 VPs generally target golgi apparatus [43]. The total system consisted of about 120,000 atoms in an orthorhombic simulation cell with a free KCl concentration of 150 mM. Equilibrium MD simulations were performed after energy minimization and 1,000 ns of equilibration with position restraints. The protein, water, and lipid components were energy minimized using Charmm-Gui Bilayer Builder [44] and Gromacs 2019.2 [45,46]. All simulations were carried out under periodic boundary condi-tions at constant temperature (T = 310°K) and pressure (P = 1 bar). Then the root mean square deviation (RMSD) and Root Mean Square Fluctuation (RMSF) of C-α carbon from wild type apo-SARS-CoV-1 E and apo-SARS-CoV-2 E with ligand bound complexes were analyzed by g_rms, g_rmsf tools.

SARS-CoV-2 E protein contains relatively higher instability index
According to NCBI BLASTp, SARS-CoV-1 and SARS-CoV-2 E proteins share 93.42% sequence identity. The selected CoV E protein contains nearly same numbers of amino acids outside and inside of the transmembrane. Their grand average of hydropathicity values are not significantly diverse. However, SARS-CoV-2 E has a higher theoretical pI (isoelectric point) and instability index than average value. Only Bat SARS-like coronavirus has this type of properties.
All of the E proteins have Non-structural protein NS3/Small envelope protein E domain (Fig. 1). Some insignificant domains such as Ellis van Creveld protein 2 like protein and FAM163 family were found in SARS-CoV-2 E. This insignificant FAM163 family domain was absent in SARS-CoV-1 E (Supplementary Data 1).

SARS-CoV-2 E protein demonstrated higher binding affinity for the ligands than SARS-CoV-1 E
The E proteins of SARS-CoV-1 and SARS-CoV-2 went under homology modeling and their structures were refined. More than 97% amino acids were in favored region of Ramachandran plot. Details about the structural qualities are given in Supplementary Data 2. To target SARS-CoV-2 E, anticancer amiloride analogs, i.e., HMA and 3A5NP2C were used in this study (Fig. 2). Their physicochemical and pharmacokinetic properties are given in (Table 2) and Supplementary Data 2. They have similar bioavailability. HMA might inhibit CYP1A2 and CYP2C19. However, this characteristic was absent in 3A5NP2C. SARS-CoV-1 E demonstrated -5.5 kcal/mol and -5.8 kcal/mol binding affinity scores with HMA and 3A5NP2C respectively. Whereas HMA and 3A5NP2C interacted with SARS-CoV-2 E with -7.3 kcal/mol and -7.1 kcal/ mol binding affinity scores, respectively ( Table 2).
Only Phe 4 of SARS-CoV-1 E interacted with both HMA and 3A5NP2C. Whereas six common residues of SARS-CoV-2 protein E Phe 4, Asn 66, Leu 12, Tyr 57, Val 62, Lys 63 participated in ligand-receptor interactions (Fig. 3). Moreover, HMA and 3A5N-P2C interacted via H-bonds with SARS-CoV-2 E which was not observed in SARS-CoV-1 E. The pyrazine ring, which is essential for the structure activity relationship of amiloride, showed various types of interactions. The extra pyrimidine group of 3A5NP2C engaged more amino acids in the receptors.

HMA and 3A5NP2C altered the wild type E proteins structures
The ligand-receptor complexes went under 1 ms simulation or 1,000 ns MD simulation. The ligand bound and ligand-free E proteins deviated in same pattern for less than 100 ns and after that their RMSD values were different (Fig. 4). Ligand-receptor complexes also showed altered mobility, especially in the cytoplasmic, transmembrane, and non-cytoplasmic regions. SARS-CoV-2-E-3A5NP2C complex also demonstrated different RMSD value at 100ns whereas SARS-CoV-2-E-HMA showed less deviations. On the other hand, HMA reduced more fluctuation than 3A5NP2C in these domains of SARS-CoV-2 E proteins. During the simulation, the RMSD of SARS-CoV-1 E (black line) and SARS-CoV-1-E-HMA (green line) did not show equal overlapping value after few seconds. This indicated that the bonded drug altered the wild type conformations hence it could not execute the IC activities [22]. Similar patterns were also observed between SARS-CoV-2 E   and drug bound SARS-CoV-2 E. RMSF analysis showed that these compounds reduced the mobility of N and C terminal region of SARS-CoV-2 E (Fig. 4). The schematic representation of the whole study is given in Fig. 5.

Discussion
VPs are involved in viral assembly and pathogenesis, which promotes ion imbalance within host cells, disrupting cellular pathways [11,43,47]. Therefore, impeding their activity by specific drugs of-  fers promising antiviral therapeutics [6,48]. Coronaviruses, including SARS-CoV-1, express small envelope (E) protein that forms VPs. These VPs act as cation-selective ion channels in lipid bilayers that are able to pump Na + and K + ions [49]. The transmembrane α-helical domain of E forming multimeric α-helical bundle is responsible for the IC activity [50]. Similar to SARS-CoV-1, SARS-CoV-2 also encodes VP forming E protein [51,52]. This E protein is 93.42% identical with SARS-CoV-1 E. SARS-CoV-2 E has a conserved NS3/Small envelope protein E domain (NS3_envE) that triggers inflammation in host cells (Fig. 1) [13].
Previous reports showed that NHE-1 blocker HMA inhibits HCoV-229E and MHV E protein IC conductance in lipid bilayers [21]. Functions of SARS-CoV-1 E got significantly reduced by HMA in human embryonic kidney 293 cells [22]. Moreover, HMA and 3A5NP2C have anticancer properties [27,28]. According to SWISS-ADME, these molecules are impenetrable through blood brain barrier and can circulate in the body with 0.55 bioavailability score. 3A5NP2C drug compound has lower gastrointestinal absorption than HMA. However, HMA might act as an inhibitory substrate of CYP1A2 and CYP2C19 enzyme. Specific modification in the side chain of amiloride analogs might make the compound a potential broad spectrum anti-microbial compound with safer ADMET properties and efficacy. Hence, further studies may contribute to the structural modification of these analogs devising a safer and more effective antiviral [53].
In this study, we have found that the structural characteristics of HMA and similar compound 3A5NP2C exhibited more binding affinity score with SARS-CoV-2 E than SARS-CoV-1 E. Phe 4 of SARS-CoV-1 E was the common interacted residue for both HMA and 3A5NP2C whereas SARS-CoV-2 E interacted via Phe 4, Asn 66, Leu 12, Tyr 57, Val 62, Lys 63. This indicated that these ligands could bind to both inside (amino acid number: 1-11), outside (amino acid number: 35-75) and within transmembrane helix (amino acid number: 12-34) regions of SARS-CoV-2 E ( Table 1).
The E proteins and drug-receptor complexes were placed in a Golgi lipid bilayer model for MD simulations. SARS-CoV-2 E showed higher mobility and deviations in bilayer than SARS-CoV-1 E. Especially the outside N and C terminal regions of the protein were more mobile than SARS-CoV-1 E. Since, SARS-CoV-2 can replicate 3 times faster than SARS-CoV-1; hence, it can release/shed more viruses from the host cells [54]. Whether this enhanced mobility and instability of SARS-CoV-2 E are contributing in higher viral release is further needed to be explored (Table 1, Fig. 4). Here, the SARS-CoV-1 E and SARS-CoV-1-E-HMA complex were run as control since HMA can inhibit the VP activities. The differences between RMSD and RMSF of the ligand-bound and ligand-free SARS-CoV-1 E indicated that these compounds interfered with the normal physiology of the viral protein. The same patterns were also observed for ligand-bound and ligand-free SARS-CoV-2 E proteins. The flexibility of the proteins was altered in cytoplasmic, non-cytoplasmic and transmembrane regions. These flexibilities are critical for viral physiology [43]. MD simulation revealed that this region would gain rigidity which strongly supports that the viral physiology will be hindered extremely [55]. Hence, these deviations and alterations of the protein E strongly suggest that these amiloride analogs will disrupt the cation-selective IC activities.
Amiloride and HMA are well-studied compounds; therefore, in the near future, their antiviral activity against SARS-CoV-2 can be evaluated to develop new treatments. However, amiloride has some rare adverse effects and side effects. The most dangerous effects include hyperkalaemia [56]. Therefore, similar to nafamostat mesylate, serum potassium values should be monitored carefully in COVID-19 patients after the administration of amiloride analogs [57]. Although HMA has a lower K + sparing diuretic effect than amiloride, these drugs might show side effects by hindering the general physiology of uPA [27,28].
Several amiloride analogs can inhibit the functions of coronavirus VPs. Our study strongly suggests their antiviral activities against SARS-CoV-2. Further in vitro screening and in vivo experiments are necessary to consider amiloride analogs as a prospective drug against this virus.

Conflicts of Interest
No potential conflict of interest relevant to this article was reported.