Energetics of Preferential Binding of Retinoic Acid-Inducible Gene-I to Double-Stranded Viral RNAs with 5' Tri-/Diphosphate over 5' Monophosphate.

Retinoic acid-inducible gene-I (RIG-I) is a cytosolic sensor protein that recognizes viral RNAs and triggers an innate immune response in cells. Panhandle-like base-paired blunt-ended 5' ppp/pp-dsRNA is a characteristic feature of viral RNAs. Structural studies of RIG-I C-terminal domain bound 5' ppp/pp-dsRNA complexes show the direct interaction between all the 5' terminal phosphates (α, β, and γ) and protein, suggesting γ phosphate might be a major recognition determinant for RIG-I binding. Biochemical studies, however, suggest that 5' pp-dsRNA is the minimal determinant for RIG-I binding and antiviral response. Despite biochemical and structural studies, the origin of viral RNA recognition by RIG-I is an unsolved problem. X-ray structures of RIG-I bound dsRNA not only provide atomic insight into the interaction network but also provide sufficiently good models for computational studies. We report structure-based molecular dynamics (MD) free energy calculations to quantitatively estimate the energetics of RIG-I binding to dsRNA containing 5' ppp, 5' pp, and 5' p. The results suggest that RIG-I weakly discriminates between 5' ppp-dsRNA and 5' pp-dsRNA (favoring former) and strongly disfavors 5' p-dsRNA with respect to the rest. Interestingly, direct interaction between γ phosphate of 5' ppp-dsRNA and RIG-I is a robust feature of the MD simulations. dsRNA binding to RIG-I is associated with Mg2+ dissociation from the 5' phosphate/s of dsRNA. The higher Mg2+ dissociation penalty from 5' ppp-dsRNA with respect to 5' pp-dsRNA offsets most of the favorable interaction between RIG-I and γ phosphate of 5' ppp-dsRNA. This leads to weak discrimination between 5' ppp-dsRNA and 5' pp-dsRNA. 5' p-dsRNA is discriminated strongly because of the loss of interaction with RIG-I.

Introduction

Retinoic acid-inducible gene-I (RIG-I) is a major pathogen recognition receptor that recognizes a broad range of viruses (e.g., influenza, rabies, dengue, hepatitis C, etc.) and triggers an antiviral response in the cytoplasm.[1−6] Panhandle-like blunt-ended dsRNA with 5′ triphosphate is a signature of viral RNA, which is preferentially recognized by RIG-I.[7−9] RIG-I is composed of a C-terminal domain (CTD), helicase domains (HEL1, HEL2, and HEL2i), Pincer domain (P), and N-terminal caspase activation and recruitment domains (CARDs) (see Figure a, CARD not shown). The CTD of RIG-I is a beta sheet bundle stabilized by Zn2+ ion (Figure a), which recognizes the viral RNA by interacting with the 5′ tri-/diphosphates of dsRNA.[10−12] Helicase core domains (HEL1 and HEL2) form an ATP binding cleft and along with HEL2i domain bind to the backbone of dsRNA[13−16] (Figure a). The Pincer domain (P) acts as a transducer of information by connecting CTD and ATPase core formed by HEL1 and HEL2 domain (Figure a).[16] N-terminal CARD domains are responsible for downstream signaling in response to viral RNA binding.[6,14] The mechanism of downstream signaling is not well-understood; however, it has been proposed[15,16] that viral RNA binding by CTD and helicase domains creates an ATP binding cleft.[17] In response to ATP binding and/or hydrolysis, the RIG-I helicase domain undergoes compaction that might potentially eject the CARDs for interaction with other proteins and subsequent signaling.[17]

Figure 1

(a) Structural overview of blunt-ended 5′ ppp-dsRNA bound RIG-I (containing helicase and CTDs) complex (PDB: 4AY2).[19] Individual domains of RIG are color-coded: CTD: cyan, HEL1: dark green, HEL2: light green, HEL2i: light blue, Pincer (P): orange. Loop spanning position 847–853 of CTD is shown as thick tube and highlighted with an arrow. Zn2+ of the zinc-binding domain of CTD is represented as a cyan sphere and ADP-Mg2+ bound to the HEL-1 domain is in the magenta sticks. Blunt-ended dsRNA is shown in yellow and the guanine base at the 5′ end with triphosphate is shown as sticks. Same color coding is used throughout this paper. The spherical region selected for MD simulation is indicated in the black broken circle. (b) Thermodynamic cycle (TC) for dsRNA:RIG-I binding. Vertical legs correspond to binding; horizontal legs correspond to the alchemical transformation of 5′ ppp-dsRNA into 5′ pp-dsRNA, either in the solvated RIG-I complex (above) or in the free state in water (below). The 5′ ppp-dsRNA/5′ pp-dsRNA binding free energy difference is ΔΔG = ΔGcomp – ΔGfree = ΔGbind(5′ pp-dsRNA) – ΔGbind(5′ ppp-dsRNA). Horizontal legs of the thermodynamic cycle (ΔGcomp, ΔGfree) were computed by MD free energy simulations. Mg2+ bound dsRNA is considered as the prevalent species in the lower horizontal leg; Mg2+ is shown as a green sphere. X-ray structure of the RIG-I:5′ppp-dsRNA complex: (c) representative structure (PDB: 3LRR)[11] revealed that RIG-I CTD of dsRNA interacts with α, β, and γ phosphates of 5′ adenine terminal of dsRNA through extensive electrostatic interactions. (d) PDB: 4AY2(19) suggests no direct interaction between γ phosphate of 5′ ppp-dsRNA and RIG-I. The network of electrostatic interactions are represented as dashed lines. Key residues of RIG-I CTD involved in RNA binding are shown as stick models.

