Molecular dynamic simulations reveal anti-SARS-CoV-2 activity of mitocurcumin by potentially blocking innate immune evasion proteins NSP3 and NSP16

The coronavirus disease 19 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is affecting human life in an unprecedented manner and has become a global public health emergency. Identification of novel inhibitors of viral infection/replication is the utmost priority to curtail COVID-19 progression. A pre-requisite for such inhibitors is good bioavailability, non-toxicity and serum stability. Computational studies have shown that curcumin can be a candidate inhibitor of certain SARS-CoV-2 proteins; however, poor bio-availability of curcumin limits its possible therapeutic application. To circumvent this limitation, we have used mitocurcumin (MC), a triphenyl phosphonium conjugated curcumin derivative, to study the ability to inhibit SARS-CoV-2 infection using molecular docking and molecular dynamics (MD) simulation. MC is serum stable and several fold more potent as compared to curcumin. Molecular docking studies revealed that MC can bind at active site of SARS-CoV-2 ADP Ribose Phosphatase (NSP3) and SARS-CoV-2 methyltransferase (NSP10-NSP16 complex) with a high binding energy of − 10.3 kcal/mol and − 10.4 kcal/mol, respectively. MD simulation (100 ns) studies revealed that binding of MC to NSP3 and NSP16 resulted in a stable complex. MC interacted with critical residues of NSP3 macro-domain and NSP10-NSP16 complex and occupied their active sites. NSP3 is known to suppress host immune responses whereas NSP10-NSP16 complex is known to prevent immune recognition of viral mRNA. Our study suggests that MC can potentially inhibit the activity of NSP3 and NSP10-NSP16 complex, resulting in compromised viral immune evasion mechanism, and thereby accentuate the innate immune mediated clearance of viral load. Graphical abstract


