Exploration of inositol 1,4,5-trisphosphate (IP3) regulated dynamics of N-terminal domain of IP3 receptor reveals early phase molecular events during receptor activation

Inositol 1, 4, 5-trisphosphate (IP3) binding at the N-terminus (NT) of IP3 receptor (IP3R) allosterically triggers the opening of a Ca2+-conducting pore located ~100 Å away from the IP3-binding core (IBC). However, the precise mechanism of IP3 binding and correlated domain dynamics in the NT that are central to the IP3R activation, remains unknown. Our all-atom molecular dynamics (MD) simulations recapitulate the characteristic twist motion of the suppressor domain (SD) and reveal correlated ‘clam closure’ dynamics of IBC with IP3-binding, complementing existing suggestions on IP3R activation mechanism. Our study further reveals the existence of inter-domain dynamic correlation in the NT and establishes the SD to be critical for the conformational dynamics of IBC. Also, a tripartite interaction involving Glu283-Arg54-Asp444 at the SD – IBC interface seemed critical for IP3R activation. Intriguingly, during the sub-microsecond long simulation, we observed Arg269 undergoing an SD-dependent flipping of hydrogen bonding between the first and fifth phosphate groups of IP3. This seems to play a major role in determining the IP3 binding affinity of IBC in the presence/absence of the SD. Our study thus provides atomistic details of early molecular events occurring within the NT during and following IP3 binding that lead to channel gating.

. Different sub-domains of IP 3 -bound NT of IP 3 R. The suppressor domain (SD, 1-224); IP 3 binding core-β (IBC-β, 225-436) and IBC-α (437-604) are shown as cyan, orange and green, respectively. The HS loop (residue 165-180) of SD that interacts with the IBC-β of a neighboring monomeric subunit in the functional tetrameric form of IP 3 R is colored in red. IP 3, which binds at the cleft formed between IBC-β and IBC-α is shown in stick representation.
www.nature.com/scientificreports www.nature.com/scientificreports/ also provided insight into the prime role of the SD in mediating functional coupling between ligand binding and channel gating 10 . However, a dynamic view of the structural events during and immediately following IP 3 binding at the NT that lead to the opening of a distal (located ~100 Å away) pore still remains to be elucidated.
As mentioned above, the SD plays a crucial role in traversing IP 3 -binding events at the NT to the channel pore. However, as the name ('suppressor domain') implies, the SD was found to suppress the IP 3 binding affinity of IBC 11,12 . Purified protein constructs comprising only the IBC-β and IBC-α binds IP 3 with higher affinity than the NT or even the whole protein. A cluster of 7 conserved amino acid residues located on one side of the SD were found to be critical for the suppression of IP 3 binding affinity 13 . IP 3 R isoforms are well known to differ in their IP 3 -binding affinity 3 and this has been attributed to the structural difference within their SD 11 . The Gibb's free energy (∆G) of IP 3 binding calculated from the thermodynamics analysis also reflected on the suppressing nature of SD, where ∆G of IP 3 binding at 296 K was −37.1 ± 0.1 kJ/mol and −43.2 ± 0.1 kJ/mol in presence (i.e., NT) and absence (i.e., IBC) of the SD respectively 14 . Suppression of IP 3 -binding affinity of the IBC by the SD is thus well established by now, though the exact molecular basis underlying such phenomenon remains elusive.
In the present study, we have used a computational approach, combining ab initio and homology modelling, all-atom MD simulations and principal component analyses to capture the landscape of IP 3 -mediated conformational dynamics in IP 3 R NT and to unearth the molecular details of SD-evoked differential IP 3 binding affinity of the NT. Analysis of essential dynamics revealed that the binding of IP 3 triggers a characteristic twist motion of SD, followed by 'clam closure' of IBC-β and IBC-α. Inspection of the atomic-level interaction profile led to the identification of a tripartite interaction involving Arg54,Glu283 and Asp444 at the SD-IBC interface, which potentially involved in triggering IP 3 -mediated receptor activation. Very intriguingly, hundreds of nanoseconds long MD simulations of IP 3 -bound IBC in presence and absence of the SD revealed that the ' Arg269-flipping' and its positional interaction either with 1 st (P1) or 5 th (P5) phosphate groups of IP 3 largely govern the affinity of IBC for IP 3 . To the best of our knowledge, this is the first report that has unveiled the dynamic picture of IP 3 -invoked conformational changes in IP 3 R NT and more importantly, the molecular mechanism underlying the SD-dependent differential IP 3 binding affinity of the receptor.

