Full Text Article

Introduction

Viruses are parasitic microorganisms whose expressed machinery is both economical and versatile, considering its multiple essential and accessory functions (many not yet understood) with a minimum of genetic material and proteins. One example of this switchblade molecular machinery are the non-structural proteins (NSPs), which only subsists as long as the virus stays inside the cell, where they actively participate in viral replication, and are not present within the inert particles [1]. Among the NSPs expressed by chikungunya virus (CHIKV), nsp3 is known to be the most difficult to understand [2]. But even ten years after the structural elucidation of its first functional domain [3], few in silico studies have been carried out to date.

The first CHIKV three-dimensional protein structures obtained from crystallographic data were nsp3 Macrodomains, also known as Macrodomain-X, a functional fold that is evolutionarily very conserved and found in all kingdoms of living beings [4]. In CHIKV genome, it is comprised from the N-terminal portion of nsp3 to its residue 160, has a topology of the α/β/α sandwich composed of 4 α-helix and 6 β-sheets, the all β-sheets are located in the middle of the protein, surrounded by three α-helix on one side, and one α-helix on the other. While the understanding of the specific role of nsp3 viral cycle is still fuzzy, the function of its macrodomain has been gradually revealed [5-7].

Evidence about nsp3 macrodomain function came from information that homologues were able to recognize ADP-ribose (ADPr) nucleotides and remove their first phosphate group, the 1″-phosphate phosphatase activity. Based on this, Malet, et al. (2009) crystallized the macrodomain (both CHIKV and VEEV) with good resolution (≤ 2Å) under three different conditions: apo form - in complex with ADP ribose - an adenosine diphosphate nucleotide linked to an additional ribose in the distal tip, at its second phosphate; and also in holo form - in complex with an RNA trimer (AAA) - respectively corresponding to the PDB ID 3GPG, 3GPO, 3GPQ codes. There was no difference between the apo and holo conformations of the protein. Furthermore, the analysis of crystals with different ligands, in principle, allows identification of two sites: ADPr site - above the 2, 4 and 5 β-sheets and surrounded by the loops that connect β2 to α1 and β5 to α3; and RNA trimer site, where nucleotide fit in the same site as the ADPr but, additionally, depended on the electrostatic attraction of the phosphate groups with a positively charged region composed of the G112, V113, and Y114 residues.

Site-directed mutagenesis outlined essential residues for phosphatase activity [8], and the analysis of the complexes allowed identification of important residues in substrate binding, which could be involved in catalysis. Thermal shift assays showed the macrodomain specificity by adenine when compared to guanine; by a di-phosphate instead of one or three; and by distal ribose, instead of glucose, phosphate, nicotinamide or nothing. Following studies have shown, however, that phosphatase activity is not very pronounced in alphaviruses [9], and that there is another, much more notable function of this target: besides binding to ADPr units, they are able to recognize these sugars in glycosylated proteins and remove them, an activity known as ADP-ribosyl hydrolase or even mono-ADP-ribosyl hydrolase “MARilase” [6, 10, 11].

Protein glycosylation with ADP ribose (ADPr), ribosylation, is a post-translational modification carried out by ADP-Ribosyltransferases (ARTs), which transfer the Nicotinamide ADPr monomer Adenine Dinucleotide (NAD+) to the side chain of amino acids, an action which may result in allosteric modulation - comparable to the phosphorylation or methylation of nucleic acids [12]. Dynamic and still poorly understood, the fact is that the reversal of ribosylation of proteins (still unknown) is an essential activity for CHIKV replication in both mammalian and mosquito cells. And it has already been demonstrated with the isolated nsp3 md the hydrolysis of ADPr in Aspartate and Glutamate residues, what does not happen in Lysine residues, and was also concluded that MARilase activity (and not just ADPr affinity) is critical for viral replication in cell lines, even though mutants with greater affinity for ADPr have been more virulent [6].