Introduction
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel strain of Severe Acute Respiratory Syndrome-related coronavirus (SARS-CoV), and it belongs to the family Coronaviridae. Coronaviruses are enveloped, positive-sense, single-stranded RNA viruses and are characterized by spiky projections that resemble crowns. Four human coronaviruses are globally endemic and known to cause mild to moderate symptoms of cold. These strains account for 10-30% of upper respiratory tract infections globally [1,2]. On the other hand, novel coronavirus strains, possibly of zoonotic origin, pose the major threats to human life due to lack of pre-existing immunity that allows the virus to escape innate immune responses.
SARS-CoV-2, which belongs to betacoronavirus strain and is likely of zoonotic origin [3], has wreaked havoc by claiming more than 5.84 million lives and infecting over 414 million individuals (February, 2022). The basic reproduction number (R 0 ) of the virus has been estimated to be between 1.4 and 3.9 [4,5], and this is responsible for rapid spread of the disease. On 30th January 2020, SARS-CoV-2 was designated a Public Health Emergency of International Concern by the WHO, and on 11th March 2020, it was declared by WHO as a pandemic. Several vaccines are being used under Emergency Use Listing approval in a bid to control the Covid-19 pandemic [6]. However, vaccination is a time-consuming process and till date, 69% of the world population hasn't even received a single dose. Moreover, the long term efficacy of vaccine has not been established [7]. Worldwide Covid-19 incidences have not subsided, and several countries are seeing a spurt in Covid-19 cases [8,9]. Several countries are 1 3 showing a high rise of Covid-19 cases even after complete vaccination of significant part of their populations. Under these unprecedented conditions, the scientific community is focused on discovering a therapeutic cure for this disease. A major roadblock in finding a cure for SARS-CoV-2 is its high infectivity and the requirement of stringent handling protocols, which makes direct testing of candidate molecules a difficult proposition. Fortunately, a number of SARS-CoV-2 protein structures have been resolved and are available in public domain. This makes in silico studies a viable alternative to test the effectiveness of candidate molecules in binding and altering the functionality of critical SARS-CoV-2 proteins. Successful hits from in silico studies can further be tested and translated to in vitro and in vivo models.
Currently, remdesivir is the only drug approved by the FDA for treatment of Covid-19 [10]. WHO has actually discouraged the use of remdisivir [11]. Dexamethasone has shown beneficial effect in severe Covid-19 cases [12], but was actually associated with longer hospital stay in mild cases [13]. Hence, there is an urgent need for identification of small molecule inhibitors that can bind to and inactivate viral proteins, and thereby halt SARS-CoV-2 infection [14][15][16]. Existing antivirals like remdisivir, favipiravir, molnupiravir target viral replication enzymes. Expanding the target horizon to other SARS-CoV-2 proteins can provide novel mechanisms of preventing the progression of SARS-CoV-2. SARS-CoV-2 is highly efficient in evading the host immune responses [17,18] by causing suppression of host antiviral cytokines and also masking of viral RNA in order to avoid detection by innate immune receptors. Targeting the viral proteins responsible for such   immune evasion activities can be an exciting alternative strategy to render the virus susceptible to immune  responses. Large scale screening studies are usually carried out for identifying novel small molecule inhibitors. Two major issues with newly identified small molecule inhibitors are their (a) serum stability (b) potential toxicity. Often, high affinity binders identified from large scale screening studies cannot be eventually used owing to their poor serum stability and/or toxicity. A sensible alternative approach is to study small molecules with known anti-viral properties that are already recognized to be stable and non-toxic. In this direction, curcumin, active ingredient of the celebrated Ayurvedic herb turmeric, which is well documented to possess antiviral properties against a host of viral pathogens, could be a starting point for a viable pharmacological strategy to combat the virus [19][20][21]. Curcumin has been reported to work through a range of anti-viral mechanisms including inhibition of binding of enveloped viruses [19], interference of viral replication machinery or suppression of cellular signaling pathways essential for viral replication [22,23]. Computational studies have demonstrated the binding of curcumin to SARS-CoV-2 spike protein and its cognate human target ACE2 [24]. Unfortunately, despite having a favorable predicted binding affinity, curcumin suffers due to its poor bio-availability [25,26], and hence, multiple strategies are being developed to enhance the bio-availability of curcumin [26].
Mitocurcumin (MC), a derivative of curcumin, was synthesized by conjugating curcumin with two triphenyl phosphonium moieties (Fig. 1a, b) [27]. This had resulted in enhanced stability and uptake inside cells as compared to curcumin. Based on molecular docking studies, we had shown that MC could bind to TrxR (thioredoxin reductase) with higher affinity as compared to curcumin. Subsequent wet-lab experiments demonstrated 25-50 times more effective killing of A549 lung cancer cells by mitocurcumin as compared to curcumin [28]. MC has also been shown to have anti-bacterial properties (superior time-kill curve compared to ciprofloxacin against Bacillus subtilis) at physiologically achievable micromolar concentrations [29], and thus, is a much more realistic alternative to the parent molecule, curcumin. This differential activity may be due to the presence of triphenyl phosphonium ring that provides a distinct structural feature and might be helpful in conferring stability by favorable lodging at hydrophobic pockets. MC was shown to be stable in serum, with 93% recovery after 120 min of incubation [30]. Mitocurcumin has been evaluated found to be safe up to a dose of 250 mg/kg body weight (oral) in mice and up to a dose of 5 mg/kg body weight (intraperitoneal) in mice [30]. Oral dosing of MC lead to marginal enhancement in lymphocyte numbers in mice, which might be an added advantage in fighting SARS-CoV2. Evaluation of drug-likeness by Lipinski rule of five (Ro5), which estimates the drug-likeness of orally delivered drugs, showed that it satisfied two out of four conditions (less than 5 h-bond donors, less than 10 h-bond acceptors) [31]. It is worth mentioning that in 2016 and 2017, the average molecular mass of FDA approved drugs has been above the Ro5 cutoff of 500 Da [32]. The outcome of Ro5 of MC is in line with finding by Benet et al. who has shown that many natural compounds and derivatives are known to violate multiple Ro5 parameters and are likely to be non-orally active drugs [33]. This suggests that MC might be more effective in a parenteral route of administration. Indeed, parenteral mode of administration of MC (1 mg/kg intra-peritoneal injection) could protect mice from E. coli induced sepsis and death (Unpublished data). In light of its safe and stable nature, in silico studies were performed to identify its potential as an inhibitor of SARS-CoV2 proteins responsible for viral entry, replication, immune evasion and propagation.
This study provides insight into the putative application of MC as a small molecule inhibitor for SARS-CoV-2 proteins as well as its ability to serve as a potential drug to inhibit SARS-CoV-2 infection and replication in the host.