Structural studies[11] suggest that the structure of CTD of RIG-I is very similar in complex and in isolation, except for the loop region (residue 847–853, Figure a). The free CTD structures in solution[10,18] have revealed that the loop is very flexible. Comparison with the dsRNA bound structures[11] suggests that dsRNA binding might stabilize the specific conformation of the loop. X-ray structures (PDB 3LRR and 3LRN)[11] of RIG-I CTD bound 5′ ppp-dsRNA have revealed extensive interactions between α, β, γ-phosphates of dsRNA and RIG-I, whereas a full RIG-I(ΔCARDs) bound 5′ ppp-dsRNA structure (PDB 4AY2) suggests[19] no direct interaction between γ-phosphate and RIG-I. The conflicting γ-phosphate–RIG-I interaction is also limited by the resolution of the X-ray structures. It can be concluded that out of three X-ray structures (PDB 3LRR, 3LRN, and 4AY2), first two structures suggest a direct interaction between γ phosphate and RIG-I. Thus, γ phosphate appears to be a major recognition factor for RIG-I binding, and RIG-I should discriminate strongly between 5′ ppp-dsRNA and 5′ pp-dsRNA, preferring the former. It is worth mentioning that none of the resolved structures of the complexes contain Mg2+ bound to the terminal 5′ tri-/diphosphates of dsRNA.

Recent biochemical studies[20] have revealed that blunt-ended dsRNA with 5′ diphosphate (5′ pp) is the minimum requirement for RIG-I binding and antiviral response. As both 5′ ppp-dsRNA and 5′ pp-dsRNA can induce antiviral response, γ phosphate is certainly not a recognition determinant. Biochemical studies further reveal that isolated RIG-I domains bind more strongly to viral RNA with respect to the full RIG-I.[21] It has been suggested[21] that covalently linked RIG-I domains reduce the overall binding affinity but increase the specificity for efficient discrimination of viral RNAs from host RNAs.

Despite the advancement of structural and biochemical studies, the atomic insight into the dynamics of these complexes is unknown and the following key questions on dsRNA binding to RIG-I remain unanswered: (a) how strongly RIG-I discriminates between 5′ ppp-dsRNA, 5′ pp-dsRNA, and 5′ p-dsRNA (i.e., relative binding affinity)? (b) What is the relationship between relative binding affinity and the 3D structures? (c) Why γ phosphate is not a major recognition factor for RIG-I binding?

Medium resolution X-ray structures[10,19] now provide sufficiently good models for structure-based computer simulations for addressing the above questions. We report structure-based molecular dynamics (MD) free energy simulations for deciphering the energetics of dsRNA binding to RIG-I, thereby linking 3D structures and energetics. Starting with X-ray[19] structure (PDB: 4AY2, full RIG-I(ΔCARDs):5′ ppp-dsRNA, resolution 2.8 Å) as our initial models, we have performed MD simulations of RIG-I:5′ ppp-dsRNA, RIG-I:5′ pp-dsRNA, and RIG-I:5′ p-dsRNA complexes. Calculations involve the change in binding affinity upon “mutation” of the 5′ terminal of dsRNA using an appropriate thermodynamic cycle (TC) shown in Figure b. The calculations quantitatively estimated the binding affinity difference between 5′ ppp-dsRNA, 5′ pp-dsRNA, and 5′ p-dsRNA to RIG-I binding. Our estimated binding free energy difference between 5′ ppp-dsRNA/5′ pp-dsRNA is small (∼2 kcal/mol in favor of 5′ ppp-dsRNA) and 5′ ppp-dsRNA/5′ p-dsRNA is large (∼9 kcal/mol favoring 5′ ppp-dsRNA). The signs are consistent with the experimental observations and the magnitudes seem biochemically sensible. dsRNA binding to RIG-I is associated with Mg2+ dissociation from 5′ terminal of dsRNA, desolvation, and protein–RNA interactions. Our simulations suggest that γ phosphate establish direct electrostatic contact with RIG-I. However, higher Mg2+ dissociation penalty from 5′ ppp-dsRNA with respect to 5′ pp-dsRNA offsets most of the γ phosphate-RIG-I interaction, leading to the loss of binding affinity for RIG-I. This seems to be responsible for weak discrimination between 5′ ppp-dsRNA and 5′ pp-dsRNA, with RIG-I weakly preferring the former. The binding of 5′ p-dsRNA is strongly disfavored because of the substantial loss of interaction between 5′ terminal of dsRNA and RIG-I.

Materials and Methods

MD Setup

The MD setup is given in Figure S1. Structure of human RIG-I bound to 5′ ppp-dsRNA was taken from the Protein Data Bank (entry 4AY2;[19] crystallographic resolution 2.8 Å). We have selected 4AY2 as a template for MD for two reasons: (1) the panhandle-like bound RNA with 5′ triphosphate is an excellent mimic of viral RNAs and (2) the RIG-I in this complex contains CTD as well as the helicase domains. A spherical region of radius 30 Å, centered at the terminal phosphate of dsRNA, was cut from the selected structure and considered for MD simulations. We retained residues that had at least one nonhydrogen atom within the 30 Å sphere. Nonhydrogen atoms of RIG-I and dsRNA in the outer region between 27 and 30 Å from the sphere’s center (“buffer region”) were harmonically restrained to their experimentally determined positions. The restraints in the buffer region were increased gradually from 3.0 to 5.0 kcal/mol/Å2 as one moves closer to the outer boundary, leaving the inner 27 Å radius shell fully flexible. A cubic water box (edge length = 80 Å) was overlaid for solvation, and waters that overlapped with RIG-I/dsRNA were removed. We deleted the terminal phosphate/s from 4AY2 and considered the resulting complex as the initial model of RIG-I bound 5′ pp-dsRNA/5′ p-dsRNA. The MD structures were compared with the X-ray structures of RIG-I CTD bound 5′ pp-dsRNA and found them to be essentially the same. The total number of atoms in our simulation model is about ∼50 100. The number of water molecules present in our MD is about 15 000. Root-mean-square deviation (rmsd) of the heavy atoms (within 27 Å of simulation sphere) of the complex with respect to the X-ray structure (PDB 4AY2) is given in Figure S2. Average rmsd was calculated for the heavy atoms within 27 Å of the simulation sphere, averaging over the 5 ns MD trajectory with 2 ps interval. The rmsd data suggested that the MD structures are very similar to the X-ray structure 4AY2. The largest rmsd for 5′ p-dsRNA bound complex is expected as the models are generated by deleting the β and γ phosphates from 5′ ppp-dsRNA bound RIG-I (PDB: 4AY2).