Results and Discussion
To elucidate the conformational dynamics associated with the activation of IP 3 R by the natural agonist IP 3 , we performed all-atom unbiased MD simulations of IP 3 R NT in its apo and IP 3 -bound states. Different systems simulated in this study are presented in Table 1 and detailed description of system preparation is provided in the methodology section. Each system was initially simulated for 300 ns of production run following 10 ns of equilibration at 300 K. The results from the simulation studies are discussed below. IP 3 driven twist motion of the suppressor domain. The root mean square deviation (RMSD) of NT with respect to the minimized starting structure as a function of simulation time shows that the IP 3 -bound NT exhibits significant structural deviation, marked by higher RMSDs (Fig. S1a). It is noteworthy that, compared to the apo form, IP 3 -bound NT fluctuates considerably throughout the 300 ns of simulation. This observation is further supported by the computed residue-wise root-mean-squared fluctuations (RMSFs), indicating increased dynamics of the NT due to the binding of IP 3 (Fig. S1b). Although RMSDs and RMSFs reflect the conformational changes at NT while IP 3 binds, it is important to understand the prominent and collective domain motions associated with IP 3 binding, which can provide insight into the IP 3 -mediated activation of IP 3 R.
We performed principal component analysis (PCA) on simulation trajectories to capture the significant and concerted motions of different NT domains, in presence and absence of IP 3 . PCA identifies the top essential modes (eigenvectors) that represent the major part of collective motions. The top two eigenvectors (principal components, PCs) have captured the significant motions of IP 3 R NT (Fig. S2). To compare the conformational space sampled by apo (system 1) and liganded NT (system 2), the two-dimensional projection of the simulation ensembles onto the plane defined by the first two PCs were plotted (Fig. S3a). Although the regions explored by the two systems overlap, notable difference in conformational sampling can be seen between the apo and IP 3 -bound states along the PC1 and PC2. A similar trend was also noted from the PCA on the second set of  www.nature.com/scientificreports www.nature.com/scientificreports/ simulations data. This observation is consistent with the results from RMSD and RMSF analyses (Fig. S1a,b), where the protein dynamics was found to be intensified in the presence of the IP 3 .
The principal motions of protein residues can be better visualized and interpreted by representing the eigenvectors as porcupine plots 15 . The porcupine plots showing the motion of SD and IBC of apo and IP 3 -bound NT, along the direction of PC1 and PC2, are presented in Fig. 2a-d respectively. In the apo state ( Fig. 2a,b), loop regions in the NT experienced more flexibility compared to rest of the protein. Particularly, only the arm loop of the SD exhibits significant fluctuations whereas the remaining secondary structures show minimal motions. Similarly, the loop regions of the IBC-β and IBC-α showed comparatively more flexibility than the other part of the domain. Notably, in the apo state, the IBC-α exhibits directional movement towards the IBC-β, where the latter experiences negligible mobility (Movie S1).
Interestingly, IP 3 binding at the IBC substantially altered the dynamics of different structural units of NT (Fig. 2c,d). The binding of IP 3 has directly influenced the dynamics of IBC. As in apo NT, IBC-α moved towards the IBC-β in the IP 3 -bound NT. Moreover, in contrast to the apo NT where IBC-β exhibits minimal displacement, the IBC-β of the IP 3 -bound NT underwent larger dynamics and moved towards the IBC-α. As shown by the vector direction in Fig. 2c,d, IP 3 induced motions in the NT brings the IBC -β and -α close to each other. This observation is consistent with the findings from earlier studies. The crystallographic experiments showed that IP 3 binding evokes 'clam closure' , where the β and α domains of IBC get close to each other 4,5,8 . Such ligand-induced domain closure is also reported in other ligand-gated ion channels that notably include the ionotropic glutamate receptors 16 . Markedly, as seen from the porcupine plots, the binding of IP 3 at IBC has allosterically influenced the dynamics of SD. The SD attains a typical twist motion, where the entire SD sweeps over the inner surface of the IBC (Fig. 2c,d, Movie S2). The SD anchored on the inner surface of the IBC through the IBC-β and the IP 3 binding causes the SD to twist along the inner surface of IBC-α. The β2/β3, β5/β6, and β10/β11 loops of SD interact with α8 and α9 helices of IBC-α at the α-interface (see Fig. 1). Such a twist motion triggers the SD to move towards the IBC-α, which in turn causes an increased translational displacement of the conserved HS loop of SD. The IP 3 -evoked change in NT dynamics is consistent with the results from RMSD and RMSF analyses, where the increase in RMSD and RMSF in system 2 can be attributed to the IBC clam closure and twist motion of SD (Fig. S1a,b). The extent of translational displacement of SD due to its twist motion in IP 3 -bound NT is further amplified by calculating the twist angle along the simulation trajectories (described in the methodology). Figure 2e shows the evolution of twist angle along the simulation time in apo and IP 3 -bound systems. It can be noted that the SD twists minimally in the apo system with the twist angle fluctuating over a range of 0°-10°, with an average of ~5°. In accordance with the observation from essential dynamics (Fig. 2c,d), the angle calculation also captured the increased twist motion of SD caused by the IP 3 binding at IBC (Fig. 2e). In contrast to the apo system, the angle underwent much larger fluctuation (2°-23°, with an average of ~13°) during simulation. Such a wide range of twist angle shows the intensive translational displacement experienced by the SD during IP 3 binding to the IBC. Recent structural studies also support the observed IP 3 driven twist motion of the SD, where superimposition of apo and IP 3 -bound NT revealed ~9° angular displacement between the SD arm helices 5 . However, in contrast to the finding from static crystal structures, MD simulations here provided a dynamic view of IP 3 mediated binding events at the IP 3 R NT. As depicted in Fig. 2e, the SD twist extensively during IP 3 binding, where the domain sampled through the conformations with lowest twist angle value of 2°, as seen in apo system, as well as explored the conformations exclusive to the IP 3 -bound state, where the SD shows highest degree of twist (~23°). The change in buried interfacial surface area at SD -IBC interface, illustrated in Fig. 2f, also support the observed twist motion of SD. The computed average buried interfacial surface area is 1999 ± 232 Å 2 in the apo system whereas the buried area reduced to 1239 ± 134 Å 2 in the IP 3 -bound system. Such a reduction in SD -IBC interfacial buried area in IP 3 -bound system implies that the interface has weakened during the IP 3 binding and many of the residues at the interface are exposed to the solvent. This can be attributed to the enhanced twist motion of SD that affected the SD -IBC interfacial packing. Taking together, the analyses of apo and IP 3 -bound systems clearly indicate that the IP 3 binding at IBC allosterically influences the SD dynamics and the latter experiences a characteristic twist motion.
It is intriguing to note that kinetic analyses of single IP 3 Rs activity recorded using the nuclear patch-clamp technique have shown rapid intra-burst flickers that, unlike the much longer inter-burst intervals, remain largely insensitive to the variation of concentrations of IP 3 as well as other cytosolic modulators (e.g. free Ca 2+ or ATP) 17,18 . Such stereotypical, ligand-insensitive flickers were also observed for a mutant rat IP 3 R1 reconstituted in planar lipid bilayers 19 . The flickers were also noticeable, albeit with much lower frequency, in spontaneously active IP 3 Rs recorded from the Xenopus oocyte at ultra-low (<5 nM) level of cytosolic free Ca 2+ concentration 20 . Based on the PCA and twist angle analysis (Fig. 2), we hypothesise that the basal twist motions observed in the SD in our simulation with the apo NT relate to conformations that probably underlie the ligand-insensitive intra-burst flickers of IP 3 Rs, which, however, requires further experimental validation.