Ligand preparation
Curcumin and MC were drawn in ACD ChemSketch. DG-AMMOS was used to generate respective 3D conformers [34]. Ligands were imported in AutoDock Tools and prepared for docking by assigning the correct AutoDock 4 atom types, subsequently adding Gasteiger charges, merging nonpolar hydrogens, detecting aromatic carbons, and setting up the 'torsion tree'.

Molecular docking analysis and visualization
Molecular docking was performed using Autodock Vina [35]. PDB structures were imported in AutoDock Tools and prepared for rigid docking. H 2 O, ions and pre-existing ligands were removed from the structure to free up the active site wherever necessary. Charges were computed, non-polar hydrogens were merged and then the file were exported in pdbqt format. Grid box module was used to set the site of docking. For spike protein, ACE2 and RNA dependent RNA polymerase, multiple sites were chosen for docking. Site of docking were chosen from pre-existing literature or based on pre-existing ligands bound to the structure (Table 1). For example, in case of NSP3, target docking site was the binding site of natural ligand ARP (ADP ribose phosphate), while for NSP16, target docking site was the binding site of natural ligand SAH (S-Adenosyl homocysteine) (Fig. 2a, b). Grid sizes were kept at maximum of 27,000 cubic Angstrom. Docking was performed at default exhaustiveness which is calculated depending on the complexity of the ligand and structure.
After docking, the docked poses were visualized, and those having very high binding affinity but poor orientation with respect to active site were discarded. Highest affinity pose with proper binding site orientation were reported. Poses with high binding affinity (− 8 kcal/mol or lower values) were visualized with LigPlot+ to elucidate the type of interaction at the binding site [36]. Docked structures having predicted binding energy value lower than − 10 kcal/mol were taken up for molecular dynamics simulation.

Sequence analysis and alignment
Protein sequences were extracted from Uniprot/PBD and aligned using pBlast module of NCBI Blast. Aligned sequences were visualized with LAST Webservice.

Molecular dynamics simulation
Molecular Dynamics (MD) simulations were performed using GROMACS software package (version 2019.4) in double precision mode [37]. GROMOS 53a6 force field [38] was used to describe to proteins. Ligand topologies were obtained from the PRODRG server [39]. SPC model was used for defining water [40]. Simulations were run for protein alone (Protein-APO), protein with native ligand (Protein-LIG) and protein with mitocurcumin (Protein-MC). The docked complexes/APO-protein structures were initially minimized using steepest descent minimization for a permitted maximum value of 50,000 steps. A time step of 2 femtoseconds (fs) was used, and NVT equilibration was carried out for 100 picoseconds (ps) with a modified Berendsen thermostat. Further, NPT equilibration was performed for 100 ps using Berendsen thermostat for pressure coupling. Particle Mesh Ewald (PME) method [41] and force-switching scheme were used for long range electrostatics calculation and van der Waals interactions, respectively. Short-range van der Waals cutoff was set at 1.2 nm. Position restraint was applied for protein and ligand in the equilibration steps. For the production MD run, Berendsen thermostat with velocity rescaling [42] and Parrinello-Rahman barostat [43] were used. Production MD was run for 100 ns, with energy and trajectories being saved every 10 ps.
Inbuilt modules of GROMACS were used to calculate root-mean-square deviation (RMSD), radius of gyration (R g ), solvent-accessible surface area (SASA), hydrogen bond (HBOND) and interaction energies (Columbic short range and Leonard-Jones short range).

MM/PBSA analysis
Molecular Mechanics-Poisson-Boltzmann Surface Area (MM/PBSA) method was used to estimate the binding  [44]. This method estimates the total binding energy of a complex as the difference between total free energy of the complex and the sum of the total free energy of the individual components. Total free energy is estimated from multiple components: molecular mechanics-based calculations using Coulomb and Lennard-Jones potential functions are used to estimate the electrostatic and van der Waals forces; polar solvation energy is estimated by solving the Poisson-Boltzmann (PB) equation; non-polar solvation energy is calculated by using SASA, Solvent accessible area (SAV) or Weeks-Chandler-Andersen (WCA) theory-based models [45].
g_mmpbsa tool was used for calculation of MM/PBSA [44,46]. Nonpolar solvation energy was estimated by SASAonly model which assumes that nonpolar solvation energy is linearly proportional to SASA.

