The nicotinic acetylcholine receptor (nAChR) and other pentameric ligand-gated ion channels are native to neuronal membranes with an unusual lipid composition. While it is well-established that these receptors can be significantly modulated by lipids, the underlying mechanisms have been primarily studied in model membranes with few lipid species. Here, we use coarse-grained molecular dynamics simulation to probe specific binding of lipids in a complex quasi-neuronal membrane. We ran a total of 50 μs of simulations of a single nAChR in a membrane composed of 36 species of lipids. Competition between multiple lipid species produces a complex distribution. We find that overall, cholesterol selects for concave inter-subunit sites and polyunsaturated fatty acids select for convex M4 sites, while monounsaturated and saturated lipids are unenriched in the nAChR boundary. We propose the “density-threshold affinity” as a metric calculated from continuous density distributions, which reduces to a standard affinity in two-state binding. We find that the density-threshold affinity for M4 weakens with chain rigidity, which suggests that flexible chains may help relax packing defects caused by the conical protein shape. For any site, PE headgroups have the strongest affinity of all phospholipid headgroups, but anionic lipids still yield moderately high affinities for the M4 sites as expected. We observe cooperative effects between anionic headgroups and saturated chains at the M4 site in the inner leaflet. We also analyze affinities for individual anionic headgroups. When combined, these insights may reconcile several apparently contradictory experiments on the role of anionic phospholipids in modulating nAChR.
The nicotinic acetylcholine receptor (nAChR) is a well-studied excitatory pentameric ligand-gated ion channel (pLGIC). nAChRs are found at high density in post-synaptic membranes and the neuromuscular junction in mammals, and the electric organ in Torpedo electric rays. The nAChR is activated by the binding of agonists such as nicotine or acetylcholine to the orthosteric site in the extra-cellular domain (ECD). When post-synaptic nAChRs are activated en-mass, they stimulate an action potential. Thus, nAChRs play a critical role in both cognition and memory1 and neuromuscular function.2,3 nAChR and the greater pLGIC superfamily play various roles in neurological diseases related to inflammation,4–7 addiction,8 chronic pain,9 Alzheimer’s Disease,3,10–12 spinal muscular atrophy,13 schizophrenia,3,14 and neurological autoimmune diseases.15,16
nAChRs are highly sensitive to their local lipid environment. nAChR poorly conducts ions in model phosphatidylcholine (PC)-only membranes but can conduct a current with the addition of cholesterol or anionic lipids,17–25 though too much cholesterol can also cause a loss of function.20,26–28 Functional studies using Xenopus oocytes29–36 require lipid additives such as asolectin20,29–31,33–36 or lipids from synaptic membranes37 to recover native levels of nAChR ion flux. Understanding the mechanism of modulation requires understanding how and where the modulating lipid interacts with the receptor, and these interactions may themselves be dependent upon the rest of the lipid composition.
Mammalian neuronal membranes38–41 have unique compositions compared to other mammalian membranes.42–46 Neuronal membranes are more similar to the membrane of the Torpedo electric ray’s electric organ47,48 than the average mammalian membrane.46 The neuronal membrane38–41 is rich in lipids in which one or both chains are polyunsaturated fatty acids (PUFAs), particularly the n − 6 PUFA arachidonic acid (AA), the n − 3 PUFAs docosahexaenoic acid (DHA), and eicosapentaenoic acid (EPA). These three PUFA’s comprise ∼20%–25% of the acyl chains of neuronal phospholipids and are involved in secondary signaling49,50 and neuronal development.51 PUFAs are linked to a number of neurological diseases and disorders that overlap nAChR related diseases. PUFAs play a role in major depressive and bipolar disorder,50,52–55 schizophrenia,51,53,55–58 and Alzheimer’s Disease.52,59–63
Functional experiments have focused on the role of anionic lipids and cholesterol as modulators of pLGICs18–25,64 (the role of polyunsaturation has received comparatively little attention due to common challenges with oxidation of polyunsaturated chains). Such experiments have been overwhelmingly consistent with a role for direct binding of lipids as a modulatory mechanism. Lipids interacting directly with the protein are called “boundary lipids” and can include lipids in the first few lipidation shells (“annular” or superficial lipids) or those buried more deeply in the protein (“non-annular” or embedded lipids). As for most membrane proteins, it is experimentally challenging to capture the boundary lipid composition of pLGICs because lipids are small molecules that may remain partly fluid even in their bound state. Numerous structures of pLGICs have revealed a conserved arrangement for both the transmembrane domain (TMD) and the extracellular domain (ECD). In the TMD, each subunit has four membrane helices (M1–M4) with the five subunits forming a “star” shape around a central pore [Fig. 1(a)]. The M2 helix lines the pore, the M1 and M3 helices form a middle ring that includes the inter-subunit cavities, and the M4 helices form the tip of the star. Structural methods have resolved potential cholesterol molecules65,66 and phospholipids67–69 bound to subunit interfaces, but crystallographic disorder introduced by lipids typically precludes identification of lipid species. Mass spectrometry has revealed specific binding of anionic lipids, with additional mutagenesis studies suggesting localized sites in the inner leaflet near the M4 helices.70
Molecular dynamics (MD) simulations are particularly useful for visualizing and characterizing microscopic interactions within a fluid system. Given a putative cholesterol or lipid binding mode, atomistic simulations can be used to probe stability of the lipid binding mode. For pentameric channels, such approaches have primarily demonstrated stability of bound cholesterol,71 particularly at inter-subunit sites.65,72 Unfortunately, the timescales currently accessible to fully atomistic simulations are too short for diffusion of lipids, which precludes spontaneous lipid sorting by proteins. While methods that circumvent diffusive barriers are particularly promising73,74 for lipid–protein interactions, they require candidate lipids and binding sites as initial inputs.
Coarse-grained MD simulations use simplified molecular models that can reveal spontaneous lipid sorting, domain formation, and protein partitioning over simulation timescales.75–78 Coarse-grained MD simulations have been used previously to probe interactions of pLGICs with propofol79 as well as spontaneous lipid binding in model membranes.70,80,81 In previous work, we found that nAChR embedded in multiple domain-forming model membranes partitioned to the PUFA-rich liquid-disordered domain,80 rather than to the cholesterol-rich liquid-ordered or “raft” domain that was suggested by cholesterol modulation. We observed that cholesterol still occupies embedded sites on the nAChR TMD, where it is shielded from the liquid-disordered domain. However, native membranes are primarily composed of heteroacidic lipids with two distinct chains, where each chain has a different domain preference; such lipids will naturally destabilize domains. In non-domain-forming model membranes composed of heteroacidic lipids, two classes of fivefold symmetric sites emerged: an inter-subunit site and the M4 site [Fig. 1(b)]. Cholesterol and saturated chains were enriched at the inter-subunit interfaces and n − 3 PUFA acyl chains were enriched around the M4 helices.81 These results were consistent with binding to minimize packing defects: the rigid lipids could fill in the concave regions at the inter-subunit sites, while the flexible chains would easily deform around the “star points” of the M4 helices. Yet, it was not clear whether these same patterns would be upheld in the more complex environment of a native neuronal membrane, which has many more options for minimizing any packing defect.
Neuronal membranes also contain a sizable fraction of anionic lipids in the inner leaflet.38–40 With collaborators in the Cheng laboratory, we recently70 showed that anionic headgroups bind preferentially to the pLGIC Erwinia ligand-gated ion channel (ELIC), when the same acyl chains are used for both headgroups. Through coarse-grained MD, we found specific binding sites for 1-palmitoyl-2-oleoyl phosphatidylglycerol POPG in the inter-subunit sites (inner leaflet); these sites contained basic amino acids that were also implicated through mutagenesis.70 In nAChR, the high-density of basic amino acids are in the M4 site (inner leaflet) rather than the inter-subunit site (inner leaflet), so we would expect a shift for nAChR even in model membranes due purely to the protein sequence. The relative roles of headgroup charge vs acyl chain saturation in driving affinity are unknown.
The use of complex quasi-realistic membranes in coarse-grained MD simulations is increasingly more feasible. In 2014, Ingólfsson et al.46 simulated an “average mammalian” membrane containing 63 lipid species, followed in 2017 by a coarse-grained neuronal membrane.41 Multiple accessible and realistic membranes have been developed for the comparison of protein–lipid interactions between model and quasi-native membranes.45,77,78,82,83 To our knowledge, no such coarse-grained MD simulations using quasi-native membranes have been used with pLGICs.
While the model membranes we used previously are useful for identifying putative sites, they have critical limitations. As stated previously, model membranes typically vary headgroup charge or acyl chain saturation, not both. Model membranes also do not allow for identification of more specific chemical variations within general saturation classes (i.e., n − 3 PUFAs such as DHA vs n − 6 PUFAs such as α-linolenic acid) or like-charged headgroups [PC vs PE, or phosphoserine (PS) vs phosphoinositol (PI)]. For this work, we embed the neuromuscular nAChR84 in a coarse-grained neuronal membrane.41 To test whether the predictions we developed from model membranes hold for native membranes, we develop a new method for quantifying affinities for partially occupied binding sites.
The remainder of this paper is organized as follows: Section II presents our simulation and analysis approach, including introduction of the density-threshold affinity. Section III presents the results and discussion of site selectivity of neutral lipids, followed by a reoriented discussion of the same data that are focused on lipid preferences of individual sites. We then consider selectivity of anionic lipids in the inner leaflet and finally consider the effects of specific headgroup differences. Section IV concludes.
A. Simulation composition
All simulations used the coarse-grained MARTINI 2.285 topology and force-field. nAChR coordinates were based on a cryo-EM structure of the αβγδ muscle-type receptor in a native torpedo membrane (PDB 2BG984). This is a medium resolution structure (4 Å) and was further coarse-grained using the martinize.py script; medium resolution is sufficient for use in coarse-grained simulation, and the native lipid environment of the proteins used to construct 2BG9 is critical for the present study. The secondary, tertiary, and quaternary structures in 2BG9 were preserved via soft backbone restraints during simulation as described below, so any inaccuracies in local residue–residue interactions would not cause instability in the global conformation.
nAChR was embedded in a coarse-grained neuronal membrane based on the work of Ingólfsson et al.41 That neuronal membrane contains phospholipids, sterols, diacylglycerol, and ceramide. Membranes presented in this paper only consider phospholipids and cholesterol for a total of 36 unique lipid species (Table SI 1).
Coarse-grained membranes were built using the MARTINI script insane.py,86 which was also used to embed the coarse-grained nAChR within the membrane. The insane.py script randomly places lipids throughout the inter- and extra-cellular leaflets, and each simulation presented in this manuscript was built separately. All simulation box sizes were 40 × 40 × 35 nm3 with ∼4500 − 5000 lipids and total ∼450 000 beads. This simulation box size corresponds to a protein area density of about 625 µm−2, which is comparable87 to nAChR density in the early stages of synapse formation and substantially lower than that found (104 μm−2)88 at the mature neuromuscular junction. The early stages of clustering are lipid sensitive,27,89,90 while the late stages are driven by linkage of nAChRs via the cytoplasmic peripheral membrane protein rapsyn.91
Molecular dynamics simulations were run using the MARTINI 2.285 force-field and GROMACS92,93 2019.2. All systems used van der Waals (vdW) and electrostatics with reaction-field and a dielectric constant of εr = 15 and electrostatic cutoff length at 1.1 nm. Energy minimization was performed for 1 000 000 steps, but energy minimization tended to concluded after ∼5000–10 000 steps.
Volume and pressure equilibrations were run with isothermal–isochoric (NVT) and isothermal–isobaric (NPT) ensembles, respectively. NVT and NPT simulations used a time step of 15 fs and were run for 0.3 ns using a Berendsen thermostat held at a temperature of 323 K and Berendsen pressure coupling with compressibility set to 3 × 10−5 bar−1 and a pressure coupling constant set to 3.0 ps for the NPT ensemble.
Production simulations were run using a time step of 20 fs for 5 µs for ten replicas. Simulations were conducted in the NPT ensemble by using the velocity rescaling to a temperature of 323 K with a coupling constant set to 1 ps. Semi-isotropic pressure coupling was set to Parrinello–Rahman with compressibility at 3 × 10−5 bar−1 and the pressure coupling constant set to 3.0 ps.
Secondary structure restraints with MARTINI recommendations were constructed by the martinize.py85 script and imposed by GROMACS.92,93 The nAChR conformation was preserved by harmonic bonds between backbone beads separated by less than 0.5 nm and calculated using the ElNeDyn algorithm94 associated with MARTINI85 with a coefficient of 900 kJ mol−1.
C. Calculation of polar density distributions
As in our previous work,70,80,81 the two-dimensional density distribution ρB of the beads within a given lipid species B around the protein was calculated on a polar grid,
where ri = iΔr is the projected distance of the bin center from the protein center, θj = jΔθ is the polar angle associated with bin j, Δr = 10 Å and radians are the bin widths in the radial and angular direction, respectively, and is the time-averaged number of beads of lipid species B found within the bin centered around radius ri and polar angle θj. In order to determine enrichment or depletion, the normalized density is calculated by dividing by the approximate expected density of beads of lipid type B in a random mixture, xBsBNL/⟨L2⟩, where xB is the fraction of lipid B in the membrane, sB is the number of beads in one lipid of species B, NL is the total number of lipids in the system, and ⟨L2⟩ is the average projected box area,
Python software for these calculations are under active development and are available online.95
This expression is approximate because it does not correct for the protein footprint or any undulation-induced deviations of the membrane area. The associated corrections are small compared to the membrane area and would shift the expected density for all species equally, without affecting the comparisons we perform here. For a given lipid species or class, analysis excluded any replicas in which fewer than five lipids of the species/class were in the leaflet at any point in the sampled simulation.
D. Calculation of the density-threshold affinity
Although lipids do occupy clearly detectable hot spots, binding to these sites are not straightforward to describe by a traditional two-state model. Lipids are chains that may partially occupy or fully occupy a site, and they may share a site with another lipid that is partly or fully occupying the site. While the standard affinity can be determined from the probability of single occupancy, the density-threshold affinity is determined from the probability that a site is occupied by more beads than would be expected based on bulk density.
For a given site, consider two probability distributions: the probability Psite(n) of finding n beads within the site and the probability Pbulk(n) of finding n beads within a region of equivalent area in a randomly mixed bulk, respectively. For a lipid that has no affinity for this binding site, we expect Psite(n) = Pbulk(n), while Psite(n) should be right-shifted for a strong affinity and left-shifted in the presence of competition. We calculate the degree of right or left shift by first finding the number of beads npeak corresponding to the peak of the density distribution in the bulk. For a randomly mixed bulk, we expect Psite(n) to be unimodal and close to Gaussian, so npeak should be approximately equal to the mean occupancy ⟨n⟩. Using the mode (npeak) rather than ⟨n⟩ preserves the discretization of the original distribution and is less sensitive to extreme outliers. For a phase-separated bulk with a multimodal distribution, ⟨n⟩ would be more meaningful than npeak. Alternatively, it would be possible to calculate a different affinity with respect to each reference phase by determining npeak using regions that are completely contained within specific phases.
As illustrated in Fig. 1(c), we then integrate Psite on both the left and right side of the threshold npeak to yield P< and P>, respectively,
While all probability distributions in the current manuscript are unimodal (Fig. SI 1), this approach is not limited to unimodal distributions of Psite. Regardless of the shape of the distribution for the bound state, the threshold separates it into two macrostates analogous to those in a classic binary binding model. The free energy difference between the two macrostates is
where R is the gas constant and T is temperature. We term this free energy difference the “density-threshold affinity.” In the special case of binary occupancy,
where KD is the dissociation constant and [L] is the ligand concentration. In a dilute solution, the volume per ligand is typically much larger than the site volume, yielding Pbulk(n) = 1 for n = 0, and vanishes for all n > 0, yielding npeak = 0. Consequently, for this special case, and . Then, Eq. (5) reduces to the classic form for the chemical potential RT ln KD − RT ln[L].
E. Binding site definition and occupancy calculations
We consider four classes of site: each of the two leaflets contains a concave inter-subunit site class and a convex M4 site class. There are five pseudosymmetric sites within each site class for a total of 20 sites [Fig. 1(b)]. The boundaries for each site were drawn to correspond to the localized binding hot spots observed for heteroacidic membranes81 and are non-overlapping. Inter-subunit sites include bins with angular components between the M1 and M3 helices of two adjacent subunits and a radial component satisfying 10 < r ≤ 32 Å. M4 sites include bins with complementary angular components (so that no sites overlap) falling within the M1 and M3 helices of a single subunit and a radial component satisfying 10 < r ≤ 44 Å. For a full description of radial and angular dependencies, please see Table SI 4.
In order to calculate Psite(n), a distribution was taken across frames at 10 ns intervals. For any frame, the beads of a given lipid or chain type were binned onto a fine polar grid with Δr = 4 Å and . The bins falling within the site boundaries were then summed to calculate the occupancy n. This approach allowed for straightforward adjustment of site boundaries if needed without needing to re-bin the whole trajectory. All density-threshold affinities reported in the main manuscript reflect at least 290 samples with n ≥ 1; most values represent at least 10 000 samples. In Tables SI 3.1–3.4, we provide a further breakdown of affinities by specific lipid species and pseudosymmetric sites, and these values may represent fewer samples.
F. Calculation of accessible area
Calculation of Pbulk requires determining the accessible site area in order to calculate the densities in a bulk region of similar area. The area A accessible to the lipids is the difference between the total site area Atot and the area Aex excluded by the protein: A = Atot − Aex. Atot is straightforward to calculate by summing over the areas of the bins i within the site boundaries: Atot = ∑iriΔriΔθi. Calculating Aex is less straightforward, and although there are many possible ways to do this, for self-consistency, we used the same tools from our primary analysis.
In a single lipid membrane, Psite(n) = Pbulk(n) as long as Pbulk(n) is calculated using the proper area A. We exploit this identity to calculate A for each site by running a single nAChR in pure di-palmitoyl phosphatidylcholine (DPPC) for ∼370 ns at 323 K and determining the value of A for each site such that Psite(n) and Pbulk(n) have the same peak. These areas are reported in Table SI 4.
III. RESULTS AND DISCUSSION
A. Effect of acyl chain on site selectivity among neutral lipids
Representative frames from a typical trajectory of boundary lipids are shown in Fig. 2(a), with representative poses shown in Fig. 2(b). In order to quantitatively compare the lipid distributions for the native system to our previous model system, we plotted the enrichment of boundary density relative to bulk density on a two-dimensional polar heat map centered on the protein. This enrichment is shown in Fig. 3(a) for cholesterol and various acyl chains grouped by saturation. Saturated and monounsaturated acyl chains are not significantly depleted or enriched in the boundary of the protein. Regions of cholesterol density are much more localized than in the model membrane [Fig. 3(c)], with pockets of high enrichment very close to the protein and weak depletion in the remainder of the boundary region. Both n − 6 and n − 3 PUFA chains yield fivefold symmetric enrichment around the M4 helices, as observed for n − 3 PUFAs in the model membrane. In the neuronal membrane, however, this enrichment is less well-defined and spreads into the inter-subunit regions. In particular, additional pockets of significant enrichment are apparent in the β − δ subunit interface in the outer leaflet. The overall area of the regions of PUFA-enrichment decreases in the inner leaflet, where n − 3 PUFAs are enriched around M4 helices, but n − 6 PUFA density is not fivefold symmetric and has weak enrichment. Overall, the loss of definition in site boundaries diverges from the well-defined fivefold enrichment for n − 3 PUFAs that we saw in model membranes.81
In order to reduce these distributions to affinities that are more straightforward to interpret, we calculated the density-threshold affinity ΔG for various lipid species as defined in Eq. (5). We organize this information in two different ways: Fig. 5 provides the “lipid’s perspective” and is organized to identify the preferred site for a given lipid type (the lipids’ “site selectivity”), while Fig. 6 provides the “site’s perspective” and is organized to identify the most favorable lipids for a given site (the sites’ “lipid specificity”).
We first consider site selectivity for neutral lipids. Affinities for neutral lipids and cholesterol are shown in Fig. 5(a), where more negative values of ΔG indicate a stronger density-threshold affinity and more positive values indicate more displacement by other lipids. Overall, as shown in Fig. 5(a), saturated lipids have similar density-threshold affinities across all sites, which is consistent with the generally flat distribution observed in Fig. 3. Yet, saturated lipids do yield a slightly stronger affinity for inter-subunit sites, at least in the outer leaflet, which may drive the high amount of saturated enrichment observed at these sites in model membranes. Outer leaflet monounsaturated lipids are slightly more unfavorable in inter-subunit sites than M4 sites, and this difference grows in the inner leaflet.
In contrast to saturated and monounsaturated lipids, PUFAs and cholesterol are highly selective for particular site classes. As shown in Fig. 5(a), neutral PUFAs have significantly stronger affinities for M4 sites than for inner-subunit sites in the same leaflet. Such selectivity is consistent with the PUFA enrichment density in Fig. 3(a), where n − 3 PUFAs can occupy most regions of the TMD but have particularly high levels of enrichment around M4. It is further consistent with our expectations from model membranes [Fig. 3(c)]. Regardless of the site class, PUFAs favor the outer leaflet site over the inner leaflet site, but both sets of M4 sites are more favorable than both sets of inter-subunit sites. Conversely, cholesterol has significantly stronger affinities for inner-subunit sites than for M4 sites, which is also consistent with the enrichment density in Fig. 3(a) and our expectations from model membranes [Fig. 3(c)]. For cholesterol, however, the leaflet is a bigger determinant of affinity than the site; cholesterol has a stronger affinity for either outer leaflet site compared to either inner leaflet site.
B. Lipid preferences of inter-subunit and M4 sites
We now switch perspectives to considering which neutral lipids are most favorable for particular sites. As shown in Figs. 6(a) and 6(b), inter-subunit sites in both leaflets prefer cholesterol to phospholipids, which is expected based on the results from model membranes. Upon visual inspection, this result may appear to diverge from the cholesterol polar density plots in neuronal membranes [Fig. 3(a)]. The present results show that while the overall footprint of cholesterol enrichment in Fig. 3(a) is small, this small region actually reflects a tight and persistently occupied binding site. The highly right-shifted distributions for cholesterol are shown in Fig. S1.
PUFA chains yield affinities for the inter-subunit site that are approximately >0.5 kcal/mol stronger than saturated lipid affinities [Figs. 6(a) and 6(b)], which was unexpected based on the results from model membranes but is consistent with the corresponding enrichment density in Fig. 3(a). More generally, neutral phospholipid affinities for inter-subunit sites obey the following trend, from strongest to weakest: n − 3 > n − 6 > saturated > monounsaturated. Thus, even though PUFA chains prefer M4 sites to inter-subunit sites and saturated chains prefer inter-subunit sites to M4 sites, PUFAs have a stronger affinity for either site type than do saturated lipids.
For inter-subunit sites, monounsaturated lipids have the weakest affinities (>0.5 kcal/mol), which may reflect a limited number of ways to pack the single kink of a monounsaturated chain in this concave site. In contrast, cholesterol and PUFAs are either small or highly flexible and may more easily pack across multiple sites. Saturated chains may pack parallel to the protein surface in these sites [Fig. 2(b)].
As shown in Figs. 6(c) and 6(d), M4 sites in both leaflets have the strongest affinity for n − 3 PUFAs, and affinity weakens as acyl chain rigidity increases; from strongest to weakest, the phospholipid affinities follow: n − 3 > n − 6 > monounsaturated > saturated. This is consistent with a role for PUFAs in minimizing unfavorable membrane deformations caused by the pLGIC’s conical-star shape.96–101 Surprisingly, cholesterol had a stronger affinity for M4 sites than any acyl chains other than n − 3 PUFAs. Cholesterol is rigid and small and has asymmetric sides (rough and smooth), which potentially allows it to embed between helices and compete with n − 3 PUFAs for binding. Any cholesterol bound within the grooves of the subunit interface (as hypothesized based on atomistic simulations71 and observed in β subunits of nAChR using coarse-grained simulations80) will also get counted within the M4 site.
All pLGICs display a high degree of structural conservation in the TMD and ECD, both across proteins and across subunits of the same protein, but sequences can vary significantly. Affinities for each of the fivefold pseudosymmetric sites within each site class are given in Tables SI 3.1–3.4, along with significance tests. Cholesterol selects for specific subunits within each of the site classes, which is consistent with its rigid structure. Three of the four site classes (inter-subunit sites in both leaflets and outer leaflet M4 sites) include deep grooves in the protein, which could contribute to sequence-dependent affinities. Saturated, n − 3, and n − 6 chains reveal significant preferences within these site classes. Affinities for individual lipid species have a higher sampling error than aggregates based on headgroup or chain saturation, and most variations at the species level do not reach statistical significance. PUPE (a PE headgroup with one saturated chain and one long-chain n − 3 polyunsaturated chain) is an exception; it reveals selectivity among all three grooved site classes.
C. Effect of headgroup charge on affinity depends on leaflet and binding site
Figure 3(b) compares the density enrichment for anionic headgroups with that of neutral headgroups. Data are shown for the inner leaflet only because anionic lipids are not present in the outer leaflet at the start of simulations and few anionic lipids flip flop to the outer leaflet. In the inner leaflet, the anionic lipids are expected to select for sites that are lined with basic amino acids, which are in different locations depending on subunits (Fig. SI 2) As shown in Fig. 4, anionic lipids are generally enriched around the M3/M4 helices for the αγ, γ, δ, and β subunits. Anionic lipids are enriched at inter-subunit sites and around M4 sites for all subunits but the α subunits. Non-α nAChR subunits have basic amino acids closer to M4 helices, as shown in Fig. SI 2A. We incorporate data from all five pseudo-symmetric sites to obtain the density-threshold affinities reported in Fig. 5(b), which suggest that anionic lipids have significantly stronger affinities for M4 sites on average. The average anionic affinity difference between inter-subunit and M4 sites is ∼−1.0 kcal/mol, as shown in Figs. 5 and 6(d) and SI Table 2. Although the magnitude of the affinity difference varies with acyl chain saturation, the sign is unchanged.
We now switch again to the “site perspective” to compare whether inner leaflet sites would prefer occupancy by anionic or neutral lipids. As shown in Fig. 6(c), lipid affinity values for inter-subunit sites are either insensitive to charge (saturated or n − 6 PUFA chains) or weaker for anionic lipids by at least 0.5 kcal/mol (monounsaturated and n − 3 PUFA chains). In comparison, at the M4 site, saturated chains in anionic lipids have significantly stronger affinities than those in neutral lipids (a difference of ∼0.5 kcal/mol). All other lipid chains attached to anionic headgroups have weaker affinities for the M4 site. The clear trend observed in neutral lipids (stronger affinities for more flexible acyl chains) is thus broken in anionic lipids because saturated anionic lipids are so favorable.
In summary, we observe that binding sites have clear preferences for a particular headgroup charge and acyl-chain saturation, suggesting nAChR lipid occupancy is driven by two steps, a “coarse-sorting” by headgroups and then “fine-sorting” by acyl chains. Zwitterionic headgroups will occupy the nAChR boundary, but acyl chain saturation drives the hierarchy within individual site classes. Anionic lipids diverge from this pattern at the inner M4 site, which selects for anionic lipids independent of saturation.
D. Role of individual lipid headgroups in determining affinity
Neutral and anionic are bulk terms that categorize numerous lipid head-groups by charge. To understand the role of the chemical distinctions between headgroups of like charge, we broke the headgroup affinities down by headgroup species in Tables I and II. In the outer leaflet, lipids contain a mixture of PE and PC headgroups. The small neutral PE headgroup has the strongest affinity across all headgroups for both inter-subunit and M4 sites, −0.2 ± 0.3 and −1.1 ± 0.2 kcal/mol, respectively. The larger neutral PC headgroups are weaker than PE by ∼>0.5 kcal/mol. In living cells, as in this neuronal membrane, PUFAs are more frequently tethered to PE than to PC or SM,38–41,45 so it is possible that this affinity simply reflects the high affinity of PUFA chains. However, even for identical chains, both experimental and simulation data80 suggest stronger PE–ELIC than PC–ELIC interactions.
|.||Inter-subunit sites .||M4 sites .|
|.||ΔG (kcal/mol) .||ΔG (kcal/mol) .|
|PE||−0.2 ± 0.3||−1.7 ± 0.2|
|PC||1.6 ± 0.3||0.9 ± 0.1|
|.||Inter-subunit sites .||M4 sites .|
|.||ΔG (kcal/mol) .||ΔG (kcal/mol) .|
|PE||−0.2 ± 0.3||−1.7 ± 0.2|
|PC||1.6 ± 0.3||0.9 ± 0.1|
|.||Inner inter-subunit .||Inner M4 .|
|.||sites ΔG (kcal/mol) .||sites ΔG (kcal/mol) .|
|PE||0.3 ± 0.2||−0.2 ± 0.1|
|PI||1.4 ± 0.3||0.3 ± 0.1|
|PS||1.4 ± 0.2||0.5 ± 0.2|
|PC||1.3 ± 0.3||0.8 ± 0.1|
|PIP3||3.1 ± 0.5||2.4 ± 0.4|
|PIP2||2.4 ± 0.3||1.3 ± 0.4|
|PIP1||2.2 ± 0.3||1.3 ± 0.4|
|PA||2.8 ± 0.3||1.9 ± 0.4|
|.||Inner inter-subunit .||Inner M4 .|
|.||sites ΔG (kcal/mol) .||sites ΔG (kcal/mol) .|
|PE||0.3 ± 0.2||−0.2 ± 0.1|
|PI||1.4 ± 0.3||0.3 ± 0.1|
|PS||1.4 ± 0.2||0.5 ± 0.2|
|PC||1.3 ± 0.3||0.8 ± 0.1|
|PIP3||3.1 ± 0.5||2.4 ± 0.4|
|PIP2||2.4 ± 0.3||1.3 ± 0.4|
|PIP1||2.2 ± 0.3||1.3 ± 0.4|
|PA||2.8 ± 0.3||1.9 ± 0.4|
Table II shows specific headgroup affinities in the inner leaflet. As in the outer leaflet, lipids with PE headgroups still have the strongest affinity of all lipids, but in the inner leaflet, we are also able to distinguish affinities for anionic species. For the inter-subunit site, PI, PS, and PC have similar affinities (within statistical error) and have significantly stronger affinities for these sites than the phosphoinositides (PIPS) PIP1, PIP2, PIP3, which have a significantly stronger affinity than phosphatidic acid (PA). Thus, from strongest to weakest, PE > PI ∼ PS ∼ PC ≫ PIP1 ∼ PIP2 ∼ PIP3 ≫ PA for the inter-subunit site. In contrast, at the M4 site, more significant differences among moderate affinity headgroups emerge. PI has significantly stronger affinity than PS (a difference of 0.3 ± 0.1 kcal/mol), and PS has a significantly stronger affinity than PC (a difference of 0.2 ± 0.1 kcal/mol). From strongest to weakest, PE > PI > PS > PC ≫ PIP1 ∼ PIP2 ∼ PIP3 ∼ PA for the M4 site.
Using coarse-grained simulations of the nAChR within a quasi-neuronal membrane containing over 30 lipid species, we have observed spontaneous lipid binding and quantified lipid specificity for two types of sites in the protein TMD. These two site classes represent the most concave (inter-subunit site) and convex (M4 site) portions of the star-shaped nAChR and were initially observed as “hot spots” in our previous simulations70,81 of model membranes. Compared to classic ligand binding sites, these sites are superficial and have a large volume. The “ligands” occupying them are also non-traditional: lipids are flexible chain molecules that may only partially occupy the site and are likely to share the site with other partially occupying ligands. While our laboratory has developed promising alchemical approaches73 for calculating traditional affinities of atomistic lipids for more highly localized, well-defined sites, these hot spots required a different approach. Here, we have proposed a softer “density-threshold affinity” for characterizing these affinities from spontaneous, unbiased coarse-grained simulations. While we restrict the use of this method here to nAChR, it should be straightforward to extend to any other transmembrane proteins with detectable regions of density enrichment.
Our results are summarized graphically in Fig. 7. Based on our results from model membranes, we had hypothesized that PUFAs would select for the convex M4 sites and that raft-forming lipids such as cholesterol and saturated lipids would select for the concave inter-subunit sites. Overall, our results were consistent with this expectation. Yet, although lipids containing PUFAs do prefer the M4 site to the inter-subunit site, their affinity for even the inter-subunit sites are stronger than that of all other phospholipids. This result underscores the reliable partitioning of nAChR to PUFA-rich, liquid-disordered domains that we observed in homoacidic, domain-forming membranes80 and suggests that PUFAs may have been absent from the inter-subunit site in heteroacidic membranes81 because of the constraints of the lipid topology. In the latter simulations, all lipids contained one saturated chain and one PUFA chain, so binding of the PUFA chain to its preferred M4 site requires the saturated chain to find the most favorable location nearby (in the inter-subunit site) and may block binding of other PUFA chains to that site. These constraints are relaxed in the native neuronal membrane, which has a more diverse lipid composition with multiple different chain pairings; about 6% of the phospholipids in our simulated membranes contain no saturated chain at all. Nonetheless, our previous results81 using simplified binary heteroacidic/cholesterol membranes played a key role in identifying the natural site boundaries.
As expected, within each leaflet, cholesterol has the strongest affinity for the inter-subunit sites, although the affinity of cholesterol for the M4 sites was second only to that of n − 3 PUFAs. When combined, these results are consistent with an overwhelming amount of evidence spanning four decades that suggests direct interactions between cholesterol and nAChRs, regardless of the phospholipid composition of the membrane. One surprise for cholesterol was the role of the leaflet in determining affinity: cholesterol has a stronger affinity for either outer leaflet site compared to either inner leaflet site. This result may reflect competition with anionic saturated lipids in the inner leaflet, which would be consistent with multiple experiments,102–105 suggesting that anionic lipids can partially or fully compensate for a loss of cholesterol. This result is also consistent with models of cholesterol embedded71 in the outer TMD (which has numerous gaps in the amino acid density) but not the inner TMD.
Based on our results using ELIC,70 we had expected that anionic lipids would select for sites on the inner leaflet lined with basic residues. In the homomeric ELIC, these residues are symmetrically arranged, while in the heteromeric nAChR, they vary by subunit (Fig. SI 2A), with the M4 site containing the most such residues on most subunits. The present results support that expectation: anionic lipids have a stronger affinity for M4 than inter-subunit sites.
For both outer and inner leaflets, lipids with a PE headgroup have stronger affinities than those with a PC head group. This result is consistent with our previous coarse-grained simulations70,80 as well as experimental data.70 The most straightforward explanation might involve the smaller size of the PE headgroup, which could facilitate acyl chains getting tangled in the protein TMD. These subtle differences in steric hindrance between PE and PC are not captured in MARTINI 2.2, however. The effective size of the PE headgroup in MARTINI is still smaller than that of PC because it has a stronger attraction to several species of “nonpolar” and charged beads in the MARTINI force-field, which mimics PE’s ability to donate hydrogen bonds. The increased affinity of PE for the protein may have a similar origin because the protein has many hydrogen bond acceptors. Further investigation of the role of hydrogen bonding between the PE headgroup and the protein will require an atomistic force-field.
Among anionic lipids in the inner leaflet, regardless of the site, PS and PI have an affinity greater than or equal to PC and much greater than the other anionic lipid headgroups (PIP1, PIP2, PIP3, and PA). The lipid headgroups PS and PI both have a charge of −1, while PA in the MARTINI force-field85 carries a charge of −2, and PIP1, PIP2, and PIP3 have charges of −3, −5, and −7, respectively. These results suggest that the inner leaflet sites select for monoanionic headgroups, while multianionic headgroups are highly unfavorable. Due to the limitations of the coarse-grained model, future atomistic calculations are required to validate and understand the apparent preference of the M4 site for PI over PS.
All results presented here relied on a single coarse-grained model. Although a specific group of predictions that we previously made using MARTINI were validated with experiments,70 all coarse-grained force-fields have limitations. Pitfalls of the MARTINI 2.2 model have been discussed extensively.106 Our particular study investigates the ligand-like properties of lipids, and ligand binding of any kind suffers from a coarse-grained treatment due to significant smoothing of the molecular surfaces. Furthermore, many of the effects we investigate here concern lipid chain flexibility, and the chain loses degrees of freedom through the 4-to-1 mapping of the MARTINI model. Consequently, we expect that coarse-graining should generally reduce selectivity among neutral lipids and that the equivalent calculations using an atomistic model could yield a much wider range of binding affinities. Here, we have identified several lipid–site pairs for testing by methods for binding affinities using atomistic lipids, such as alchemical free energy perturbation.73
The present results highlight the utility of model membranes for developing hypotheses of specific lipid–protein interactions and the need to test those hypotheses within more complex native membranes. They could also be tested and aid in interpretation of experiments carried out in more complex membranes. For instance, we would expect that mutations of the basic residues facing the inner leaflet would reduce binding of saturated phospholipids with anionic headgroups, which would be replaced with bound cholesterol. We would also predict that if PUFAs cause gain of function via binding to the inter-subunit site, this gain would be enhanced by replacing some heteroacidic lipids with homoacidic lipids while keeping the total fraction of PUFA chains constant. In general, the present results provide valuable insight into how to predict local lipid competition, which is one of the primary challenges of interpreting experiments in complex membranes.
See the supplementary material for neuronal membrane composition, lipid probability distributions, additional binding affinity tables, and details of the site boundaries.
G.B. and L.S. were supported by the Busch Biomedical Foundation. The authors acknowledge the Office of Advanced Research Computing (OARC) at Rutgers, The State University of New Jersey for providing access to the Amarel and Caliburn clusters and associated research computing resources that have contributed to the results reported here. We are grateful to Dr. Jérôme Hénin for helpful input and suggestions.
The data that support the findings of this study are available from the corresponding author upon reasonable request. Scripts for polar density analysis and plotting scripts can be found on github: https://github.com/BranniganLab/densitymap.