IP 3 binding breaks the triad of interaction at SD -IBC interface. In the earlier section, we have
shown that the SD sweeps over the inner surface of IBC at the SD -IBC interface when IP 3 binds at the IBC pocket. As observed in earlier structural studies 4,5,13,21 , our simulations also identified crucial residue contacts at the SD -IBC interface in apo and IP 3 -bound NT. Several residues from the SD and IBC involved in interface communication are shown in Fig. S4. A salt bridge between Lys225 and Asp228 positions the SD relative to the IBC and stabilizes the β-interface by helping in arranging the residues at the interface. Notably, a network of hydrophobic and electrostatic interactions via Val33, Leu32, Arg54, and Lys127 from SD and Asp444, Phe445, Asp448, Ala449, Val452, and Leu476 from IBC-α are involved in maintaining the α-interface.
Besides identifying residue contacts seen in static crystallographic structures, MD simulations of the NT provided further insight into the allosteric rearrangements in interfacial residue network during IP 3 binding. The detailed analysis of the apo system revealed the formation of a bidentate hydrogen bond (H-bond), involving residues from all the three NT domains -SD, IBC-β and IBC-α. In the apo NT, Arg54 of SD forms a bidentate H-bond with Glu283 and Asp444 of IBC-β and IBC-α, respectively (Fig. 3a). It can be noted from the figure that Arg54 is situated at the center of the triangulated structure formed by the three NT domains. The guanidinium side chain of Arg54 H-bonds with the side chain carboxylate group of Glu283 and Asp444, with an average distance of 2.8 Å and 3.2 Å respectively. The −NH group of guanidinium interacts with the −COO group of Glu283 while the guanidinium −NH 2 group H-bonds with the −COO group of Asp444. Interestingly, IP 3 binding at the IBC breaks this triad of interaction at the SD -IBC interface (Fig. 3b). The binding of IP 3 at the IBC allosterically disrupts the interaction of Arg54 with Asp444 and the former then rearranges its side chain in such a way that the guanidinium −NH 2 group H-bonds with the −COO of Glu283. To further evaluate the stability of this bidentate interaction, we calculated the probability distribution of H-bond distance in Arg54 -Glu283 and Arg54 -Asp444 interaction during the entire span of simulation. In apo NT (Fig. 3c), Arg54 -Asp444 distance distribution was nearly Gaussian, peaking at the average distance of 2.8 Å, whereas the Arg54 -Glu283 interaction distance was distributed with two peaks, first at 3.2 Å and the second one at 5.1 Å. IP 3 binding shifts this equilibrium interaction of Arg54 towards Glu283, with the measured Arg54 -Glu283 average distance of 2.8 Å (Fig. 3d). On the other hand, Arg54 was found to be H-bonding with Asp444 (peak at a distance of 2.8 Å) during only about 20% of the simulation time while the Arg54-Glu283 interaction persist in about 80% of the simulation trajectories. Moreover, the distribution clearly indicates that, in the apo NT, though Arg54 shows comparatively equal affinity towards Glu283 and Asp444, the interaction with Asp444 is found to be stronger than that with Glu283. The energy of Arg54 -Glu283 and Arg54 -Asp444 interaction was calculated as −5.8 ± 2.8 kcal/mol and −15.3 ± 0.6 kcal/mol respectively. Notably, Arg54 -Asp444 interaction energy was decreased to −4.7 ± 0.9 kcal/ mol in IP 3 -bound NT (system 2) while the Arg54 -Glu283 interaction energy was −6.6 ± 1.3 kcal/mol, further supporting the IP 3 -driven residue rearrangements at the SD -IBC interface. The IP 3 mediated shift in Arg54 -Glu283/Asp444 bidentate interaction can also support the IP 3 -evoked twist motion of SD. In the apo state, Arg54 of SD H-bonds with Glu283 and Asp444 and such a bidentate interaction at SD -IBC interface can lock and restrict the conformational flexibility of SD. The IP 3 binding breaks this triad of interaction by disrupting the Arg54 -Asp444 interaction at the α-interface and stabilises the Arg54 -Glu283 H-bond at the β-interface. Such www.nature.com/scientificreports www.nature.com/scientificreports/ a residue rearrangement weakens the α-interface allowing the SD to slide over the inner surface of IBC-α and at the same time, anchor the SD immensely on the inner surface of β-interface. The absence of Arg54 -Glu283 H-bonding in the IP 3 -bound crystal structure (PDBID: 3UJ0) 5 suggests that the crystal architecture may be representing one of the IP 3 -bound NT conformation explored during the simulation where Arg54 failed to make H-bond with Glu283. Biochemical studies in the past revealed that the mutation in Arg54 13 or Asp444 22 affected IP 3 -binding affinity of the receptor. However, the functional consequences of all these three residues are yet to be known to the best of our knowledge and future studies may look into this aspect.
To summaries, the present all-atom MD simulation of IP 3 R NT provides a dynamic view of conformational landscape of NT during and immediately following IP 3 binding and sheds light on the mechanism of IP 3 -mediated IP 3 R activation leading to the pore opening. In a biologically relevant, functional tetrameric form of IP 3 R, the monomeric units are arranged in such a way that the SD of one monomer interacts with the IBC-β of the neighboring monomer through the HS loop of SD (Fig. S5) 6 . It is proposed that the HS loop -IBC-β interaction holds the tetrameric IP 3 R in a closed state. The IP 3 binding evokes the IBC clam closure and causes the SD to twist towards IBC, disrupting the HS loop -IBC-β interaction that allosterically opens the ion conducting pore, which is located ~100 Å away from the IP 3 binding site 6 . In the present study, analyses of essential dynamics of the apo and IP 3 -bound NT domains clearly indicate that the IP 3 binding at IBC allosterically affect the SD dynamics and the latter experiences a characteristic twist motion. The twist motion causes the SD to swing towards IBC-α, which eventually leads to the translational displacement of the conserved HS loop of SD. The residue Tyr167 of the HS loop was found to play a prime role in coupling the IP 3 binding to gating and mutation of Tyr167Ala had shown to attenuate the IP 3 -evoked Ca 2+ release 10 . It is noteworthy that, during simulation, Tyr167 experienced a maximum displacement of ~13 Å in the IP 3 -bound system while the displacement was about ~8 Å in the apo system (Fig. S6a,b) and such dynamics of the HS loop can break the SD -IBC interaction at the inter-monomeric interface (Fig. S5). www.nature.com/scientificreports www.nature.com/scientificreports/ IP 3 stabilizes IBC domain in the absence of the SD. So far, with the aid of extensive MD simulations, we have explored the IP 3 driven conformational dynamics in IP 3 R NT. To further investigate the effect of hitherto identified SD twist motion on the IBC dynamics, we simulated two additional systems of apo and IP 3 -bound IBC in the absence of the SD (system 3 and 4, Table 1). As the RMSD in Fig. S7 shows, the apo IBC (system 3), in the absence of the SD, underwent significant conformational fluctuation from its starting configuration during the course of simulation (average RMSD of 8.4 ± 1.9 Å). However, in the presence of the SD (system 1), IBC gets stabilized after an initial equilibration phase of 40 ns with minimal structural deviation during rest of the simulation (average RMSD of 4.8 ± 1.2 Å). This suggests that, in the apo state, IBC experiences conformational constraints in presence of the SD (i.e., as a part of NT, system 1) whereas the IBC becomes increasingly dynamic in the absence of the SD. Interestingly, IP 3 binding reverses this dynamic behavior of IBC. As shown in Fig. S7, the IP 3 -bound IBC of NT (i.e., in presence of SD, system 2) fluctuates significantly along the 300 ns of simulation, whereas the binding of IP 3 stabilized the IBC in the absence of the SD (system 4 in Table 1). The RMSD values from the figure clearly indicate that, the liganded IBC is highly dynamic in presence of the SD (average RMSD of 7.1 ± 1.6 Å). On the other hand, in the absence of the SD, RMSD of IP 3 -bound IBC reached a stable plateau after an initial 40 ns of equilibration phase suggesting that the domain is stabilized during rest of the simulation (average RMSD of 5.2 ± 1.4). To further examine the local structural transformations of IBC in the presence and absence of the SD, we compared the residue-wise C α RMSF of systems 1-4. The calculated RMSF values also corroborate the observation from RMSD analyses. As illustrated in Fig. 4a, in the absence of the SD, the apo and IP 3 -bound IBC exhibit distinct dynamic properties. The IBC residues were comparatively more flexible in the apo state (system 3) than the IP 3 -bound state (system 4). However, as seen from RMSD analyses, residue-wise RMSF also shows that the presence of the SD has influenced the dynamic behavior of IBC residues (Fig. 4a, black and red lines). In the presence of the SD (ie., NT), IBC residues experienced increased thermal fluctuations in the IP 3 -bound state as compared to the apo state.
Furthermore, to understand the domain movements associated with differential dynamics of IBC in the presence and absence of the SD, we resorted to PCA to identify the essential motions in the protein. As shown in Fig. S2, the top two eigenvectors (PCs) have captured significant motions of apo and IP 3 -bound IBC from each simulation ensemble. The projection of MD ensemble onto the 2D plane defined by the top two PCs plotted in Fig. S3b provides a comparison of conformations sampled by IBC of system 3 and 4. As the figure shows, although conformational overlap is observed, apo IBC spans a larger space in contrast to that of the IP 3 -bound IBC. This suggests that the dynamics of IBC is restricted in the liganded state as compared to the apo state. This is consistent with the results from the RMSF analysis, indicates that residues with high RMS fluctuations contribute significantly to the dominant motions of the IBC.
The porcupine plots showing the directional motion of apo and IP 3 -bound IBC along the direction of PC1 and PC2 in presence and absence of the SD are presented in Fig. 4b-e. Noticeably, in the apo state (Fig. 4b,c), all the secondary structure elements of IBC were found to be profoundly dynamic in the absence of the SD (system 3), whereas the IBC exhibits comparatively more restricted motions in presence of the SD (system 1, Fig. 2a,b). In the absence of the SD, apo IBC underwent a bending motion where the βand αdomains of IBC moved towards each other (Movie S3). However, presence of the SD (i.e., in apo NT, system 1) has hampered the bending motion of IBC by severely attenuating the dynamics of IBC-β (Fig. 2a,b, Movie S1). Interestingly, IP 3 binding reverses this characteristic motion observed in the IBC in presence and absence of the SD. In the absence of the SD, binding of IP 3 has restricted the overall dynamics of IBC (system 4, Fig. 4d,e, Movie S4). The IP 3 in the binding site prevented the bending motion of IBC domains, which was observed in the apo IBC (system 3, Fig. 4b,c). This is expected as the binding site of IP 3 located at the IBC-β -IBC-α interface and the bound IP 3 imparts local stability by preventing the collapse of the interface that was observed in apo IBC. A similar observation was also made from an earlier simulation study with IBC, where the dissociation of IP 3 conferred larger domain fluctuations to the protein 9 . However, in contrast to the reduced dynamics of IBC in system 4, IP 3 -bound IBC in presence of the SD (i.e., IBC of IP 3 -bound NT in system 2, Fig. 2c,d) shows increased domain movements, potentially due to the clam closure dynamics and twist motion of SD in the IP 3 -bound NT. The binding of IP 3 at the IBC allosterically weakens the SD -IBC interface that frees the IBC whereas the stronger SD -IBC interface in the apo NT limits the IBC motions. Together, differential dynamics of IBC observed here suggests that a strong inter-domain communication prevails between SD and IBC where both of them influence each other's dynamic behavior and the binding of ligand have long-range influence on domain dynamics.