Mitocurcumin binds to active site of NSP3 with high affinity and may prevent NSP3 mediated suppression of anti-viral cytokines
NSP3 is the largest protein encoded by the coronavirus genome with its functions ranging from replication and transcription to nucleocapsid binding and immune evasion [47]. The macrodomain1, also known as NSP3b is known to play essential role in countering host innate immune responses.
Molecular docking results show that MC can bind to "ADP ribose phosphate (ARP)" binding site of NSP3 with a high binding affinity of − 10.3 kcal/mol (Table 2). This corresponds to a predicted K i of 28 nM. MC sits in the ARP binding pocket (Fig. 3a, b) and interacts primarily through hydrophobic interactions (Fig. 3c). The triphenyl groups at the two ends are stabilized by interaction with Gly47, Val 49, Gly 58, Pro 125, Gly 130, Val 155 at one end and by interaction with Asn 99, Val 77, Val 41, Asn 40, Tyr 42, Pro 74 at the other end. The core part is stabilized by Lys44, Phe 132, Ile 131, Gly 46, Ala 50, Ala 38. The presence of triphenyl groups seems to confer the added stability as compared to curcumin (− 8.4 kcal/mol), because curcumin lacked a large proportion of the hydrophobic interactions due to lack of triphenyl groups (Figs. 3c and 4a). On the other hand, Curcumin-NSP3 interaction showed a number of h-honds (Gly130, Ala 129, Asn40, Glu104) to stabilize the interaction (Fig. 4a).
MD simulation was carried out to estimate the stability of the complex between MC and NSP3 (Fig. 5). RMSD was calculated for the backbone of the protein relative to initial structures (Fig. 5a). NSP3-APO, NSP3-MC, NSP3-ARP had average RMSD of 0.160, 0.176 and 0.158 nm, respectively, suggesting that the complexes were stable. Both NSP3-MC and NSP3-ARP showed an initial increase in RMSD up to 60 ns, after which they stabilized to ~ 0.181 nm and ~ 0.180 nm, respectively. Radius of gyration provides an estimate of the compactness of the protein complex. Radius of gyration remained stable throughout the MD run, and average values were estimated to be 1.52, 1.52 and 1.51 nm for NSP3-APO, NSP3-MC and NSP3-ARP, respectively  . 5b). This shows that the compactness of the protein complex is not hampered by binding of either the natural ligand ARP or test molecule MC. Number of internal hydrogen bonds in the protein was measured. A similar pattern was noted for all the 3 forms, providing an average 126, 126 and 128 internal hydrogen bonds in NSP3-APO, NSP3-MC and NSP3-ARP, respectively (Fig. 5c). Average solventaccessible surface area (SASA) was 87.7, 87.1 and 86.4 nm 2 for NSP3-APO, NSP3-MC and NSP3-ARP, respectively (Fig. 5d). No significant alteration in number of hydrogen bonds and SASA was observed suggesting that the complex between MC and NSP3 was a stable one, and comparable to the one between NSP3 and its natural ligand, ARP. The average number of h-bonds between NSP3 and natural ligand ARP was 6.36 while the average number of h-bonds between NSP3 and MC was 1.48. This suggests that the interaction between MC and NSP3 was more dependent on hydrophobic, coulombic and other electrostatic interactions, unlike the interaction of ARP with NSP3 which was more dependent on h-bonds. Short range Lennard-Jones and Coulombic interaction energy between NSP3 and MC/ARP has been reported in Table 3. Energy estimates suggest that MC and ARP can bind equally well at NSP3 active site, with MC having a slight edge. Further, MMPBSA calculations were performed on both the complexes. Net binding energy of interaction between NSP3 and MC was − 42.93 ± 1.8 kcal/mol, which indicates that MC and NSP3 can form a strong complex. The contributions of different energy terms toward the binding energy have been detailed in Table 4. Net binding energy of interaction between NSP3 and ARP stood at − 73.86 ± 0.95 kcal/ mol.
Studies have shown that ADRP activity deficient NSP3 macrodomain displays an increased sensitivity to the antiviral effect of alpha interferon [48]. NSP3's ADP-ribose-1′phosphatase activity is essential for evading host immune responses. Mutant forms deficient in ADRP activity elicited early and enhanced interferon (IFN), interferon-stimulated gene (ISG) and inflammatory cytokine release. It has also been proposed that de-poly-ADP-ribosylation or de-mono-ADP-ribosylation activities of NSP3 might also be playing an important role in immune evasion. The loss-of-function mutations in the ADRP domain of SARS-CoV NSP3 protein have resulted in reduced viral load in mice [49]. Hence, blocking the active site is expected to have the same effect as a loss-offunction mutation. We had also carried out pairwise sequence alignment of the macrodomain of NSP3 from SARS-CoV-2 and SARS-CoV (Fig. 6). Sequences showed 124/168 (74%) identity and 144/168 (85%) positive matches. The catalytically essential residues were conserved between SARS-CoV-2 and SARS-COV. Hence, it is expected that NSP3 of SARS-CoV2 will respond in a similar fashion. Binding of MC to NSP3 can potentially block its enzymatic activity and thereby preventing SARS-CoV-2 from suppressing the anti-viral interferon response. This can lead to enhanced viral load clearance by innate immune system and reduce the extent of viral pathogenesis in individuals without antigen memory i.e. first-time affected individuals.