Periodic boundary conditions were used to run MD, using the Particle mesh Ewald method[22] for long-range electrostatics, with tinfoil boundary conditions.[23,24] A cutoff distance of 16 Å was used to truncate van der Waals interaction. Temperature and pressure were maintained at 310 K and 1 bar, respectively. The temperature was controlled by using Langevin dynamics for nonhydrogen atoms with a coupling coefficient of 5 ps–1, and the pressure was controlled by Langevin piston using the Nose–Hoover method. The CHARMM36 force field[25,26] with the TIP3P water model was used.[27] Simulations were done with CHARMM[28,29] and NAMD[30] programs. We performed 275–300 ns of production dynamics for each simulation model. The simulations involve 250–300 ps of equilibration through a series of short MD runs. In the first 40 ps of equilibration, the system was heated up to 310 K and then kept fixed throughout the production MD trajectories. At the initial stage of equilibration, heavy atoms of the inner region (within 27 Å) were harmonically restrained at their experimentally resolved position with a force constant of 4.0 kcal/mol/Å2, and at the final stage of equilibration, the restraint from the inner region was completely removed. The overall charge of the simulation was neutralized by scaling down the partial charges of the phosphate backbone of dsRNA. It should be noted that 5′ ppp-dsRNA → 5′ pp-dsRNA and 5′ pp-dsRNA → 5′ p-dsRNA transformation alters the overall charge of the complex. The use of tinfoil boundary conditions implemented in NAMD[30] ensures that, as the overall charge changes, a compensating charge density is spread uniformly throughout the simulation box and does not contribute to the forces.[31,32]

Protocol for Binding Free Energy Calculation

Relative binding free energies (ΔΔG) of 5′ ppp-dsRNA/5′ pp-dsRNA/5′ p-dsRNA binding to RIG-I were calculated by alchemically transforming a terminal phosphate into a ghost following the horizontal legs of the TC described in Figure b. The vertical legs correspond to RNA binding. On the other hand, horizontal legs correspond to the alchemical transformation of dsRNA which cannot be realized experimentally. We computed the free energy change associated with the phosphate deletion (horizontal arms of Figure b) and calculated the relative binding free energy as ΔΔGbind = ΔGcomp – ΔGfree = ΔGbind(5′ pp-dsRNA) – ΔGbind(5′ ppp-dsRNA). A hybrid energy function[33−36] (U) was used to represent a mixture of two endpoint states for a particular horizontal leg (Figure b). Coupling coordinate λ was used to connect two endpoints. Two coupling coordinates, λelec, and λvdW, were used to modify the electrostatic and van der Waals energy terms, where U = U(λelec, λvdW) = U1 + U2. First MD simulations gradually remove the atomic charges of the γ-phosphate of 5′ ppp-dsRNA with simultaneous modification of the β-phosphate charges to model 5′ pp-dsRNA according to the following equationwhere Uelec(5′ ppp-dsRNA) and Uelec(5′ pp-dsRNA) represent the Coulomb interactions involving 5′ ppp-dsRNA and 5′ pp-dsRNA charges, respectively. Next, we varied the van der Waals interaction parameter of γ-phosphate of 5′ ppp-dsRNA into those of 5′ pp-dsRNA by modifying the coupling coordinate λvdW from 1 to 0, and leaving a ghost phosphate, usingFree energy derivative from Boltzmann statistics was calculated as ∂G/∂λ = ⟨∂U/∂λ⟩λ, where λ = λelec or λvdW and the brackets “⟨ ⟩” represent averaging over MD trajectory for a particular value of λ. For the alchemical transformation of the charges, we used 11 equally spaced λelec values between 1 and 0 (1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, and 0.0). The free energy derivative ⟨∂U/∂λ⟩λ at each window was computed from a finite-difference estimate. Similarly, the van der Waals interactions were also transformed; the successive λvdW values were 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01, and 0.001. Each λ window simulation lasted for 1–2 nanoseconds, and the data from last 600 to 1600 ps of each simulation were used for averaging. Free energy change was calculated using numerical integration method. The standard trapezoidal method for both the complete electrostatic stage and the van der Waals stage, down to a λvdW of 0.05, was used to integrate the derivatives. van der Waals free energy derivative between λvdW = 0.05 to 0, was fitted to the function A0 λvdW– and then analytically integrated, where A0 and A1 are adjustable parameters. Uncertainties of the free energy derivatives at each λ value were estimated by dividing the trajectory segment (which was used for averaging) into two batches and taking the deviation of the batch averages. The same is reported in the manuscript as statistical error. Multiple runs were performed for each state and the average results obtained from those different runs are reported in the main text of the manuscript. Each free energy calculation was based on 24–48 ns of data collection averaged over 6–12 replicas with different initial velocities (Table S2). Overall, a total of about 1.6 μs of MD free energy simulations has been done to get good convergence and reasonable statistical error (1–2 kcal/mol), comparable to the earlier reported force field uncertainty.[36] The λ versus free energy derivatives are shown in Figure S3. Different runs are in excellent agreement with each other (Figure S3 and Table S3). The uncertainty of overall free energy change (Table S2) and free energy derivatives at each of the λ points (Table S3) is well within the acceptable statistical uncertainty.

Results and Discussion

Biochemical and Structural Views of RIG-I Binding to dsRNA