SD suppresses the IP 3 affinity by shuffling the mode of interaction between the IBC and IP 3 .
Though the SD is critical for coupling the IP 3 -induced events to channel-opening, it is well known to suppress the affinity of the IBC or the full IP 3 R protein for IP 3 binding 11,12 . The IP 3 binding affinity of the NT or the whole receptor is reduced by more than one order of magnitude compared to that of IBC alone. Experimentally measured dissociation constants (K d ) of IP 3 from different studies showed the higher affinity of IBC for IP 3 in the absence of the SD 11,12,14,23 . The binding free energy (∆G) of IP 3 at 296 K is calculated to be 37.1 ± 0.1 kJ/mol for the NT whereas the IBC, in the absence of the SD, shows higher binding affinity with ∆G of −43.2 ± 0.1 kJ/mol 14 . Furthermore, study with different IP 3 R isoforms revealed that IP 3 binding affinity of the IBC is almost identical in the absence of the SD suggesting the critical role of SD in allosterically modulating the IP 3 binding to the receptor 24 . Though previous studies revealed the inhibitory action of the SD, the precise molecular mechanism underlying the reduced binding affinity of IP 3 in presence of the SD is not yet understood clearly. As discussed in the following sections, the present study with the aid of sub-microsecond MD simulations unveiled the molecular details of 'suppressing action' of the SD.
Firstly, to get a qualitative comparison of binding affinity and to complement the experimental observation, we computed the free energy of binding of IP 3 to NT (i.e., presence of the SD, system 2) and IBC alone (i.e., www.nature.com/scientificreports www.nature.com/scientificreports/ absence of the SD, system 4) using MM/PBSA approach (see methodology section). The calculated binding ∆G was found to be −70.4 ± 1.2 kcal/mol and −79.5 ± 1.1 kcal/mol for NT and IBC respectively. In agreement with the experimental findings, it is evident from the binding free energy that the IBC, in the absence of the SD, shows comparatively higher affinity for IP 3 than the whole NT. In order to unravel the molecular mechanism underlying such differential binding affinity, we further analyzed the interaction of individual protein residues with IP 3 in terms of pairwise interaction energy. The pairwise interaction energy of residues lining the IP 3 -binding pocket is presented in Fig. 5a. Surprisingly, most of the pocket residues interact with IP 3 with relatively similar strength except Arg269. The interaction energy of Arg269 with IP 3 was about 15 units more in the absence of the SD as compared to that in the presence of the SD. In presence of SD (system 2), Arg269 -IP 3 interaction energy was calculated to be −17.3 ± 0.6 kcal/mol. Markedly, in the absence of the SD (system 4), Arg269 of IBC interacts strongly with IP 3 with energy of −33.8 ± 0.8 kcal/mol. In brief, the energy calculations suggest that the Arg269 -IP 3 interaction plays a crucial role in determining IP 3 R's affinity for IP 3 in presence and absence of the SD.
To obtain a better understanding of the receptor -IP 3 interaction, we compared the ensemble-averaged structures of IP 3 -bound NT (system 2) and IBC (system 4) complexes. This is shown in Fig. S8a,b. For clarity, only the ligand and the adjacent protein residues that are involved in direct and consistent interactions are shown. It can be noted from the figure that the sugar moiety and the phosphate groups of the IP 3 efficiently interact with the www.nature.com/scientificreports www.nature.com/scientificreports/ residues spanning the IBC -β and -α domains 25 . As expected, the binding pocket is lined by a number of basic and polar residues to accommodate a highly negatively charged molecule like IP 3 . Three phosphate groups -P1, P4, and P5 -of the IP 3 molecule are involved in an extensive network of H-bonds with five arginine residues -Arg265, Arg269, Arg504, Arg511 and Arg568. The Lys508 that lies at the interior of binding pocket also forms H-bond with the P5 of IP 3 . During the entire span of simulation, Thr267 and backbone amino group of Gly268 was found to interact with P4 phosphate group. The C6′-OH of the IP 3 forms H-bond with either Arg504 or Tyr567 in both system 2 and system 4, whereas other hydroxyl groups fail to make any direct interaction with protein residues. Such a less significant role of IP 3 hydroxyl groups in its binding with IP 3 R is reported from the biophysical studies with IP 3 derivatives 26 .
As discussed above, most of the pocket residues show similar pattern of interaction with IP 3 in the presence and absence of the SD. However, very interestingly, the SD seemed to influence the mode of interaction of Agr269 with IP 3 . In the presence of the SD (system 2), the guanidinium group of Arg269 was found to interact effectively with the P1 (phosphate group attached to the C1′ of the inositol sugar ring) of the bound IP 3 molecule (Fig. S8a). On the other hand, in the absence of the SD (system 4), Arg269 of the IBC makes multiple interactions with the bound IP 3 , where the guanidinium group of the residue H-bonds with P5 (phosphate group attached to the C5′ of the inositol sugar ring) and C6′-OH of the IP 3 molecule (Fig. S8b). In both systems, the guanidinium −NH/− NH 2 groups were involved in H-bond with the PO 3 2− /C6′-OH groups of the IP 3 molecule. As evident from the figure, Arg269, which interacts with the P1 of the IP 3 in the presence of the SD (ie., in NT), is repositioned in the absence of the SD and forms H-bond with the P5 and C6′-OH groups of IP 3 that lies deep inside the pocket. We further examined the stability of Arg269 -IP 3 interaction by measuring the Arg269 -PO 3 2− distance along the 300 ns simulation trajectories (Fig. 5b,c). In the presence of the SD, the guanidinium group of Arg269 interacts with the P1 having a H-bonding distance of ~3.8 Å throughout the simulation whereas the P5 was found to be located ~8 Å away from the Arg269 (Fig. 5b). Strikingly, Arg269 flipped its mode of interaction in the absence of the SD (system 4). As shown in Fig. 5c, rather than interacting with the P1 (distance of ~7 Å) as seen in the presence of the SD, the guanidinium group of the residue H-bonds with P5 in the absence of the SD, maintains a distance of ~3.8 Å. A similar analysis of time evolution of Arg269 -C6′-OH distance during simulation is depicted in Fig. S9. The guanidinium −NH 2 -C6′−OH average distance was measured to be 5.8 ± 1.0 Å in the IP 3 -bound NT (system 2) as compared to the H-bond distance of 3.9 ± 0.5 Å in the IP 3 -bound IBC (system 4).
A closer look at the initial configurations from the simulation trajectories has also provided further insight into the SD-invoked Arg269 flipping. In the crystal structures, Arg269 was found to interact with the P5 of the IP 3 in both IP 3 -bound IBC (PDB ID: 1N4K) as well as NT (PDB ID: 3UJ0) with a distance of about 2.7 Å. It is to be noted that the region spanning residues 268-270 that includes Arg269 is not resolved in any of the apo-NT crystal structures (PDB IDs: 3UJ4, 3T8S apo, 5 × 9Z, 5XA0) as well as an IP 3 -bound NT structure (PDB ID: 3T8S). It is, however, apparently resolved in the IP 3 -bound NT structures of the Cys-less rat IP 3 R1 (PDB ID: 3UJ0) 5 as well as www.nature.com/scientificreports www.nature.com/scientificreports/ the mouse IP 3 R1 (PDB ID: 5GUG, 5XA1) 7 . But, for all these structures, the crystal structure of mouse IP 3 R1 IBC (PDB: 1N4K) that is clearly devoid of the SD, effectively served as the reference for the initial molecular replacement 5,7 . We suspect that this region containing Arg269 is very flexible and this is supported by our simulation as well as close inspection of the available crystal structures of the NT and full length IP 3 R1 4,5,7 . Intriguingly, within the initial 3 ns of the simulation in our study, Arg269 within the IP 3 bound NT breaks its interaction with P5 and flips towards P1, measuring a reduction in Arg269 -P1 distance from 8 Å to ~3.8 Å (Fig. 5b inset) and retained the distance of ~3.8 Å during rest of the simulation. It can be postulated that, in the starting NT configuration, the SD influenced the overall IBC architecture where the loop bearing the Arg269 was unable to occupy in its initial conformation, forcing the Arg269 to reposition itself to interact with the P1 of the IP 3 .
To strengthen our finding further that the SD modulates Arg269 flipping and interaction with IP 3 , we created a new simulation system (system 5) by removing the SD from the final conformation of IP 3 -bound NT obtained after the 300 ns simulation of system 2. Hence, the newly generated system consists of IP 3 -bound IBC of the previously simulated NT where the Arg269 interact with the P1 of the IP 3 . Here, our aim was to examine whether the IBC can sense the absence of the SD and reorient the Arg269 in such a way that it interacts with P5 as seen in IP 3 -bound IBC (i.e., in the absence of the SD, system 4). As in the above section, we analysed the Arg269 -IP 3 interaction in system 5 by monitoring the guanidinium -PO 3 2− distance along the simulation time. In Fig. 6, we combined the calculated distance from system 2 and 5. To our surprise, IBC in system 5 senses the absence of the SD during the course of simulation that eventually leads to the necessary structural rearrangements in the IBC with flipping of the Arg269 towards P5. As shown in Fig. 6, after ~280 ns of simulating the SD-removed NT (i.e., IBC), Arg269 -P5 distance decreased from about 10 Å to the H-bonding distance of ~3.8 Å and, at the same time, Arg269 -P1 distance increased to about 8 Å. This suggests that Arg269 breaks its interaction with P1 and flips towards P5 after ~280 ns of simulating the SD-removed IP 3 -bound IBC. Remarkably, the flipping of Arg269 was also reflected in the IP 3 binding affinity. As expected, the IP 3 binding affinity of IBC for IP 3 was increased in the absence of the SD, with a rise in ∆G of binding from −71.9 ± 2.1 kcal/mol, in the IP 3 -bound NT, to −80.9 ± 0.9 kcal/mol, after the removal of the SD. The flipping has also altered the Arg269 -IP 3 interaction energy (Fig. 5a). By shifting the site of interaction from P1, in presence of the SD, to P5, after the removal of the SD, Arg269 has strengthened its interaction with IP 3 with an increase in energy from −17.3 ± 0.6 kcal/mol to −32.6 ± 0.1 kcal/mol. It is noteworthy that the calculated energies are in good agreement with that calculated from IP 3 -bound IBC in system 4 (Fig. 5a, blue bar). A movie showing the flipping mechanism of Arg269, in the absence of the SD, in system 5 is presented in Movie S5. Together, these analyses suggest that the SD allosterically modulates the IBC interaction with IP 3 and Arg269 plays a crucial role in determining the binding affinity of IBC for IP 3 in presence and absence of the SD. Higher IP 3 binding affinity in the absence of the SD can be ascribed to the multiple interactions of Arg269 with IP 3 , involving Arg269 -P5 and Arg269 -C6′-OH, in contrast to the Arg269 -P1 interaction in the presence of the SD (Figs 5b,c, S8 and S9).