Mitocurcumin was also found to bind to the active site of NSP10-NSP16 complex and may prevent NSP16 mediated masking of viral RNA strand
Viral RNA recognition through RIG-I (retinoic acid-inducible gene I) and IFIT (Interferon Induced proteins with Tetratricopeptide repeats) family of proteins is the major mechanism through which the anti-viral pathways are induced in response to foreign RNA in cytoplasm [50,51]. RIG-I is a major stimulant of anti-viral interferon production. The basis of detection of viral RNA as 'foreign' lies in the characteristic 5′ cap structural modifications. Specifically, N-7 and 2′-O positions of the 5′ guanosine cap are crucial as these remain methylated in eukaryotic RNA. IFIT1 is an interferon stimulated gene which senses lack of 2′-O-methylation and inhibits viral RNA [52]. Viruses can escape this detection by modifying their 5′ cap ends using 2′-O-methyltransferase activity, and this allows them to evade the immune system [53]. NSP10-NSP16 complex of SARS-CoV2 is responsible for 2′-O-methyltransferase activity and subsequent evasion from IFIT1 mediated    (Table 2). It is stabilized by a mix of hydrogen bonds and hydrophobic interactions (Fig. 7c). Hydrogen bonds through Cys6913, Tyr6930   (Figs. 7c and 4b), and this could be responsible for comparatively lower binding affinity (− 7.5 kcal/ mol) of curcumin.
Further, RMSD analysis of docked structures revealed that NSP16-APO and NSP16-MC complex had comparable RMSD values of 0.23 and 0.21 nm, respectively (Fig. 8a). Beyond 40 ns time point, RMSD became stable for APO and NSP16-MC. Interestingly, the SAM bound NSP16 backbone showed relatively higher RMSD fluctuations with peak value of ~ 0.38 nm at 70.78 ns, following which it stabilized at ~ 0.3 nm. The values were within an acceptable range and shows that the ligand bound form of NSP16 was stable. After a slight decline in initial 10 ns of the MD run, radius of gyration displayed minimal fluctuation and produced comparable average values of 1.89, 1.89 and 1.9 nm for NSP16-APO, NSP16-MC and NSP-SAM, respectively (Fig. 8b). Low radius of gyration values suggests that the ligand bound form of NSP16 continued to stay compact. Average number of internal hydrogen bonds was 223, 219 and 218 for NSP16-APO, NSP16-MC and NSP-SAM, respectively (Fig. 8c). Hydrogen bond count showed a slight increase in the initial 10 ns of the run, following which it became stable. SASA analysis revealed slight dip in the solvent-accessible surface area in the first 10 ns of the run. Then it became stable and produced average SASA values of 152.2, 151.3 and 150.3 nm 2 for NSP16-APO, NSP16-MC and NSP-SAM, respectively (Fig. 8d). The concomitant decrease in radius of gyration, decrease in SASA and increase in internal hydrogen bonds in the first 10 ns of the run suggests that the new hydrogen bonds caused a compaction of the protein structure and subsequently decreased its exposure to the solvent. The comparable results of APO and ligand bound forms suggest that binding of MC to NSP16 is a stable one. In between NSP16 and natural ligand SAM, the average number of h-bonds was 4.8, while between NSP16 and MC, the average number of h-bond was 0.27. The interaction between MC and NSP16 seems to be heavily reliant on hydrophobic, coulombic and other electrostatic interactions, in contrast to the h-bond dependence between SAM and NSP16. Comparison of short-range interaction energies for NSP16-MC and NSP16-SAM complex shows that they had similar values, with SAM holding a slight edge over MMPBSA analysis estimated net binding energy for interaction between NSP16 and MC to be − 118.51 ± 1.74 kcal/ mol. This indicates a very strong binding between MC and NSP16. The energy terms have been detailed in Table 4. Interaction between NSP16 and its natural ligand SAM was slightly weaker with an estimated net binding energy of − 88.86 ± 1.451 kcal/mol.
Coronaviruses possess 2′-O-methyltransferase activity in the NSP16 protein [54] which allows the viral RNA to evade innate immune system detection by modifying their 5′ cap ends. NSP16 works in concert with NSP10 and uses S-Adenosylmethionine as the methyl group donor to modify their 5′ cap by 2′-O-methylation. The binding of MC to S-Adenosylmethionine binding site would deprive NSP16 of the essential methyl-group donor and thus prevent the 2′-O-methyltransferase activity. In absence of 2′-O-methyl group attachment to the 5′cap of viral RNA, it would get detected by host immune system, thereby leading to enhanced clearance of viral load.