RIG-I recognizes 5′ ppp-dsRNA or 5′ pp-dsRNA and discriminates between host and viral RNA. Biochemical studies have showed that RIG-I has similar affinity for both 5′ ppp-dsRNA and 5′ pp-dsRNA, but binding is the weakest with 5′ p-dsRNA.[20] Earlier studies indicated that 5′ p-dsRNA fails to induce an antiviral response, which suggests that ligand binding is necessary but not sufficient for RIG-I activation.[11,21] Structural studies[11,12,18,19] have revealed that 5′ terminal of dsRNA is recognized by CTD of RIG-I, primarily by forming a network of electrostatic interactions. Direct interaction between α and β phosphates of dsRNA with Lys and His residues of RIG-I has been confirmed in all the resolved structures[11,12,19] (Figure c,d). Conflicting literature, however, exists regarding the interaction between γ phosphate of dsRNA and RIG-I. While crystal structures of RIG-I CTD bound 5′ ppp-dsRNA complexes [PDB codes: 3LRN (resolution 2.6 Å),[11]3LRR (resolution 2.15 Å)[11]] reveal extensive electrostatic contact between all three phosphates of dsRNA with RIG-I (Figure c), another X-ray structure (PDB: 4AY2 (resolution 2.8 Å)][19] of full RIG-I (without CARD domain) bound 5′ ppp-dsRNA suggests no interaction between γ phosphate and RIG-I (Figure d), suggesting γ phosphate may not be a major determinant for viral RNA detection. Still, the possibility of interaction between γ phosphate and RIG-I has not been ruled out.[19] The controversy related to γ phosphate interaction with RIG-I is based on the slightly different orientation of γ phosphate in the complex (Figure c,d).

We aligned some of the high resolution X-ray structures of RIG-I bound 5′ ppp/pp-dsRNA, and the results are given in Table S1. The CTD of RIG-I in all the structures is very similar within an rmsd of 0.8–1.1 Å, including the loop region. Comparison of RIG-I CTD in the free and 5′ ppp/pp-dsRNA bound structures (Table S1) does show noticeable conformational changes for the loop region (rmsd 3.2–9.9 Å). The loop region of unbound RIG-I CTD (dsRNA-free) is known to be highly flexible in solution,[10] and dsRNA binding might stabilize specific conformation of the loop. 5′ p-dsRNA bound RIG-I structure has not yet been resolved experimentally. It is worth mentioning that none of the resolved structures of the RIG-I:dsRNA report the presence of divalent metal ion, for example, Mg2+ bound to the terminal phosphate/s of dsRNA. It should be noted that water molecules deposited in the PDB files may result from misinterpretation of the electron-density maps of Mg2+. Thus, we analyzed the water molecules present in the X-ray structures 3LRR(11) (RIG-I CTD:5′ ppp-dsRNA) and 3NCU(12) (RIG-I CTD:5′ pp-dsRNA). Interestingly a water molecule (with an average beta factor > 45 Å2) is present in close proximity to an oxygen of the α phosphate (2.7 Å in 3LRN, 3.4 Å in 3NCU) in both the structures. The large beta factor of the reported water and its monodentate coordination with the α phosphate suggest that the binding site is most likely to be Mg2+-free. It is worth mentioning that both 4AY2(19) and 3NCU(12) were crystallized in the buffer containing MgCl2. Mg2+ has been isolated in the ATPase domain (Figure a) along with ADP in 4AY2;[19] however, none of the reported structures resolved Mg2+ in the 5′ terminal of dsRNA in its RIG-I bound complex.

Mutagenesis of key residues involved in interaction with the 5′-phosphates has shown to affect RNA binding and signaling by RIG-I.[11] Structural studies[11,12,18,19] further revealed that the CTD of RIG-I primarily interacts with the 5′ terminal and the backbone of the dsRNA, demonstrating that RIG-I can bind to dsRNA in a sequence-independent manner.

dsRNA with 5′ ppp, 5′ pp, and 5′ p in Solution

Consideration of unbound state is essential for understanding the dsRNA binding to RIG-I. Experiments[37] suggest that ATP4– prefers to exist mainly in [ATP:Mg]2– form with fully deprotonated γ-phosphate along with trace amounts of [ATPH:Mg]−, [ATP:Mg2], and [ATP]4–. The same is also true for ADP3– and AMP2–. We may assume the same situation, that is, Mg2+ bound dsRNA is the prevalent species when free in water. The X-ray structures of RIG-I:dsRNA, however, have ruled out the possibility of Mg2+ in the complex. This suggests that binding of 5′ ppp/pp-dsRNA involves Mg2+ dissociation from dsRNA. It is worth mentioning that the calculations of ΔGcomp (without bound Mg2+) and ΔGfree (with bound Mg2+) may introduce force field errors in relative binding free energies (ΔΔG’s, Figure b). Alchemical transformation of 5′ ppp-dsRNA to 5′ pp-dsRNA with bound Mg2+ involves the removal of a single Mg2+:phosphate interaction. Mg2+ polarizes the electronic cloud,[36] and the fixed charge force fields are inaccurate as polarization effects are only implicitly included. The effect of electronic polarization by Mg2+ is known[36] to cancel out substantially if Mg2+ is present in both the horizontal legs of the TC. In this study, however, force field error might be significant as the Mg2+ bound 5′ ppp/pp/p-dsRNA is expected to be the prevalent form in the RIG-I free form. Thus, we first computed relative binding of Mg2+ to dsRNA using an appropriate TC (Figure ) and compared with experimental binding data.[38] The results are summarized in Table . The results clearly indicate that the calculated relative binding free energies (column 4 of Table ) are larger compared to that of free adenosine nucleotide in water (column 5 of Table ). Experimentally reported binding constants[38] for Mg2+ binding to ATP and ADP (ADP and AMP) gives a standard binding free energy difference of +2.1 kcal/mol favoring ATP (+2.5 kcal/mol favoring ADP). These binding free energies are expected to be close to the calculated ΔΔG’s (column 4, Table ). Thus, the force field gives large error by overestimating Mg2+:phosphate binding because of the lack of explicit electronic polarizability.[39−41] Fixed charge force field is also known to make errors in accurately capturing GTP and GDP binding when electrostatic interactions with a divalent ion are involved.[36] Therefore, to compute the ΔGfree (Figure b), we have taken the experimental Mg2+ binding data (Table ), such that ΔGfree = ΔG2 + ΔΔGexpt. The MD simulation show that three water coordinated Mg2+ forms tridentate coordination with β,γ phosphates of 5′ ppp-dsRNA and α,β phosphates of 5′pp-dsRNA (Figure S4). Mg2+ forms bidentate coordination with the phosphate group of 5′ p-dsRNA, and the coordination sphere of Mg2+ is completed by coordination with four water molecules (Figure S4).