Conclusions
As with most other ligand-gated ion channels, a detailed understanding of the activation gating of IP 3 Rs remains a holy grail in the field. Although, the currently available structural and biochemical information has significantly advanced our understanding of the structure-to-function relationship in IP 3 R, details of the dynamics of IP 3 -mediated domain movements within IP 3 R NT remained largely unexplored. In the present study, extensive MD simulations of apo and IP 3 -bound NT and analysis of essential dynamics have provided a detailed picture of correlated domain movements during and immediately following IP 3 binding. Compared to the apo form, SD of IP 3 -bound NT showed characteristic twist motion, where the SD move towards the IBC-α. During the binding process, IP 3 has brought together the β and α domains of IBC, leading to the 'clam closure' . Our results demonstrated a dynamic view of unique structural plasticity acquired by the IP 3 R NT as a consequence of IP 3 Figure 6. Removal of the SD from the IP 3 -bound NT flips the Arg269 interaction with the IP 3 . Time evolution of Arg269 guanidinium -PO 3 center of mass distance combined from system 2 and system 5 (see Table 1 for system definition). The SD was removed from IP 3 -bound NT (system 2) after the initial 300 ns of simulation and thus generated IP 3 -bound IBC (system 5), which was simulated for another 400 ns. Corresponding free energy of IP 3 binding before and after the removal of SD are shown. Color scheme: Arg269 -P1, green and Arg269 -P5, red.
www.nature.com/scientificreports www.nature.com/scientificreports/ binding and support the observations inferred from crystallographic studies 5 . The H-bond and energetic analyses revealed the IP 3 -induced molecular rearrangements at the SD -IBC interface, where Glu283-Arg54-Asp444 triad of interaction breaks during the IP 3 binding. The Arg54-Asp444 interaction breaks while the Arg54 maintains its interaction with Glu283. MD simulations further show the existence of inter-domain dynamic correlation in the IP 3 R NT and reveal that the SD is critical to the conformational dynamics of IBC. Very interestingly, the study unraveled the SD-dependent flipping of interaction of Arg269 between P1 and P5 phosphate groups of IP 3 , which seems to play a major role in governing the IP 3 binding affinity of IBC in the presence/absence of the SD. It is to be noted that the residues identified in this study (Arg54, Arg269, Glu283 and Asp444) are conserved across different IP 3 R isoforms and in different species (Fig. S10). To conclude, the present study has brought out the importance of alternate conformations of IP 3 R NT needed for its function and the comprehensive understanding of the IP 3 -binding mechanism of IP 3 R derived here could help to understand the physiological activation and pharmacological regulation of IP 3 Rs.