All of this reflects the importance of studying macrodomain structure for the development of anti-CHIKV drugs. Considering the studies with CHIKV proteins, certainly the potential of nsp3 macrodomain in drug development has been neglected, but recently more attention has been given due to the great structural similarity of this protein in different viruses of families Togaviridae, Hepeviridae and Coronaviridae , of such relevance to public health as the pandemic SARS-Cov2. Therefore, it has been suggested that macrodomain is an strategic molecular target for the development of broad-spectrum antivirals, and the large amount of data generated due to SARS-Cov2 boosted studies involving the macrodomain also of other viruses, such as CHIKV [13].

It is noteworthy the fact that there is a lack of inhibitors directly associated with CHIKV targets, especially nsp3 md. Beyond the substrate (ADPr), the unique in vitro validated compounds with inhibitory activity against CHIKV nsp3 md are pyrimidone derivatives obtained very recently by a fragment-based drug design which proposed the 2-pyrimidone-4-carboxylic acid as a scaffold for nsp3 modulation, from which a derivative showed a viral replication inhibition of IC50 of 23 μM in human fibroblasts. Thus, there is still a hindrance to the implementation of ligand based virtual screening evaluations. On the other hand, there are a series of virtual screens already carried out based on the structure of CHIKV macrodomain, but that were carried out in different software and approaches that make it difficult to compare them. Because of this, in this work we aim to investigate the previously reported virtual hit putative CHIKV nsp3 md inhibitors and the in vitro CHIKV inhibitors (mainly those that are not associated with targets) in relation to their affinity to CHIKV and SARS-Cov2 macrodomains.

Here we 1) review the structural aspects of CHIKV nsp3 macrodomain “nsp3 md” as well as the contributions of in silico and in vitro assays related to this target/virus, 2) evaluate the affinity of most promising compounds from all over these studies with CHIKV and SARS-Cov2 nsp3 mds by the same method, and 3) use the protein ligand interaction profile of the virtual hit compounds to understand the specificities of each modulation pattern and 4) build a virtual screening strategy for the search of new putative antiviral compounds.

Material and Methods

Subset selection and preparation

Initially, an extensive search was accomplished for studies using the in silico molecular docking and molecular dynamics methods with CHIKV nsp3 macrodomain, both in indexed repositories and grey literature. Then a new search was performed in the same literature sources for in vitro inhibitory compounds for CHIKV. The information on the in silico methodologies used was collected and presented in the form of tables (with corrected and standardized numbering or residues), and the 16 top virtual hit molecules from those studies were drawn and had their smiles generated on the Marvin Sketch software [14], for “virtual hits” subset preparation. After the selection of CHIKV in vitro inhibitors (IC/EC 50 < 0.1 mM = 139 compounds) (see supplementary Table 1), were also selected the scaffold pyrimidine fragment proposed as the minimal essential structure for CHIKV nsp3 inhibitory action; and the inhibitor of nsp3 macrodomain of SARS-Cov2.