Figure 2

Thermodynamic cycle for Mg2+:dsRNA binding. Vertical legs correspond to binding. Horizontal legs correspond to the alchemical transformation of 5′ ppp-dsRNA to 5′ pp-dsRNA in the presence (upper horizontal leg) and absence (lower horizontal leg) of Mg2+. Free energy change for the horizontal legs is computed using MD free energy simulations. Experimentally reported standard Mg2+ binding free energies to ATP and ADP have been considered (expected to be very close to Mg2+ binding to 5′ ppp-dsRNA/5′ pp-dsRNA) to quantify the force field error. ΔG1 – ΔG2 = ΔGbind[5′ pp-dsRNA:Mg2+] – ΔGbind[5′ ppp-dsRNA:Mg2+] ≈ ΔGbind[ADP:Mg2+] – ΔGbind[ATP:Mg2+].

Table 1

Relative Binding Free Energies of Mg2+ Binding to dsRNA (with 5′ ppp/5′ pp/5′ p)a

alchemical transformation	ΔG1	ΔG2	ΔΔGcalc = ΔG1 – ΔG2	ΔΔGexpt	ΔGfree = ΔG2 + ΔΔGexpt
5′ ppp-dsRNA → 5′ pp-dsRNA	440.2 (1.7)	436.6 (1.6)	3.6 (2.3)	ΔGbind0[ADP:Mg]− – ΔGbind0[ATP:Mg]2– = 2.1	438.7 (1.6)
5′ pp-dsRNA → 5′ p-dsRNA	397.6 (1.4)	387.7 (1.2)	9.9 (1.8)	ΔGbind0[AMP:Mg] – ΔGbind0[ADP:Mg]− = 2.5	390.2 (1.2)

Free energies are in kcal/mol. Statistical uncertainty is given in the parentheses. The MD trajectories were divided into two equal halves and the difference between the computed ΔG’s from the two-halves is reported as uncertainty in the parenthesis. The uncertainties for ΔΔG’s were calculated by propagating the uncertainties of individual ΔG’s.

Table 2

Relative Binding Free Energies: 5′ ppp-dsRNA/5′ pp-dsRNA/5′ p-dsRNA Binding to RIG-Ia

alchemical transformation	ΔGcomp	ΔGfree	ΔΔG
5′ ppp-dsRNA → 5′ pp-dsRNA	441.0 (1.8)	438.7 (1.6)	2.3 ± 2.4
5′ pp-dsRNA → 5′p-dsRNA	397.4 (2.0)	390.2 (1.2)	7.2 ± 2.3
5′ ppp-dsRNA → 5′ p-dsRNA			9.5 ± 3.3

Free energies are in kcal/mol. Uncertainties are calculated in the same way described in Table .

Structure-Based Energetics of dsRNA Binding to RIG-I

To compute the dsRNA discrimination by RIG-I and elucidate the effect of different dsRNAs (i.e., 5′-ppp, 5′-pp and 5′-p) on the selectivity, we carried out extensive MD free energy (MDFE) simulations of RIG-I:dsRNA complexes (5′-ppp or 5′-pp or 5′-p) using the X-ray structure (PDB 4AY2)[19] as the template. We computed the change in binding affinity for dsRNA upon 5′ ppp → 5′ pp → 5′ p mutations in the RIG-I bound complex. Relative binding free energies are summarized in Table . The results suggest that RIG-I imposes a very high energetic penalty of about 9 kcal/mol for binding 5′ p-dsRNA with respect to 5′ ppp-dsRNA, which corresponds to a probability of 10–7 relative to binding 5′ ppp-dsRNA. Interestingly RIG-I selectivity between 5′ ppp-dsRNA and 5′ pp-dsRNA was predicted to be weak with a binding free energy difference of 2 kcal/mol, favoring 5′ ppp-dsRNA. This corresponds to a read-through frequency as high as 10–2 for dsRNA(5′-pp) with respect to dsRNA(5′-ppp). Biochemical studies indeed suggested RIG-I could be activated by 5′ ppp-dsRNA or 5′ pp-dsRNA but not by 5′ p-dsRNA.[20]

Comparison between MD and X-ray Structures

The rmsd of the heavy atoms of about 1.3–1.8 Å with respect to the X-ray structure suggests that the simulated structures agree well with the X-ray structure (Figure S2). Robust features of our MD simulations are as follows:

γ-Phosphate of 5′ ppp-dsRNA forms direct electrostatic interaction with the positively charged lysine side chains of RIG-I (Figure a). The MD structures of the binding pocket are almost identical to their corresponding X-ray structures (Figure and Table ).

Figure 3

Binding pocket: MD snapshot (cyan) compared to the crystal structure (gray; PDB code: 4AY2 for 5′ppp-dsRNA:RIG-I, PDB code: 3NCU for 5′pp-dsRNA:RIG-I). Hydrogens are omitted for clarity. (a) 5′ ppp-dsRNA:RIG-I. (b) 5′ pp-dsRNA:RIG-I. (c) 5′ p-dsRNA:RIG-I. MD structures are very similar to the X-ray structure, except the γ phosphate of dsRNA which forms direct interaction with Lys858 and Lys849 throughout the MD trajectory (shown in “a”). Aromatic loop residue Phe853 stacks on the terminal blunt-ended base pair and show considerable movement parallel to G1-C20 base pair (indicated by a double-headed arrow). Key residues involved in RNA terminal recognition are shown as stick models and loop (residue 847–853) is shown as a cartoon. Color code same as Figure .

Table 3

Selected Interatomic Distances Averaged over the MD Trajectoriesa

		5′ ppp-dsRNA complex			5′ pp-dsRNA complex		5′ p-dsRNA complex
		X-ray			X-ray