Mitocurcumin binds to SARS-CoV2 targets with a higher affinity as compared to curcumin
Molecular docking data showed that MC can also bind to main protease of SARS-CoV-2 with a predicted binding affinity of − 8.4 kcal/mol, and hence may also act as an inhibitor of the same (Fig. 9a, b). The binding is stabilized by hydrogen bond with Ser 46 and hydrophobic interactions with Thr 25, Met 49, Leu 26, His 41, Cys 145, Asp 187, His 163, Gln 189, Leu 50, Glu 166, Leu 167, Pro 168, Ala 191, Thr 190, Glu 47. The predicted binding energy was higher as compared to curcumin (− 6.6 kcal/mol). Importantly, interaction with catalytically active residues Cys 145 and His 41 can potentially inhibit the activity of Main Protease (Fig. 9c). MC was also found to interact with active site of other SARS-CoV-2 proteins with varying degrees of affinity (Table 2). MC could bind to spike protein salt bridge site (− 5.9 kcal/mol), ACE2 salt bridge site (− 5.8 kcal/mol), nucleocapsid (− 7.2 kcal/mol), NSP15 (− 6.8 kcal/mol) and Rdrp RNA binding site (− 6.7 kcal/mol). It must be noted that in most cases, the affinity of MC was greater than the parent molecule curcumin, probably due the tendency of distinguishing triphenyl groups to get lodged in hydrophobic pockets and confer extra stability compared to curcumin (Figs. 3c and 4a; Figs. 7c and 4a). Based on these observations, MC appears to be a better choice than curcumin as an anti-viral molecule against SARS-CoV-2 infection.

Conclusions
Our molecular docking coupled with MD simulation results suggest that MC can stably bind at the active site of NSP3 and NSP16. By virtue of blocking the active site, this binding has the potential to down-regulate the ADP Ribose phosphatase activity and 2′-O-methyltransferase activity of NSP3 and NSP16, respectively. It is an exciting finding as this demonstrates the potential of a single molecule to target two crucial pathways simultaneously. These interactions with NSP3 and NSP16 hold the possibility to counteract two crucial immune evasion strategies of SARS-CoV-2, rendering it vulnerable to host innate immune response mediated clearance. To the best of our knowledge, this is the first MD simulation report on targeting the two arms of SARS-CoV2 immune evasion using a single molecule. The predicted K i values for NSP3 and NSP16 stand at 28 and 23 nM, respectively. The stable, nontoxic nature of mitocurcumin coupled with the physiologically achievable predicted K i values makes it a realistic candidate molecule against SARS-CoV-2 infection.