Differential dynamics of the serotonin1A receptor in membrane bilayers of varying cholesterol content revealed by all atom molecular dynamics simulation.

Abstract The serotonin1A receptor belongs to the superfamily of G protein-coupled receptors (GPCRs) and is a potential drug target in neuropsychiatric disorders. The receptor has been shown to require membrane cholesterol for its organization, dynamics and function. Although recent work suggests a close interaction of cholesterol with the receptor, the structural integrity of the serotonin1A receptor in the presence of cholesterol has not been explored. In this work, we have carried out all atom molecular dynamics simulations, totaling to 3 μs, to analyze the effect of cholesterol on the structure and dynamics of the serotonin1A receptor. Our results show that the presence of physiologically relevant concentration of membrane cholesterol alters conformational dynamics of the serotonin1A receptor and, on an average lowers conformational fluctuations. Our results show that, in general, transmembrane helix VII is most affected by the absence of membrane cholesterol. These results are in overall agreement with experimental data showing enhancement of GPCR stability in the presence of membrane cholesterol. Our results constitute a molecular level understanding of GPCR-cholesterol interaction, and represent an important step in our overall understanding of GPCR function in health and disease.


Introduction
The G protein-coupled receptor (GPCR) superfamily comprises the largest and most diverse group of proteins in mammals, and is involved in information transfer (signal transduction) from outside the cell to the cellular interior (Chattopadhyay, 2014;Pierce et al., 2002, Rosenbaum et al., 2009. GPCRs regulate cellular responses to a variety of stimuli, and mediate multiple physiological processes. Due to this reason, GPCRs have emerged as major drug targets in all clinical areas (Heilker et al., 2009). It is estimated that $50% of clinically prescribed drugs target GPCRs (Schlyer & Horuk, 2006).
GPCRs are integral membrane proteins with seven transmembrane passes. A large portion of these receptors therefore remains in contact with the membrane. This raises the obvious possibility that the membrane could be an important regulator of receptor structure and function. Exploring the interaction of GPCRs with membrane lipids assumes significance in light of the fact that cells have the ability of varying the lipid composition of their membranes in response to a variety of stress and stimuli, thereby changing the environment and the activity of the receptors in their membranes. The organization and function of GPCRs have been shown to be dependent on the membrane lipid composition, specifically cholesterol (Burger et al., 2000;Jafurulla & Chattopadhyay, 2013;Oates & Watts, 2011;Pucadyil & Chattopadhyay, 2006). Cholesterol is one of the important components of eukaryotic cell membranes and plays a crucial role in membrane dynamics and organization (Chaudhuri & Chattopadhyay, 2011;Mouritsen & Zuckermann, 2004;Simons & Ikonen, 2000). The role of membrane cholesterol in the organization and function of GPCRs represents an exciting area of research (Burger et al., 2000;Jafurulla & Chattopadhyay, 2013;Oates & Watts, 2011;Pucadyil & Chattopadhyay, 2006). Nonetheless, the detailed mechanism underlying the effect of membrane cholesterol on the structure and function of GPCRs is not clear and appears to be complex (Paila & Chattopadhyay, 2009. A possible mechanism by which membrane cholesterol has been proposed to modulate GPCR organization and function is through a direct (specific) interaction, which could induce a conformational change in the receptor. A second possibility proposes an indirect way by altering the membrane physical properties. Yet another possibility could be a combination of both. An important common feature observed in recently solved high resolution crystal structures of GPCRs [such as b 1 -adrenergic receptor (Warne et al., 2011), b 2 -adrenergic receptor Hanson et al., 2008), A 2A adenosine receptor (Liu et al., 2012), and metabotropic glutamate type 1 receptor (Wu et al., 2014)] is the close association of cholesterol molecules with the receptor. Cholesterol binding motifs such as the cholesterol recognition/interaction amino acid consensus (CRAC) motif have been identified in GPCRs (Jafurulla et al., 2011), but the exact role of these motifs in preferential cholesterol binding to GPCRs is still not clear.
Molecular dynamics simulations have emerged as important tools in understanding GPCR-membrane interactions (Grossfield, 2011;. Previous simulations with cholesterol-rich membrane bilayers have suggested that cholesterol affects local structural properties of rhodopsin (Khelashvili et al., 2009) and human A 2A adenosine receptor (Lyman et al., 2009). Recent microsecond time scale simulations of A 2A adenosine receptor and the serotonin 1A receptor have pointed toward a close, but highly dynamic receptor-cholesterol association (Lee & Lyman, 2012;Sengupta & Chattopadhyay, 2012). In the b 2adrenergic receptor, cholesterol has been recently shown to be implicated in decreased tryptophan side-chain dynamics leading to a conformational ''lock'' (Cang et al., 2013) and in oligomerization of the receptor (Prasanna et al., 2014). The presence of cholesterol in the membrane therefore appears to modulate the structural plasticity of GPCRs, although the exact role of cholesterol remains to be explored.
The serotonin 1A receptor is an important member of the GPCR superfamily and is implicated in the generation and modulation of various cognitive, behavioral and developmental functions (Kalipatnapu & Chattopadhyay, 2007;Müller et al., 2007;Pucadyil et al., 2005). The agonists (Blier & Ward, 2003) and antagonists (Griebel, 1999) of this receptor represent major classes of molecules with potential therapeutic applications in anxiety-or stress-related disorders. As a result, the serotonin 1A receptor serves as an important drug target for neuropsychiatric disorders such as anxiety and depression (Celada et al., 2013). In our previous work, we comprehensively demonstrated using a variety of experimental approaches that membrane cholesterol plays a crucial role in the organization, dynamics and function of the serotonin 1A receptor (reviewed in Jafurulla & Chattopadhyay, 2013;Pucadyil & Chattopadhyay, 2006). Although cholesterol binding ''hot-spots'' on the receptor surface have been identified (Jafurulla et al., 2011;Sengupta & Chattopadhyay, 2012), the effect of cholesterol binding on receptor structure and dynamics remain still unexplored.
In this work, we have carried out all atom molecular dynamics simulations to analyze the effect of cholesterol on the structure and dynamics of the serotonin 1A receptor. We have simulated the serotonin 1A receptor in cholesterol-rich 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane bilayers, chosen to mimic the natural biological environment in which the serotonin 1A receptor resides. Simulations in POPC bilayers alone were carried out as a control. Our results show that the serotonin 1A receptor exhibits differential dynamics in POPC/cholesterol bilayers relative to that in POPC bilayer. Specifically, we show that transmembrane helix VII has lower conformational flexibility in the presence of membrane cholesterol in the time scales of the simulations. These results are in overall agreement with experimental results which show that membrane cholesterol increases the stability of the serotonin 1A receptor (Saxena & Chattopadhyay, 2012) and the b 2 -adrenergic receptor (Yao & Kobilka, 2005), a closely related GPCR which has a considerable sequence similarity with the serotonin 1A receptor in the transmembrane region (Paila et al., 2011b).

System setup
All atom molecular dynamics simulations of the serotonin 1A receptor were carried out with the receptor embedded in POPC and POPC/cholesterol bilayers. The concentration of cholesterol was 30 mol% in POPC/cholesterol bilayers. Two different systems were considered. In the first set, simulations were performed with a previously developed homology model of the serotonin 1A receptor (Paila et al., 2011b). The model was chosen since it has been shown to be stable and was able to reproduce cholesterol interactions in a coarse-grain model (Sengupta & Chattopadhyay, 2012). An appropriate number of Clcounter ions were added to neutralize the charge of the system. The membrane embedded receptor was simulated at two different ionic strengths (in the presence of 0.1 M NaCl and with only Clcounter ions). These systems are referred to as model 1 (0.1 M NaCl) and model 2 (0 M NaCl) throughout the rest of the paper. In the second set, a new homology model was built based on the recent crystal structure of the serotonin 1B receptor . This is referred to as model 3 and is described below. An appropriate number of Clcounter ions were added to neutralize the charge of the system. All the systems were generated using CHARMM-GUI (Jo et al., 2008). A comprehensive list of all simulations performed is shown in Supplementary Table 1 (all Supplementary material available online). The residues corresponding to different regions of the receptor consisting of seven transmembrane helices (I-VII), the carboxy terminal helix VIII and seven loop regions are shown in Supplementary Table 2.

Homology model
The structure of the serotonin 1A receptor was modeled based on the serotonin 1B crystal structure (PDB Code: 4IAR ). The amino acid sequence of the serotonin 1A receptor exhibited increased consensus with that of serotonin 1B receptor (query cover: 98%, identity: 40%, E-value: 3 e-103). The cytochrome b562 insert between transmembrane helices V and VI was removed prior to building the model. To model loop 6 (intracellular loop 3), an independent BLAST search was performed. Chain A from the crystal structure of the N-terminal core of Bacillus subtilis inorganic pyrophosphatase (PDB: 1WPN, query cover: 44%, identity: 26%, E-value: 1.3) was chosen to model loop 6. Subsequently, a replica exchange molecular dynamics simulation was performed of the homology model of loop 6 in water. A large number (208) of replicas were considered spanning the temperature range from 200-580 K. Position restraints were applied on the backbone atoms of the terminal amino acid residues. Exchange time among the replicas was set to 100 ps with the acceptance ratio greater than 10%. Each replica was simulated for $4 ns with a total simulation time of 800 ns. The conformation with highest secondary structure similarity to previous NMR studies (Chen et al., 2011) was chosen for further analysis. The homology model of the serotonin 1A receptor (without loop 6) generated from the crystal structure of the serotonin 1B receptor and the loop 6 conformation obtained from replica exchange molecular dynamics were used as templates to generate the complete homology model of the serotonin 1A receptor. EasyModeller ver 4.0 GUI package using MODELLER ver. 9.14 was used for model building. Quality check on the model assessed using ProQM server showed a global quality of 0.545 (0: worst, 1: best).

Simulation parameters
All atom molecular dynamics simulations were performed of the membrane embedded receptor. The first set of simulations (model 1 and 2) were performed using the program NAMD 2.7b1 (Kalé et al., 1999) with the PARAM27 version of the CHARMM force-field (MacKerell et al., 1998). The second set of simulations (model 3) were performed using GROMACS version 4.5.5 (Pronk et al., 2013). The CHARMM force-field version 36 was used for the study. The main difference between the two versions is the back-bone parameters to correct for the over stabilization of the helical stretches in the protein. The system was first minimized for 1000 steps using conjugate gradient method with the protein and lipid headgroup atoms fixed to allow reorganization of water molecules and lipid acyl chains. In the subsequent step, the temperature of the system was increased from 0-300 K for 0.5 ns, keeping the simulation box constant. During the simulation, the protein and lipid headgroup atoms were kept fixed to allow the melting of the lipid tails consistent with the liquid crystalline (fluid) phase. Subsequently, 1000 steps of minimization were carried out with all protein atoms constrained with harmonic constraints with a force constant of 1 kcal/molÅ 2 . This was followed by simulations for 0.5 ns with the same positional restraint on protein atoms. Finally, dynamics was performed for 0.5 ns with no positional restraint to equilibrate the system. Periodic boundary conditions were used in all three directions during the simulation. Production simulation was carried out using methodologies developed for the CHARMM 27 and CHARMM 36 parameter set, respectively. In model 3, the pressure was maintained using the Parinello-Rahman coupling scheme, such that the lateral and perpendicular pressures are coupled independently at 1 bar (compressibility 4.5 Â 10 À5 bar À1 ). The temperature of the system was maintained at 300 K through Langevin damping with a coefficient of 1 ps À1 . Particle-mesh Ewald (PME) (Batcho et al., 2001) was used to compute the electrostatic interaction with a real space cut-off of 10 Å . We used the SHAKE algorithm to constrain all bonds involving hydrogen and consequently an integration time step of 2 fs was used. Van der Waals energies were calculated using 10 Å cut-off and pair list distance was calculated at 12 Å .
Model 1 and 2 represent simulations for 250 ns and model 3 for 1 ms.

Analysis
The analysis of the results was carried out using VMD 1.8.7 (Humphrey et al., 1996) and GROMACS analysis tools. For the RMSD calculation, the protein structures were aligned to the first frame of the simulation in respect to the particular region for which the RMSD was being measured. RMSDs were calculated with regard to both all atoms, including sidechains (all atom) and backbone atoms only, as mentioned. RMSF was calculated from the root mean fluctuations of the Ca atoms along the trajectory.

Results and Discussion
The serotonin 1A receptor exhibits differential dynamics in the presence of cholesterol Molecular dynamics simulations were performed with membrane embedded serotonin 1A receptor in POPC and POPC/ cholesterol bilayers. The concentration of cholesterol was 30 mol% in POPC/cholesterol bilayers in all cases. Three systems were considered to avoid prior bias due to the homology model. The first two models were based on a previously developed homology model (Paila et al., 2011b), which has been shown to be stable and able to reproduce cholesterol interactions (Sengupta & Chattopadhyay, 2012). The equilibrated homology model of the receptor was simulated in the presence (model 1) and absence (model 2) of 0.1 M NaCl. In addition, a new homology model was developed (model 3) based on a recent crystal structure of the serotonin 1B receptor .
To understand the influence of the membrane lipid bilayer environment on the structure and stability of the serotonin 1A receptor, we calculated the root mean square deviation (RMSD) of the serotonin 1A receptor in membrane bilayers of POPC and POPC/cholesterol. Figure 1 shows that, on an average, the RMSD of the serotonin 1A receptor in POPC bilayers is higher than that in POPC/cholesterol bilayers in the absence of salt. This is true when RMSD was calculated for all atoms (panels a, c and e) as well as backbone atoms (panels b, d and f). At longer time scales (model 3), the deviations are consistently higher in the absence of cholesterol. The same trend was also observed when only the transmembrane helices were considered ( Figure 2). However, the deviations along the trajectory are lower since the transmembrane helices exhibit reduced dynamics. The large jumps observed in Figure 1, but not in Figure 2, indicate that these correspond to structural rearrangements in the loop regions. The difference in RMSD could be indicative of the fact that the serotonin 1A receptor is less flexible in the cholesterol-rich membrane bilayer in the time scales of the simulations.

Dynamics of individual transmembrane helices in the absence of cholesterol
To further characterize the decrease in receptor dynamics in presence of cholesterol, we analyzed the RMSD of the individual structural components of the receptor. These include the transmembrane helices I-VII, the extracellular and intracellular loops connecting the helices, and the carboxy terminal helix VIII (see Supplementary Table 2). Figure 3 shows the RMSD of the individual transmembrane helices I-VII and helix VIII of model 1 in presence of salt. The RMSD profiles calculated for the intra-and extracellular loops in POPC and POPC/cholesterol bilayers for the three models are shown in Supplementary Figures 4-6. The trend in stability for the loop regions in the presence of cholesterol displayed considerable variation and large flexibility, irrespective of salt concentration. Relatively high dynamics was observed for loop 6 (intracellular loop 3). Interestingly, the contribution of this region to the total RMSD was highest. The RMSF of the residues in this region were also relatively high (Supplementary Figure 7). Superposition of the structures at the start and end of the simulation in POPC and POPC/cholesterol bilayers of the three models is shown in Figure 6, which shows that there is no major unfolding of the transmembrane helices during the course of the simulation in POPC or POPC/ cholesterol bilayers. Transmembrane helix VII appeared to show higher structural perturbation at longer time scales, but it is not clear if it is dependent on the initial homology model. Taken together, this supports the application of coarse-grain approaches, specially the MARTINI forcefield, in which the transmembrane helices are constrained to their starting structures with a harmonic potential (Prasanna et al., 2014;Sengupta & Chattopadhyay, 2012). However, loop 6 (third intracellular loop) which consists of $120 residues (see Supplementary Table 2) is mainly unstructured (disordered) and shows considerable rearrangement. Interestingly, this is the region of GPCRs that poses considerable challenge toward crystallization efforts due to its inherent conformational flexibility. In recent crystallographic efforts, this has been addressed by stabilizing this region using monoclonal antibody, or nanobody, or replacing with lysozyme Day et al., 2007;Rasmussen et al., 2011;Rosenbaum et al., 2007).

Proximity of cholesterol to the receptor
To analyze the probable occupancy sites of cholesterol around the transmembrane helices of the receptor, we calculated the radial distribution function, g(r), of cholesterol around each of the helices. As a representative example, the values for model 3 are shown in Figure 7. The transmembrane helices were subdivided into two parts corresponding to the upper and lower leaflets of the membrane. In general, the first peak of g(r) is found at the closest distance for transmembrane helix I (black), indicating a close interaction of cholesterol with these helices. Additionally, high cholesterol peaks were also observed for transmembrane helices IV (blue) and V (yellow). The interaction of cholesterol in models 1 and 2 was similar (Supplementary Figure 8). Close interactions with cholesterol was observed in transmembrane helices I, IV, V and VII. However, due to the shorter time scales of the simulations, the sampling of the cholesterol was limited (see Supplementary Table 3), and it is possible that the sampling of the interaction modes was not adequate. In general, the interaction with cholesterol was stochastic in agreement with differential cholesterol association observed in coarse-grain and atomistic simulations . Interestingly, in our previous work using coarse-grain MARTINI approach (Sengupta & Chattopadhyay, 2012), we showed stochastic cholesterol interactions sites, pointing toward higher occupancy of cholesterol at transmembrane helices I, V and VII.

Lipid binding site at the groove of transmembrane helices I and VII
We analyzed the presence of POPC molecules at the groove of transmembrane helices I and VII, which we have earlier shown to correspond to a lipid binding site in the b 2adrenergic receptor by coarse-grain simulations (Prasanna et al., 2015). Interestingly, this site has also been identified as a putative lipid binding site from the crystal structure of A 2 adenosine receptor (Liu et al., 2012). We identified a POPC molecule that was bound at the site, both in the presence and absence of cholesterol (Figure 8a and b). However, the absence of cholesterol at transmembrane helix I, in the vicinity of this lipid molecule reduces the acyl chain order of the lipid molecule at that groove ( Figure 8c). This effect can be also visualized from the snapshots (Figure 8a and b) where the POPC molecule is much more ordered in the presence of cholesterol. The increased fluctuations in the POPC molecule appear to lead to increased interactions with transmembrane helix VII and its subsequent structural deviations. The effect of cholesterol on transmembrane helix VII is therefore suggested to arise from a combination of its direct and indirect effects.

Cholesterol-mediated effects on GPCRs are moderate but crucial
There is a lack of consensus on the mechanism of cholesterol interaction from results of previous studies on molecular dynamics simulations of rhodopsin. Cholesterol has been suggested to be mainly excluded from the receptor surface (Grossfield et al., 2006;Pitman et al., 2005), although recent data suggest the presence of high cholesterol density sites (Horn et al., 2014;Khelashvili et al., 2009). Among these reports, the only work focusing on protein structural integrity (Khelashvili et al., 2009) reports structural differences in the kink regions in transmembrane helices I, II and VII. However, Lyman and co-workers (Lee et al., 2013) analyzed the effect of cholesterol on the A 2A adenosine receptor and reported a force-field-dependent dynamics of the kink regions. Another study (Shan et al., 2012) correlated the functionally relevant conformational flexibility of the serotonin 2A receptor with the interaction of cholesterol, but it was not clear whether this conformational flexibility was cholesterol-dependent. It is therefore evident that despite several studies attempting to understand the role of cholesterol, the effect of cholesterol on the structural integrity of GPCRs in general, and the serotonin 1A receptor in particular, is still not clear. Viewed from this perspective, our work constitutes one of the first attempts to monitor receptor dynamics in the presence and  absence of cholesterol, and address the role of the individual structural elements. We have used two different models and two force-fields in an attempt to reduce the bias of the initial structure and force-field-related structural changes reported earlier. Our study is limited to sub-microsecond time regimes and longer time scale studies may reveal further finer differences. Although reduced dynamics is observed in POPC/cholesterol bilayers in the absence of salt, the reverse trend was observed in the presence of 0.1 M NaCl. While we propose that the differential dynamics observed in these models needs to be probed further, these results nonetheless validate the major finding that cholesterol affects receptor dynamics.
Our results assume additional significance in the context of the long time scale coarse-grain simulations that are being used to analyze receptor organization and protein-lipid interactions (Johner et al., 2014;Johnston et al., 2012;Khelashvili et al., 2009;Mondal et al., 2013;Periole et al., 2007Periole et al., , 2012. GPCR plasticity is important for receptor activity and signaling, and a priori, it is not clear whether a coarse-grain model would be suitable to study protein-lipid interactions. In particular, in our previous work regarding the homology model of the receptor (Paila et al., 2011b) and coarse-grain analysis of lipidprotein interactions (Sengupta & Chattopadhyay, 2012), the secondary structural elements were constrained and it was not possible to study the flexibility and dynamics of the receptor itself. The absence of large conformational changes observed in our current results, and the correspondence between cholesterol ''hot-spots'' between the rigid coarse-grain model and the flexible atomistic models, corroborate the use of coarse-grain models for lipid-dependent large-scale oligomerization studies.
In conclusion, our results show that the serotonin 1A receptor is, on an average, less flexible in the presence of membrane cholesterol. In particular, transmembrane helix VII appears to contribute to the reduced flexibility. It is clear that these results represent only a first step toward understanding the subtle effects of cholesterol on receptor conformation. The interaction of cholesterol with GPCRs is highly plastic and differential binding modes are sampled in long time scale simulations . It is therefore possible that the flexibility difference observed in the three models is dependent on this highly plastic and stochastic cholesterol interaction. These results support experimental data which show that membrane cholesterol is crucial in increasing serotonin 1A receptor stability under conditions of thermal deactivation, extreme pH, and proteolytic digestion (Saxena & Chattopadhyay, 2012). In general, membrane cholesterol has been shown to be important in improving the stability of other GPCRs such as the b 2 -adrenergic receptor (Hanson et al., 2008;Yao & Kobilka, 2005;Zocher et al., 2012), and appears to be necessary in crystallization of the receptor . Importantly, the b 2 -adrenergic receptor enjoys $48% amino acid similarity with the serotonin 1A receptor in the transmembrane region (Paila et al., 2011b). Such similarity in the transmembrane region could contribute to cholesterol sensitivity in the function of the b 2 -adrenergic receptor, similar to what was observed with the serotonin 1A receptor (Jafurulla et al., 2014;Paila et al., 2008;Pucadyil & Chattopadhyay, 2004Shrivastava et al., 2010). In agreement with this, we recently showed that adrenergic signaling is enhanced upon cholesterol depletion in cardiac myocytes (Paila et al., 2011a). Our results constitute one of the early reports on a molecular level understanding of GPCR-cholesterol interaction, and could be useful in future crystallization efforts of the serotonin 1A receptor. Figure 8. Characterizing the interaction of POPC at the groove of transmembrane helices I and VII: a schematic representation of a POPC molecule at the site in the (a) presence and (b) absence of cholesterol. The transmembrane helices I, II and VII and helix VIII are represented in blue, magenta, red and light pink, respectively. The cholesterol molecule is shown in green and POPC molecule is colored by atom type (carbon: cyan, phosphate: brown, hydrogen: white, nitrogen: blue and oxygen: red). (c) The order parameters calculated for the carbon atoms of the palmitoyl chain for the two POPC molecules in the presence (black) and absence (red) of cholesterol (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).