complex interacting pair		4AY2	3LRN	MD	3NCU	MD	MD
Lys849:NZ	RNA:02G	8.30	×	2.95 (0.67)	×	×	×
Lys849:NZ	RNA:03G	×	8.60	3.10 (0.20)	×	×	×
Lys851:NZ	RNA:03G	5.70	×	×	×	×	×
Lys851:NZ	RNA:01G	×	×	2.70 (0.23)	×	×	×
Lys851:NZ	RNA:02G	×	3.70	3.03 (0.40)	×	×	×
Lys858:NZ	RNA:02B	2.70	2.60	2.68 (0.10)	2.90	2.63(0.09)	×
Lys858:NZ	RNA:O1G	×	3.00	2.60 (0.40)	×	×	×
Lys858:NZ	RNA:N7	3.00	3.00	2.89 (0.18)	3.20	2.92 (0.13)	3.69(0.71)
His847:NE2	RNA:02B	2.90	2.50	2.75 (0.21)	2.80	2.73 (0.19)	×
Lys861:NZ	RNA:01B	3.00	2.90	2.77 (0.13)	×	2.62 (0.32)	×
Lys861:NZ	RNA:03B	×	×	×	2.90	×	×
Lys861:NZ	RNA:02A	2.60	2.30	2.72 (0.10)	×	2.75 (0.17)	×
Lys861:NZ	RNA:O1A	×	×	×	2.90	×	×
Lys888:NZ	RNA:O1A	2.90	3.30	2.67 (0.10)	3.10	2.96 (0.66)	4.03(0.69)
Lys888:NZ	RNA:02A	3.30	2.50	4.06 (0.30)	3.20	3.28 (0.49)	3.40(0.38)
His830:ND1	RNA:02′	2.70	2.40	2.80 (0.18)	2.70	2.80 (0.12)	2.82(0.19)

Standard deviations are in the parentheses. Distances are in angstrom. Residues/groups absent in the X-ray/MD structures are indicated with a cross.

The loop region (residue 847–853) is highly flexible (Figure a); however, the interactions between loop residues and 5′ terminal phosphates of dsRNA remain intact throughout the simulations (Figure ). Phe853 stabilizes the complex by hydrophobic interaction with either G1 or C20, moving parallel to the 5′ terminal base pair G1-C20.

Figure 4

RIG-I:5′ ppp-dsRNA complex from MD: (a) flexibility of the loop region (residue 847–853): overlaid 25 snapshots with a 200 ps spacing from a 5 ns MD trajectory. The loop is highly flexible for RIG-I complexes with 5′ pp-dsRNA and 5′ p-dsRNA (Figure S5). Color code same as Figure . (b) Water density around the terminal γ phosphate of 5′ ppp-dsRNA. Results are shown for the RIG-I bound complex (thick line) and free dsRNA with bound Mg2+ (dash line). The water numbers are averaged in radial shells over the MD trajectories. The same trend is observed for 5′ pp-dsRNA and 5′ p-dsRNA (Figure S6).

5′ ppp/pp/p-dsRNA binding to RIG-I is associated with desolvation of 5′ terminal of dsRNA (Figure b).

The above structural features are insensitive to the initial structural models used in these MD simulations (see Materials and Methods).

Structure and Dynamics of the RIG-I:5′ ppp-dsRNA Complex

The structures from MD simulations agree well with the corresponding 4AY2 crystal structure (Figure a and Table ). The rmsd of main chain and side chain heavy atoms is 0.82 ± 0.1 and 1.68 ± 0.15 Å, respectively, with respect to 4AY2. The 5′ triphosphate binding pocket is located at the positively charged patch of RIG-I (Figure a), and the structural details are summarized in Table . Five lysines, two histidines, and one phenylalanine form the binding pocket for 5′ triphosphate (Figure a). Direct interaction between α phosphate and Lys888 was observed with a mean O–NZ distance of 2.67 ± 0.1 Å. Lys861 interacts simultaneously with α and β phosphate with a mean O–NZ distance of 2.72 ± 0.1 and 2.77 ± 0.13 Å, respectively. β phosphate is further stabilized by interacting with the side chains of Lys858 (mean O–NZ distance of 2.68 ± 0.1 Å) and His847 (mean O–NE2 distance of 2.75 ± 0.21 Å). Lys858 interacts with 5′ G1 with an average distance of 2.89 ± 0.18 Å. Deviation of γ phosphate of 5′ ppp-dsRNA with respect to the X-ray structure is shown in Figure a; it is forming direct electrostatic contact with the side chains of Lys849, Lys 851, and Lys858. γ phosphate of dsRNA forms bidentate coordination with Lys849 (an average O–NZ distance of 3.1 ± 0.2 and 2.95 ± 0.67 Å) and Lys851 (an average O–NZ distance of 2.7 ± 0.23 and 3.03 ± 0.4 Å), and a monodentate coordination with Lys858 (an average O–NZ distance of 2.6 ± 0.4 Å). Phe853 stacks over the terminal G1-C20 base pair, and His830 forms hydrogen bond with the ribose −OH of G1 with an average distance of 2.8 ± 0.18 Å. MD trajectory shows the movement of Phe853 parallel to the G1-C20 base pair and stacking with both G1 and C20 bases (Figures a and S5). The water density around the γ phosphates of free and RIG-I bound dsRNA is shown in Figure b. On average, the γ phosphate of 5′ ppp-dsRNA is solvated with six and eight water molecules in RIG-I bound and free forms, respectively. Thus, in the binding pocket, the availability of water is reduced by two water molecules.

Structure and Dynamics of the RIG-I:5′ pp-dsRNA Complex