All these structures were drawn and had their smiles generated on the Marvin Sketch software. Finally, a subset of 324 molecules bounded in hydrolases structures from Protein Data Bank (https://www.rcsb.org/) which were similar to ID 3GPO – CHIKV nsp3) and a subset of 652 molecules retrieved in the same database bound to nsp3 viral structures were selected (crystallographic residues and ions were discarded).

The SMILE strings from all molecules were converted into SYBYL-X 2.0 mol2 3D format [15] and were prepared for docking procedures by addition of the hydrogen atoms and calculation of the Gasteiger-Hückel charges in Open Babel v. 3.2.1. [16]. The macrodomain structures from CHIKV and SARS-Cov2 used in calculations were obtained, respectively in chain A from PDB IDs 3GPO 1.9Å and 7KQP 0.88Å in Protein Data Bank. The receptor preparation with Hydrogen atoms and Gasteiger partial atomic charges (ff14SB) was set in the DockPrep module at Chimera 1.10.1 program [17].

Docking procedures

Docking was carried out in DOCK 6.9 software package, using the molecular surface file generated by Display Midas System (DMS) [18] in Chimera 1.10.1 program to access receiver solvent and the negative image of the receptor. The orthosteric site was delimited by crystallographic ligand (ADPr) position, as the closest spheres (radius: 10Å) were selected, using SPHGEN and SPHERE_SELECTOR programs Then, a calculation box 8Å distant from these selected spheres was set in both CHIKV and SARS-Cov2 macrodomain structures. Finally, the physical/chemical properties were calculated by the GRID program [19] in its standard configuration, and the GridScore function was used to rank the structures.

The docking parameters were evaluated by performing the root-mean-square deviation (RMSD) between the crystallographic ligand and the best classified pose by the program.

Considering a satisfactory value an RMSD of less than 2 Å [20]. Accordingly the GridScore value, the top scored molecules were selected for molecular dynamics simulation.

Molecular Dynamics

CHIKV macrodomain complex with top1 anti-chikv compound was submitted to molecular dynamics. First, the protonation of ionizable protein residues was calculated using PROPKA server [21], to reproduce the cell pH of 7.2 and the protonation states were manually set in Gromacs environment [22] using PDB2GMX module, that in addition to generating the topology also converted receptor coordinates into Gromacs gro file using PDB2GMX module. Ligand topology was generated in PRODRG server, as well as the conversion of its coordinates demanded by GROMOS united atom force field.

The system was placed in the center of a dodecahedron simulation box of size: 4.105/ 3.886/ 4.010 (nm) and solvated with 9054 SPC-E water models, and 1 Cl counterion was added to neutralize it. Subsequent minimization steps were carried out by Steepest Descent and Conjugate Gradient algorithms until the convergence in the energy. Then using the Verlet integrator a nanosecond of molecular dynamics was carried out for carefully heating the system, with restriction of protein atoms, until the 300K were reached. Finally, the restraints were removed in a run of 100ns of molecular dynamic simulation. For each atom, long-range the closest interactions in 14Å were considered, and electrostatic interactions were treated using the PME - Particle Mesh Ewald method.

The stability of system was evaluated through analyses of RMSD (Root Mean Square Deviation) and RMSF (Root Mean Square Fluctuation) of heavy chain atoms, RG (Radius of Gyration) calculation and the calculation of h-bonds permanence between ligand and receptor over the simulation. Also, the module of mm-gbsa was employed to calculate the free energy binding of the complex.

Intermolecular interactions analysis

The central pose from the most populated cluster of the molecular dynamics production phase was selected for Protein-Ligand Interaction Profile calculations in PLIP server [23]. And the visual inspection of complexes was made using Pymol [24].

Results and Discussion

Docking studies with CHIKV nsp3 Macrodomain (selected from literature review)

The structural features of protein-ligand interactions of virtual hits from nsp3 macrodomain docking studies, selected from a literature review, are presented in Table 1 (Brinda, et al. 2019; Kumar, et al. 2020; Nguyen, et al. 2014; Oo, Hassandarvish, Chin, Lee, Bakar, et al. 2016; Puranik, Ninad V.; Rani, Ruchi; Singh, Vedita Anand; Tomar, Shailly; Puntambekar, Hemalata M.; Srivastava, 2019; Sangeetha, et al. 2017; Santhanam Vijayasri, 2017; Seyedi, et al. 2016; Singh, et al. 2019; Vora, et al. 2019) (Table 1) [25-34].

Although the selection of starting structure for virtual tests is crucial for the reliability of results, there are incompatibilities within these choices. It is known that molecular targets when crystallized in a complex provide important information about the interacting residues, in addition, a holo structure reflects the transition state of an enzyme, in which both central and lateral chains assume a conformation which allows to “receive” the ligand, the induced fit. There was generally a consensus on the choice of the holo structure for in silico calculations [35], but almost half of the studies with the macrodomain used the apo structure. This choice was probably based only on the criterion of crystallographic resolution, lower in 3GPG (Apo 1.65Å) than in 3GPO (in complex with ADPr - 1.9Å), a minimal difference, which would not have justified the choice of 3GPG, since the two resolutions are sufficiently detailed for virtual screenings. Even so, the structural overall similarity between both structures (0.19Å rmsd) may have influenced the choice for 3GPG. The structural alignment of apo and holo states shows that among the ADPr interacting amino acids residues cited by the author, there is variation only in G30 and D31.

The structure in complex with RNA (3GPQ) was not chosen, because even though it is able to bind to polynucleotides, the macrodomain does not show PARylase activity (Poli ADP ribosylase). The same justification should apply for 4TUO, although a single study used this structure [25].

In order to consider other putative modulable regions, some studies have used the blind docking procedures, in this case the simulation box is large enough to cover the entire protein to consider the possible existence of other interaction sites. However, in a study which used the Diversity Set II Bank for screening, a chemical and structural diverse set of compounds selected by a therapeutic development program of the North American National Cancer Institute, blind strategy did not prove to be advantageous. The best scored virtual hits were located at the ADPr site and the other proposed site presented only hydrophobic interactions with ligands attached to it [27]. Also in studies which selected two regions for docking calculations, all ligands bound with more affinity within the large catalytic cavity [31] detected by CavityPlus.

For directed docking strategies, most studies with nsp3 md at least suggested the use of the crystallographic ligand (ADPr) as defining the simulation region, or relied on the help of software to predict the binding site [33], however, in most of them no details were given about the docking region. The amount of ligands evaluated diverged according to the objective of each study. Virtual screenings were carried out mostly with natural compounds, however relatively few compounds were evaluated in relation to the available databases [36]. Other in silico studies proposed which would be the most likely target to known inhibitors of chikv replication, yet the majority of CHIKV inhibitors remain with no associated biological target of action [37,38].

The diversity of computer programs and methods employed for molecular docking studies makes it difficult to compare results from different authors, and often the same research groups optimize the compounds they have prioritized. Each study showed in table 1 mapped the interactions of virtual hits complexes, important residues are marked. This analysis allows to establish which ones are important in a desired modulation, and can help to define the proper region for next calculations. The authors highlighted the important role of polar and aromatic interactions, and especially the hydrogen bonds, in which the macrodomain residues participated predominantly as donors, which for example, in ADPr (substrate) complex stabilizes the phosphate groups and distal ribose.

Until recently there were no known inhibitors for the chikv nsp3 macrodomain, thus the validation of scoring functions was not possible through the discrimination between decoys (false positives) and true ligands. Therefore, the only option available to assess the predictive capacity was to re-dock the crystallographic ligand to the active site after submitting it to minimization. The smaller the variation is between crystallographic ligand and the re-docked one, the greater the accuracy of the technique and parameters used is supposed to be. However, the re-docking values were only indicated by Nguyen (2014). In many cases, in silico calculations were employed just to complement in vitro assays, but it is worth noting that to publicize detailed information on the methodology used is fundamental for reproducibility, comparison and reliability of the results.

Molecular Dynamics Simulations with CHIKV nsp3 Macrodomain (selected from literature review)

Some limitations of docking technique (for example, the fact that the receptor structure is generally kept fixed and the non-consideration of a solvated system) can be overcome through molecular dynamics (MD) simulations of prioritized virtual hits. MD simulations are extremely useful to assess the permanence of interactions identified in docking, as well to assess the adaptations of the complex induced fit. Among the ten studies that used accomplished docking, only six implemented studies with MD simulations (Kumar, et al. 2020; Nguyen, et al. 2014; Oo, Hassandarvish, Chin, Lee, Abu Bakar, et al. 2016; Santhanam Vijayasri, 2017; Singh, et al. 2019; Vora, et al. 2019), and one more which evaluate 3GPO crystallographic complex (with ADPr) (Rungrotmongkol, et al. 2010) [39].

As in the docking studies, the methodologies used for MD were as varied as possible (Table 2). The description of system preparation was not clear in the majority of these studies. Authors did not mention how protonation of the ionizable residues was calculated and simulated, only Santhanam (2017) provided information on how ligand topology was generated (by PRODRG server). MD studies employed the most popular water models, and in general, used cubic simulation boxes. Several MD rmsd trajectories with macrodomain showed that most ligands were able to stabilize the protein. An intriguing exception occurred in Vora (2018) results, whose rmsd values were discrepant in relation to all the others. This could indicate that Vora et al probably showed a wrong description of the system, or Mulberroside ligand caused an accentuated conformational change. In any case, the disorganization of the graphs and the tiny description of results prevent a reliable interpretation (Table 2).

The authors were mainly concerned in describing the polar interacting residues. Rungrotmongkol (2010), for example, highlighted the importance of other residues for ADPr binding such as I11, N31, V33, C34 and S110, whose interaction was not noticed by Malet (2009) in the crystal complex. The study with the longest simulation periods provided interesting information, for example, residues L108, A22, V69 and R144 (two of them with permanent interactions throughout the simulation), respectively, had the lowest rmsf values, which corroborates their importance in stabilizing the ligands within the site, and correspond to the those pointed out by several authors as important for binding, in showed an extremely similar pattern both in analysis of 500ns and 2μs (Kumar, et al. 2020). This makes it possible to say that simulations shorter than 500ns may be sufficient to evaluate the macrodomain complexes interaction.

Another highlight derived from MD studies was that D10, a conserved residue in all macrodomains and important for adenine binding (Malet, et al. 2009) appeared few times in the list of other ligands interacting residues, both in the docking procedures and in after MD simulations, as well as the identification of residue D145, recognized only in two docking evaluations (Puranik, Ninad V.; Rani, Ruchi; Singh, Vedita Anand; Tomar, Shailly; Puntambekar, Hemalata M.; Srivastava, 2019; Seyedi, et al. 2016). The other residues pointed out by MD results were recurrent.

In all simulations performed in the Amber software, the binding free energy for the complex was calculated, using the method that considers molecular mechanics of Born’s generalized surface area (MM-GBSA). The obtained values were -40.05 Kcal/mol for CMPD-104 (Kumar, et al. 2020); -28.18 Kcal/mol for Hesperetin (Oo, Hassandarvish, Chin, Lee, Bakar, et al. 2016); -46.18, -34.47, -34.22 and -39.85 kcal/mol for the binders derived from Noscapine (Singh, et al. 2019); and -6.2 Kcal/mol for ADPr (crystallographic ligand – substrate) (Rungrotmongkol, et al. 2010) [39], which shows that these molecules have a high potential to inhibit the activity of macrodomain in vitro. Finally, taking into account only the best in silico scored molecules, the molecules which comprise the subset of “virtual hits”. Some of these molecules were shown in Table 3.

Selected In vitro studies with CHIKV infected cells (selected from literature review)

More than fifty inhibitory CHIKV studies in mammalian cells have already been carried out to date. The majority did not specify the protein target of compound modulation, at most, they indicated what is the replication stage of the inhibitory activity (e.g.: entry, post-entry events, virus assembly, vesicle budding, early and late virus cycle events, transcription and translation, PKC metabolic via, PKC independent via, AKT-phosphorylation, cholesterol trafficking). Information such as EC50, IC50, CC50 and SI not always was indicated, for example, as in the evaluation of plant extracts. Few of them used in silico analysis to indicate, within the available viral structures and its proposed active sites, which has more affinity to the active compounds.

Docking calculations

The redocked crystallographic ligand poses (figure 1) showed the following RMSD: CHIKV = 0.34Å; SARS-Cov2 = 0.29Å, which indicate very satisfactory results, as values are much less than 2Å (Verli, 2014), thus the DOCK program was able to reproduce the conformation of ADPr in the active site of crystal complexes with minor variation, thus it can generate reliable poses in virtual screening. The quality of the structural models has a great influence on this. The 2D protein-ligand map from both complexes (CHIKV and SARS-Cov-2 nsp3 macrodomains) shows the subtle differences between their active sites and the substrate (Figure 2), and although the active site is not entirely conserved, there is high similarity and structural positioning correspondence of the residues, the only difference is an additional interaction with the central ribose by CHIKV nsp3.

Marked residues: Adenine site = yellow; distal ribose site = red; central ribose and phosphates site = green

The structural similarity of the active sites is also evidenced through the remarkable resemblance of the results from all subsets evaluated. This reinforces the importance of drug design directed for Macrodomains in a perspective of a broad antiviral range. Some of these molecules need to be evaluated against nsp3, for confirmation of its putative target modulation. Here, the 5 top scored molecules from each subset are presented below in table 3, and the interacting residues from each one with both Macrodomains (CHIKV and SARS-Cov2) are presented in supplementary Table 2.

Considering the selected molecules from in silico studies involving CHIKV macrodomain, none of them showed affinity for the Macrodomains active site than the substrate, though Nguyen14-T1 presented H-bond interactions with some of the same ADPr interacting residues in both proteins.

Within the best scored selected anti-CHIKV compounds, three are natural diterpene esters compounds which were able to inhibit CHIKV and HIV replication in vitro at the nanomolar level and have been proposed to act through modulation of protein kinase C isoenzymes (PKC) (Nothias-Scaglia, et al. 2015). In addition to antivirals which have been associated with other mechanisms of action, such as Akt-phosphorylation and at early stages of infection, as Suramin (Salgado-Benvindo, et al. 2020). The top1 against CHIKV macrodomain, with reference “NothiasScaglia15-18”, a phorbol triester (12,13-O,O′-dinonanoylphorbol-20-homovanillate) of 0.6 uM/L EC50 anti-chikv activity, could bind very tightly to the large macrodomain cavity, comprising most of the regions accessed by other ligands and substrate, with a energy scoring value (-83.04 kcal/mol), stronger than the obtained in crystal complex with the substrate (-79.58 kcal/mol).

The analyzes of CHIKV and SARS-Cov2 nsp3 unique inhibitors served as a basis for interacting residues evaluation of other ligands, beyond the substrate, and it is important to highlight that even with differences in binding energies, they were also able to reduce viral replication in vitro.

The results showed a diversity of scaffolds structurally different from the substrate, many of them, considering the in silico analysis, even showed greater affinity to the active site of both macrodomains (from CHIKV and SARS-Cov2) or very close energy values to the value found for substrate. Some of these compounds (detailed in figure 3) have been already confirmed as putative antivirals by in vitro assays, which reinforces the importance of further evaluating this potential modulative capacity against viral nsp3 macrodomains.

Considering the 100ns of MD simulation of CHIKV macrodomain in complex with top1 anti-chikv (NothiasScaglia15-18) compound, the production phase after 30ns (average of SD = 2.3) was selected based on rmsd trajectory analyses, because the rmsd values showed variations less than 2Å (Figure 4), besides that all over the simulation the ligand stayed in the macrodomain cavity. Thus, the free binding energy complex calculated using the production phase was -133.49 kJ/mol, which corroborates that this ligand showed a good affinity with the macrodomain. Finally, through a cluster analysis the frame at 36050ns was selected for mapping the ligand-receptor interactions.

When the RMS (root mean square deviation) is evaluated by residue, loop residues present variations above 3Å. However, the most important residues for ligand stabilization, which interacted through their heavy chains (R144 and D145) are 1Å around (Figure 5). This also is evidenced, when analysing their interaction frequency throughout the simulation (Figure 6).

Also, ILE11, ASN14, ASP31, ALA36, VAL37, TYR142, ASP145

residues were important for hydrophobic interactions, whose majority are present in other ligand-receptor interacting residues obtained by previous in silico evaluations.

After the simulation, the compound was able to form new hydrophobic and polar interactions, which indicates a better accommodation than previously indicated on docking result. The hydrogen bond formed between ligand and ARG144 residue was very permanent over the simulation (Table 4), and it hadn’t even been indicated by docking. In fact, some of the relevant residues for this interaction (figure 6 - right) are different from those found in the substrate active site. But, previous MD simulation studies found in literature show that ARG144 residue was indeed important in the stabilization of the substrate at the active site (Rungrotmongkol, et al. 2010), and had been pointed out by MD studies with other compounds (including an evaluation throughout 2000ns of simulation) (Kumar, et al. 2020; Nguyen, et al. 2014). Considering the residues which contributed to stabilizing the ADPr phosphates in the substrate active site, VAL113 and ASP31 also were detected in our simulation, but not so frequently. So the putative displacement of the substrate in a competitive inhibitory action of NothiasScaglia15-18 compound should be in vitro validated.

Conclusions

The lack of treatment and vaccines for chikungunya virus infection is a challenge for medicinal chemistry. To date, most inhibitors of viral replication have not been correlated to their respective biological targets, either in the viral or host proteome. On the other hand, the fact that CHIKV presents a biological target so similar to the pandemic SARS-Cov2 and other viruses makes the macrodomain previous studies an important source of evidence for structural dynamic understanding and selection of putative broad-spectrum inhibitors.