Materials and Methods
A series of unbiased all-atom MD simulations of IP 3 R NT (residues 1-604) in its apo and IP 3 -bound states ( Table 1) were performed. The crystal structure of the IP 3 R NT was used to model the full-length NT with residues 1-604. Thus, the atomic coordinates of the apo (PDB ID: 3UJ4) and IP 3 -bound (PDB ID: 3UJ0) IP 3 R NT obtained from the protein data bank were used in this study. The missing residues; 1-6, 76-84, 204, 269-270, 294-298, 373-381, 407, 426-427, 486-500, 525-550 and 578-604 were incorporated with the help of MODELLER9v16 27 using crystal structures of IBC (PDB ID:1N4K) and SD (PDB ID: 1XZZ) as templates. Additionally, the crystal configurations were also missing a sequence of residues 319-351. Since we noted that a very low sequence similarity with available structures in databases, the 3D model of this sequence was constructed using the ab initio protein structure predictor, QUARK 28 . Quark builds the models based on structures of small fragments of 1-20 residues, which was later assembled by replica exchange Monte Carlo simulations using an atomic-level knowledge-based force field. The modeled structure of the amino acid sequence 319-351 was merged with the crystal conformation of the apo (system 1) and liganded NT (system 2) and optimized by the conjugate gradient energy minimization method in Modeler program 27 . Further, the apo (system 3) and IP 3 -bound (system 4) IBC-only structures of IP 3 R NT were generated from the full-length NT by deleting the N-terminal SD of residues 1-224.
In each of the above protein structures, hydrogen atoms were added using H++ server 29 , maintaining the ionization state of the residues at the pH of 7.4. A set of partial atomic charges for IP 3 was obtained via quantum calculations. A B3LYP geometry optimization procedure was performed using Gaussian 09 with the 6-311 + G* basis set 30 . The atom-centered RESP charges were calculated via fits to the electrostatic potentials obtained from the calculated wave functions 31 . The interaction parameters for the IP 3 were adopted from a previous work, where the forcefield parameters for phosphorylated inositols were developed based on OPLS-AA/AMBER framework 32 . Subsequently, the structures were energy minimized for 1000 steps using the steepest descent algorithm and another 1000 steps by the conjugate gradient method, using AMBER 14.0 simulation package 33 . After relaxing the added atoms in gas phase, each structure was solvated in cubic periodic box of explicit water with water molecules extending 14 Å outside the protein complex on all sides. The 3-site TIP3P water model was chosen to define the water molecules. In system 1 and 2, the simulation box contained ~45,000 water molecules, while ~22000 water molecules were present in the systems 3 and 4. Charge neutrality was maintained by adding one Cl − ion in system 1 and five, three and four Na + ions in system 2-4 respectively. The solvated systems were subjected to extensive energy minimization and thermalization with raising the temperature gradually to 300 K in canonical ensembles, maintaining harmonic restraints on crystallized heavy atoms of protein and ligand. Afterwards, solvent density was adjusted under isobaric and isothermal conditions at 1 atm and 300 K and the harmonic restraints were slowly reduced to zero. To enable the volume variation, simulations were performed in an NPT ensemble using the Berendsen thermostat and barostat. The systems were equilibrated for 10 ns, with a 2 fs simulation time step. The long-range Coulombic interactions were treated using Particle Mesh Ewald sum with a cut off of 12 Å applied to Lennard-Jones interactions and all bonds involving hydrogen atoms were constrained using the SHAKE algorithm. Each system was later subjected to 300 ns of production run, saving the trajectories at an interval of 2 ps for further analysis. All simulations were performed using the AMBER 14.0 molecular dynamics simulation package with the AMBER ff14SB force field 33 . Additionally, a new simulation system (system 5, Table 1) was generated from the final conformation of IP 3 -bound NT obtained from 300 ns of simulation of system 2. This 'SD-knockout' system was prepared by removing the N-terminal SD of 224 residues from the IP 3 -bound NT. Like other systems, simulation of system 5 was performed for 400 ns following a similar protocol as stated above. Together, a total of 1.6 µs of atomistic simulations were carried out in this work.
Principal component analysis (PCA) is employed to explore the dominant motions of apo and IP 3 -bound IP 3 R NT region. PCA tool has been used in the past to locate the essential motions in complex systems with high dimensionality, like in proteins, and the detailed methodology can be found elsewhere 34 . PCA was performed in this study using pyPcazip program package 35 . The eigenvectors from PCA were visualized as porcupine plots using VMD 36 , where the direction of the porcupine corresponds to the direction of the motion and the length of the porcupine denotes the magnitude of the motion. The twist of the SD was calculated by measuring the angle defined by the two structures, the crystal structure and a snapshot from the simulation. A plane was defined on each structure by three points that is located on the center of mass of each domain (SD, IBC-β and IBC-α). The 'twist angle' defined as the dihedral angle between the two planes was calculated. PyMOL (The PyMOL Molecular Graphics System, Version 1.7.4, Schrodinger, LLC) and VMD 36 was used to generate all structural figures.