We deleted the γ phosphate of the RIG-I:5′ ppp-dsRNA complex (4AY2) and considered that as our model for 5′ pp-dsRNA bound RIG-I. In our model, the helicase domain is present which is absent in the X-ray structure of 5′ pp-dsRNA bound RIG-I complex (PDB 3NCU).[12] MD structures show that G1-RIG-I interactions are almost identical to the crystal structure (PDB: 3NCU).[12] Structural details are summarized in Figure b and Table . With respect to the 4AY2, the deviations are 1.06 ± 0.14 Å for the main chain and 1.76 ± 0.08 Å for side chains. Direct interactions of α and β phosphates of dsRNA with Lys888, Lys861, His847, and Lys858 have been observed (with an average distance ranging between 2.62 and 3.28 Å; standard deviation of less than 0.7 Å). The overall charge of −3 of the terminal G1 is neutralized by side chains of three lysine residues. His830 and Lys858 form hydrogen bonds with the ribose −OH and N7 of G1 with an average distance of 2.8 ± 0.12 and 2.92 ± 0.13 Å, respectively. Movement of Phe853 parallel to the 5′ terminal base pair G1-C20 has also been observed in the MD trajectory. Desolvation of 5′ pp-dsRNA (similar to 5′ ppp-dsRNA) has been observed upon binding to RIG-I (Figure S6).

Structure and Dynamics of the RIG-I:5′ p-dsRNA Complex

The RIG-I:5′ p-dsRNA structure is not known experimentally. We have deleted β and γ phosphates from 4AY2 and considered that as the MD model. The average rms deviations for main chain and side chain heavy atoms are 1.1 ± 0.28 and 1.97 ± 0.16 Å, respectively, with respect to the starting MD model. Structural parameters of the RIG-I:5′ p-dsRNA complex are given in Figure c and Table . The absence of β and γ phosphates leads to the loss of interaction with RIG-I. The only interaction between His830 and ribose −OH of G1 remains intact as seen in RIG-I:5′ ppp/pp-dsRNA complexes. The average distance of 3.4 ± 0.38 Å between Lys888 and α phosphate of 5′ p-dsRNA suggests weak interaction. Similar loop flexibility (Figure S5) as seen in RIG-I bound 5′ ppp/pp-dsRNA has been observed. Binding of 5′ p-dsRNA to RIG-I is associated with desolvation (Figure S6) as also seen in 5′ ppp/pp-dsRNA. MD suggests that binding of 5′ p-dsRNA to RIG-I is unfavorable because of the loss of direct interactions between 5′ p-dsRNA and RIG-I.

Understanding the RIG-I:dsRNA binding is a difficult challenge. Multiple RIG-I conformations (bound/unbound states) and multiple dsRNA species (e.g., different protonation states of dsRNA with bound/unbound ions and its RIG-I bound/unbound conformations) contribute to the overall binding process. Despite recent biochemical studies and medium resolution structures, we are far away from understanding the structure-based detailed energy landscape associated with viral RNA recognition. Measuring the binding free energy between a specific nucleotide and a particular protein conformation is an extremely difficult task. MD simulations can only fill the gap to some extent. Using MDFE simulations, we calculated 5′ ppp-dsRNA/5′ pp-dsRNA/5′ p-dsRNA binding free energy differences. The calculated strength of discrimination, ΔΔG ≈ 9 kcal/mol, suggests that 5′ ppp-dsRNA binding to RIG-I is strongly favored with respect to 5′ p-dsRNA. On the other hand, small ΔΔG ≈ 2 kcal/mol suggests that 5′ ppp-dsRNA binding to RIG-I is weakly favored with respect to 5′ pp-dsRNA. Because the magnitude of the relative preferences (ΔΔG) is not known experimentally, the calculated MDFE values cannot be confirmed or disproved. Certainly, the signs are correct and corroborates the experiment.[20]

X-ray and NMR structures suggest that the loop region (residue 847–853) of free RIG-I is highly flexible,[10] and dsRNA binding might lead to stabilization of specific conformation of the RIG-I loop.[11] The average value of B-factor of the loop region (residue 847–853) of our template PDB 4AY2 is 85 Å2, indicating considerable movement. MD simulations suggest that the loop region (residue 847–853) is flexible and can have multiple conformations (Figures a and S5) even in the dsRNA bound state. The key interaction between the loop side chains and 5′ terminal of dsRNA are, however, retained throughout the MD trajectory (Figure a,b and Table ). The direct interaction between Lys858 and N7 of 5′ terminal base of dsRNA (Figure , Table ) is very specific to purines and might disfavor pyrimidines at the 5′ end. Almost all the interaction distances in the MD structures agree with their corresponding X-ray structures (Table ), except for the fact that γ phosphate of 5′ ppp-dsRNA establishes a direct contact with the RIG-I (absent in template PDB 4AY2). MD results are in agreement with the X-ray structures of RIG-I CTD bound 5′ ppp-dsRNA, revealing direct electrostatic interactions between RIG-I and γ phosphate of dsRNA.[11,42] It might appear that γ phosphate is a crucial recognition determinant and RIG-I should strongly favor 5′ ppp-dsRNA binding with respect to 5′ pp-dsRNA. Experiments,[20] however, have clearly shown that RIG-I discriminates weakly between 5′ ppp-dsRNA and 5′ pp-dsRNA. 5′ pp-dsRNA is the minimum requirement for RIG-I binding and antiviral response. To understand dsRNA binding, consideration of unbound state is essential. Free 5′ ppp/pp-dsRNA with bound Mg2+ is expected to be the prevalent species. Structural studies suggest that 5′ ppp/pp-dsRNA in the RIG-I bound form does not have Mg2+, indicating that Mg2+ detachment might be required for dsRNA binding to RIG-I. It is evident from the structures that the negatively charged 5′ terminal of dsRNA is neutralized in the RIG-I binding pocket through extensive electrostatic interaction, justifying the absence of Mg2+ in the binding pocket. Structures of the binding pocket from MD simulations are very similar to their X-ray structures (Figure , Table , and Figure S2), suggesting that the starting structure of the complex (without Mg2+) is a true minimum in the potential energy hypersurface. MD further revealed that dsRNA binding to RIG-I is associated with desolvation. The RIG-I binding site is water-exposed, and the extent of desolvation seems to be similar for 5′ ppp/pp/p-dsRNA binding (Figures a and S6). The results suggest that desolvation is crucial for binding but may not be significant toward relative binding preference. It should be noted that the dissociation constants for the Mg2+:AMP2–, Mg2+:ADP3–, and Mg2+:ATP4– complexes are Kd = 1.622 × 10–3, 2.239 × 10–5, and 6.607 × 10–7 M at room temperature (298.15 K) and pressure (1 atm); the standard binding free energies (ΔGbind0) are −3.8, −6.3, and −8.4 kcal/mol, respectively.[38] Clearly the energetic cost for Mg2+ dissociation is the highest from ATP4– > ADP3– > AMP2– for electrostatic reason. We may expect the same is true for dsRNA containing 5′ phosphate/s; hence, interaction between dsRNA and RIG-I in the binding pocket needs to offset the Mg2+ dissociation penalty from the dsRNA. Both 5′ ppp-dsRNA and 5′ pp-dsRNA form an extensive interaction network in the RIG-I binding pocket. The specificity between 5′ ppp-dsRNA and 5′ pp-dsRNA, however, was predicted to be weak with a binding free energy difference of ∼2 kcal/mol, favoring 5′ ppp-dsRNA. Interactions between the γ phosphate of 5′ ppp-dsRNA and lysine side chains in the RIG-I binding pocket offset the largest Mg2+ dissociation penalty (ΔGbind0 = −8.4 kcal/mol); this might be the reason for weak specificity between 5′ ppp-dsRNA and 5′ pp-dsRNA. Strong RIG-I binding specificity of ∼9 kcal/mol in favor of 5′ ppp-dsRNA with respect to 5′ p-dsRNA is due to the loss of protein-RNA interaction in the latter.

The close proximity of His847 to the negatively charged β phosphates of dsRNA might stabilize the protonated state of His847. MD simulations with the RIG-I:5′ ppp-dsRNA complex (with protonated His847) suggest alternate conformation of His847 (see Figure S7), in which the interaction with β phosphate is lost and interaction with γ phosphate is formed. It should be noted that in all the X-ray structures (PDB 3LRR, 3LRN, and 4AY2),[11,19] the position of His847 is very much similar and His847 interacts only with the β phosphate. The alternate conformation seen in MD simulation with protonated His847 is not supported by X-ray structures. To the best of our estimate, RIG-I with protonated His847 imposes a higher energetic penalty of about 3.7 ± 3.2 kcal/mol (relative to 2.3 ± 2.4 kcal/mol for deprotonated His847) for binding to 5′ pp-dsRNA with respect to 5′ ppp-dsRNA. Biochemical studies[20] have shown that both 5′ ppp-dsRNA and 5′ pp-dsRNA can induce an antiviral response. Hence the magnitude of relative binding affinity ΔΔG (RIG-I binding to 5′ ppp-dsRNA vs 5′ pp-dsRNA) is expected to be small. The deviation of protonated-His847 (in MD) with respect to its experimentally determined position, along with the larger ΔΔG (computed), probably suggests the deprotonated state of His847.

Conclusions

MD simulations of RIG-I:5′ ppp-dsRNA, RIG-I:5′ pp-dsRNA, and RIG-I:5′ p-dsRNA complexes provide crucial insights into the energetics of viral RNA binding to RIG-I. Extensive interaction between 5′ terminal of dsRNA and RIG-I is favorable for binding, whereas desolvation and dissociation of Mg2+ from dsRNA are unfavorable. The extent of desolvation is very similar for 5′ ppp-dsRNA, 5′ pp-dsRNA, and 5′ p-dsRNA, but the Mg2+ binding free energy is the largest for 5′ ppp-dsRNA and the smallest for 5′ p-dsRNA. RIG-I binding discriminates weakly between 5′ ppp-dsRNA and 5′ pp-dsRNA with a binding free energy difference of ∼2 kcal/mol, favoring 5′ ppp-dsRNA. The favorable interaction between RIG-I and γ phosphate of 5′ ppp-dsRNA in the binding pocket is mostly offset by the higher Mg2+ dissociation penalty, leading to a weak discrimination between 5′ ppp-dsRNA and 5′ pp-dsRNA. 5′ p-dsRNA is discriminated strongly by ∼9 kcal/mol with respect to 5′ ppp-dsRNA due to the loss of interaction between 5′ terminal of 5′ p-dsRNA and RIG-I. In spite of all the limitations (sampling, convergence, force fields, etc.) associated with MD simulations, the signs obtained for relative binding free energies (ΔΔG) are correct since RIG-I preferentially binds 5′ ppp-dsRNA/5′ pp-dsRNA with respect to 5′ p-dsRNA. The strength of discrimination (ΔΔG) appears to be biochemically plausible. The precise relative binding free energies are not known experimentally and we hope that our study will encourage experimental verification. The calculations provide an insight into the RIG-I:5′ p-dsRNA complex, which has not been characterized experimentally. Our MD simulations are valuable in linking microscopic structures and free energies, illustrating RIG-I selectivity for blunt-ended dsRNA with 5′ ppp/5′ pp/5′ p. Binding affinity differences between 5′ ppp-dsRNA and 5′ OH-dsRNA[21] is beyond the scope of this paper.

Our preliminary data (on 5′ ppp/OH-dsRNA binding to RIG-I) calculated using the same strategy are in good agreement with the experiment,[21] and a detailed structure-based energetics study will be published in another paper. The calculations suggest that 5′ ppp-dsRNA binding is favored with respect to 5′ OH-dsRNA by about 6.5 kcal/mol, and the loss of protein–ligand interaction is mainly responsible for the relative binding strength. Furthermore, in an attempt to quantify the energetic effect due to the presence of Mg2+ in binding, we will be introducing Mg2+ into the RIG-I binding pocket. Recent findings[43,44] suggest that dsRNA and dsDNA both can induce antiviral response but the latter does not bind to RIG-I. In future, the methodology could be used to quantify dsRNA versus dsDNA binding to RIG-I.