Untangling direct and domain-mediated interactions between nicotinic acetylcholine receptors in DHA-rich membranes

TL;DR

This study uses molecular dynamics simulations to explore how lipid domains influence nicotinic acetylcholine receptor clustering and lipid preferences in membranes, revealing that DHA-rich domains promote receptor dimerization and affect lipid access.

Contribution

It provides new insights into the role of lipid domains and DHA in nAChR clustering and membrane organization, using detailed simulation approaches.

Findings

Dimerization of nAChRs is more likely in domain-forming membranes.

nAChRs preferentially partition into DHA-rich liquid-disordered domains.

Protein concentration influences lipid boundary composition and access.

Abstract

At the neuromuscular junction (NMJ), the nicotinic acetylcholine receptor (nAChR) self-associates to give rise to rapid muscle movement. While lipid domains have maintained nAChR aggregates in-vitro, their specific roles in nAChR clustering are currently unknown. In the present study, we carried out coarse-grained molecular dynamics simulations (CG-MD) of 1-4 nAChR molecules in two membrane environments: One mixture containing domain-forming, homoacidic lipids, and a second mixture consisting of heteroacidic lipids. Spontaneous dimerization of nAChRs was up to ten times more likely in domain-forming membranes; however, the effect was not significant in four-protein systems, suggesting that lipid domains are less critical to nAChR oligomerization when protein concentration is higher. With regard to lipid preferences, nAChRs consistently partitioned into liquid-disordered domains occupied by the omega-3 ($\omega$-3) fatty acid, Docosahexaenoic acid (DHA); enrichment of DHA boundary lipids increased with protein concentration, particularly in homoacidic membranes. This result suggests dimer formation blocks access of saturated chains and cholesterol, but not polyunsaturated chains, to boundary lipid sites.

Tables

Table 1. Table 1: Composition of nAChR Boundary Lipid chains in domain forming (homoacidic) and non-domain forming (heteroacidic) membranes. Embedded and Annular chains are determined as described in Methods. Averages are across systems, and across proteins in multi-protein systems. Each protein is treated as a separate replica (n=30) and standard errors are shown in parentheses.

	mean embedded fraction		mean annular fraction		bulk
	homoacid	heteroacid	homoacid	heteroacid
cholesterol	0.29 (0.02)	0.42 (0.04)	0.23 (0.02)	0.31 (0.01)	0.30
PUFA	0.61 (0.01)	0.39 (0.04)	0.52 (0.02)	0.32 (0.01)	0.35
saturated	0.10 (0.05)	0.19 (0.01)	0.25 (0.02)	0.37 (0.01)	0.35

Equations

$$f_{\mathrm{emb}}^{\alpha }$$

$$f_{\mathrm{ann}}^{\alpha }$$

$${\rho }_B(r_i,{\theta }_j)=\frac{⟨n_B(r_i,{\theta }_j)⟩}{r_i\Delta r\Delta \theta }$$

$${\tilde{\rho }}_B(r_i,{\theta }_j)=\frac{{\rho }_B(r_i,{\theta }_j)}{x_B{s}_B{N}_L/⟨L^2⟩}$$

$$F_{x,y}=\frac{25}{n_d}\left(n_{x,y}+n_{y,x}\right),$$

Full text

††journal: Journal of Membrane Biology

1 Introduction

The muscle-derived nicotinic acetylcholine receptor (nAChR) (PDB 2BG9) (Unwin, 2005) is the most abundant neurotransmitter receptor at the neuromuscular junction (NMJ) in most vertebrates, including humans (Albuquerque et al., 2009). Within the postsynaptic membrane, nAChRs cluster in high densities (10,000 per $\mu$m2) to properly activate the skeletal muscle (Ramarao and Cohen, 1998; Breckenridge et al., 1972). As a major transmembrane protein, nAChR depends upon a highly specific lipid environment to maintain functionality. Lipids influence nAChR activity by affecting both function and organization. It is essential to understand how changes in lipid environment impact nAChR’s structure and activity, given that lipids can change in response to aging and disease (Yadav and Tiwari, 2014), and also vary across tissue and organism.

Over the past few decades, considerable progress has been made in uncovering lipid sensitivities associated with nAChR (Criado et al., 1982). In a majority of experiments, researchers have prioritized studying cholesterol over other membrane lipids. Early studies revealed that, when reconstituted into membrane mixtures, nAChR failed to conduct cations across the lipid bilayer unless cholesterol was present (Fong and McNamee, 1986; Sunshine and MG, 1992; Butler and McNamee, 1993; Fong and McNamee, 1987; daCosta et al., 2001). More recently, researchers have examined the effects of membrane dynamics on nAChR organization (Baenziger et al., 2015; Brusés et al., 2001; Marchand et al., 2002; Oshikawa et al., 2003; Pato et al., 2008; Zhu et al., 2006; Baenziger et al., 2017; Barrantes, 2007; Barrantes et al., 2000; Bermudez et al., 2010; Barrantes et al., 2010; Perillo et al., 2016; Wenz and Barrantes, 2005; Borroni et al., 2016; Unwin, 2017). In-vitro studies (Barrantes, 2007; Barrantes et al., 2010) indicate that nAChRs form larger aggregates upon cholesterol depletion. Experimental evidence suggest that cholesterol-rich lipid domains, known as lipid rafts, facilitate clustering of nAChRs (Campagna and Fallon, 2006; Marchand et al., 2002; Pato et al., 2008). More specifically, after disrupting lipid raft formation, Zhu et al. observed a significant loss of nAChR clusters in-vitro (Zhu et al., 2006). In the mature neuromuscular membrane, nAChRs are linked by the intracellular anchoring protein, rapsyn, which bridges receptors together at their bases (Zuber and Unwin, 2013). According to fluorescent studies (Marchand et al., 2002), lipid rafts mediate the association between rapsyn and neighboring nAChR molecules. In the mid-2000’s, Willmann et al. and Stetzkowski-Marden et al. proposed that lipid rafts can stabilize receptor networks and may even provide a localized environment for nAChR (Willmann et al., 2006; Stetzkowski-Marden et al., 2006).

Domain formation occurs when a membrane is comprised of at least three lipid types: cholesterol, unsaturated fatty acids, and a molecule that interacts closely with cholesterol such as saturated fatty acids or sphingomyelin (Feller, 2008; Yeagle, 2016). The membrane separates into at least two domains, a liquid-ordered or “raft” phase containing cholesterol and saturated lipids/sphingomyelin, and a liquid-disordered phase containing unsaturated lipids. Polyunsaturated phospholipids make these domains more well-defined (Levental et al., 2016). Some results from atomistic simulations (Sodt et al., 2014; Iyer et al., 2018) indicate that cholesterol may preferentially partition to the boundary between liquid-ordered phases and liquid-disordered phases composed of monounsaturated lipids. In neuromuscular membranes, intrinsic domain formation is dependent upon several lipid species, including the widely influential omega-3 ($\omega$-3) fatty acids. One $\omega$-3 in particular, Docosahexaenoic acid (DHA), is prevalent in the native neuromuscular membrane and is strongly associated with flexible and well-defined domains (Turk and Chapkin, 2013; Shaikh et al., 2004). Additionally, DHA is a major contributor to brain functioning, motor activity, and cardiac health; however, its specific effects on neuromuscular health are poorly understood (Lavandera et al., 2017; Wassall and Stillwell, 2008; Georgieva et al., 2015).

Characterizing such complex lipid-protein interactions requires detailed information on each molecular structure. The nAChR is a member of a family of ion channels, known as pentameric ligand-gated ion channels (pLGICs), which has a number of recent structures (Laverty et al., 2019, 2017; Masiulis et al., 2019; Althoff et al., 2014; Hibbs and Gouaux, 2011; Morales-Perez et al., 2016a; Baenziger and Corringer, 2011; Corringer et al., 2012; Nemecz et al., 2016; Prevost et al., 2012; Sauguet et al., 2014); however, only single structures have been determined, rather than dimers. pLGICs are composed of five subunits; each subunit is composed of an extracellular domain for ligand binding, a transmembrane domain (TMD), and an intracellular domain (ICD) (Figure 1A). The TMD is the embedded portion of pLGICs that interacts most often with surrounding membrane lipids; this region of pLGICs is composed of four alpha-helices, M1-M4, with the innermost, M2 helices forming the ion pore. The ICD, or cytoplasmic domain, is highly disordered, making it challenging to obtain information on its structure. Currently, there are two structures of nAChR pentamers resolved from x-ray crystallography (neuronal type) and cryo-electron microscopy (muscle-type) (Morales-Perez et al., 2016b; Unwin, 2005).

Other closely related channels have been resolved with bound lipids. In 2017, the inhibitory neurotransmitter receptor, the gamma-aminobutyric acid (GABAA) receptor, was crystallized with a cholesterol molecule protruding between its M1 and M3 helices (Laverty et al., 2017), similar to our previous prediction (Hénin et al., 2014). Interestingly, x-ray structures of the glutamate-gated chloride channel (GluCl) revealed embedded POPC phospholipids in the same location (Althoff et al., 2014). Together, these data provide evidence for lipid-based modulation of pLGICs, which can potentially be applied to the gating of nAChR.

Structural information alone, however, is insufficient for answering questions about native cell membranes. For one, most structures are obtained under artificial conditions, using detergents or nanodisks with non-native lipid composition. Additionally, current structures do not represent oligomers, particularly in a liquid state. Fully atomistic molecular dynamics (MD) simulations (Brannigan et al., 2008; Cheng et al., 2009; Hénin et al., 2014; Carswell et al., 2015) have complemented experimental investigations of pLGICs interacting with lipids, but they are limited in their ability to capture domain formation since atomistic lipids cannot diffuse over typical simulation microsecond time scales (Ingólfsson et al., 2014; Bond and Sansom, 2006; Parton et al., 2013; Goose and Sansom, 2013; Scott et al., 2008). Furthermore, it is unfeasible to incorporate multiple nAChRs in simulations with atomistic resolution. Coarse-grained MD (CG-MD) simulations are run over longer length and time scales, making them suitable for exploring complex model membranes. CG-MD is widely applied in simulations of both lipid-protein binding and domain formation (Bond and Sansom, 2006; Scott et al., 2008; Parton et al., 2013; Goose and Sansom, 2013). Additionally, CG-MD can capture large-scale membrane phenomenon such as protein self-assembly and lipid-mediated oligomerization (Gahbauer and Böckmann, 2016; Baaden and Marrink, 2013).

Recently, we (Sharp et al., 2019) conducted CG-MD simulations of a single nAChR from the electric ray Torpedo in mixed membranes. Contrary to expectations, nAChR consistently preferred a local lipid environment rich in PUFAs rather than cholesterol, especially long-chained $\omega$-3s, such as DHA. While cholesterol occupied the transmembrane gaps of nAChR, PUFAs were even more likely to be embedded, regardless of their phospholipid headgroup.

The present study adopts a similar approach, with a particular focus on nAChR-associated clustering. Through molecular dynamics simulations, we investigate nAChR lipid preferences and clustering behavior in membranes with and without domains. For this study, we tested three major hypotheses: 1) Membrane organization affects nAChR boundary lipid specificity: when PUFA chains are prevented from forming PUFA rich domains, their prevalence among nAChR boundary lipids will be significantly reduced. 2) Domain formation will indirectly facilitate the clustering of nAChRs, by inducing partitioning preferences and restricting diffusion within the membrane. 3) Within a dimer, we will observe sequence preferences in facing subunits.

2 Methods

2.1 System setup

In order to isolate the role of lipid domain formation on nAChR-based lipid sorting and clustering, our simulations compared two different membranes with distinct phospholipid topology. Each membrane had 30% cholesterol and 70% phospholipids. The phospholipids all had phosphatidylcholine (PC) headgroups, and half the total number of acyl chains were saturated (16:0) chains, while the other half were polyunsaturated (22:6) chains. For homoacid membranes, all phospholipids had either two polyunsaturated acyl chains (didocosahexaenoylphosphatidylcholine or dDHA-PC) or two saturated chains (dipalmitoylphosphatidylcholine or DPPC), facilitating domain-formation. In the heteroacid membranes, all phospholipids were hybrid lipids with one polyunsaturated acyl chain and one saturated chain (1-palmitoyl- 2-docosahexaenoyl- phosphatidylcholine or PDPC), which are topologically unable to separate into PUFA-rich and PUFA-poor domains.

Lipids and proteins were modeled using the coarse-grained (CG) Martini 2.2 force field (Marrink et al., 2007a). Systems included between 1-4 nAChR molecules, derived from the Torpedo electric organ (Unwin, 2005) (PDB 2BG9). This structure is the only pLGIC structure obtained in a native membrane.

We converted protein structures into CG models using the Martini script ”martinize.py”, mapping four non-hydrogen atoms to one CG interaction. We constructed and assembled our protein-bilayer systems using the Martini script, ”insane.py”, using a box sizes of 29x29x21 nm3 for one-to-two proteins, 40x40x20 nm3 for three proteins, and 44x44x22 nm3 for four protein systems, respectively (Marrink et al., 2007b). Initially, the proteins were in a circle of about 13 nm. The receptors were evenly spaced, and their $\delta$ subunits were facing the same direction. Once simulations started, nAChRs shifted from their initial orientations. We ran 24 CG-MD simulations containing 1-4 nAChRs (3 replicas per system).

2.2 Simulation details

Simulations were run using the Martini 2.2 force field parameters and the Gromacs 5.1.2 simulation package, (Marrink et al., 2007b; Pronk et al., 2013) as in our previous work, (Sharp et al., 2019). Each simulation consisted of two steps: energy minimization and molecular dynamics. For each system, we ran two consecutive energy minimizations for 10,000 steps. Harmonic restraints between backbone atoms were imposed to preserve nAChR conformation. More specifically, we applied an elastic force constant of 750 kJ/mol and set lower and upper bounds on the bond with using a bond length of 0.7 nm (Marrink et al., 2007b; Pronk et al., 2013). The molecular dynamics simulations ran for 10-20 $\mu$s at a 0.025 ps time-step. Simulation temperature and pressure were kept constant at values of 323 K and a reference pressure of 1 bar. The isotropic pressure coupling compressibility constant was maintained at 3.0 $\times{10^{-5}}$ bar*-1*.

2.3 Analysis

For a given nAChR, $n_{\mathrm{emb}}^{\alpha}$ is the total number of embedded lipids of lipid species $\alpha$, where embedded lipids satisfy the following criteria: the headgroup (in the case of cholesterol) or the terminal bead on the acyl chains (in the case of phospholipids) are within 10 Å of the M2 helices. Visual inspection indicated any lipids outside of this range were also outside of the TMD bundle.

Similarly, $n_{\mathrm{ann}}^{\alpha}$ is the total number of annular lipids of lipid species $\alpha$, where annular lipids satisfy the following criteria: headgroup or the terminal bead on the acyl chains is between 10 Å and 35 Å of the M2 helices. This range of lipids corresponded to the third lipidation shell around nAChR.

(For phospholipids, each acyl chain was counted separately as a half-lipid, so it was possible for e.g. the sn-1 chain to be embedded and the sn-2 chain to be annular. )

Boundary lipid fractions for a given species $\alpha$ are defined as

[TABLE]

where $n_{\mathrm{emb}}$ and $n_{\mathrm{ann}}$ are the total number of embedded and annular lipids, respectively.

Two dimensional density distributions of boundary acyl chain species (B), $\tilde{\rho}_{B}(r_{i},\theta_{j})$ were calculated, as a function of radius, $r_{i}$, and angle $\theta_{j}$ projected onto the membrane.

[TABLE]

where $r_{i}=i\Delta r$ is the distance from the origin to the center of bin i, $\Delta r$ is the radial bin width, and $\Delta\theta$ is the angular bin width. $\left\langle n_{B}(r_{i},\theta_{j})\right\rangle$ is the time-averaged number of acyl chain beads of species B within bin i.

To quantify enrichment or depletion of a given acyl chain with respect to random distribution, the normalized density, $\tilde{\rho}_{B}(r_{i},\theta_{j})$, was calculated.

[TABLE]

where $s_{B}$ is the number of beads of lipid species B, $N_{L}$ is total number of lipids in a system, and $\langle L^{2}\rangle$ is the average projected box area, and $x_{B}$ is lipid B concentration. The expression does not take into consideration the protein footprint or undulations present within the system, and as such is an approximation.

The radial distribution function $g(r)$ of multiple nAChRs was calculated using the three dimensional distance $r$ between the centers of mass of the receptors for each of the 8000 frames, evenly distributed through the simulation. It was derived from the distribution of pairwise separations, $P(r)$, by dividing by the expected separations in a random distribution : $g(r)=\frac{P(r)~{}dr}{r~{}dr}\times\int_{0}^{R}rdr$. The height of the $g(r)$ peak corresponds to the amount of enrichment in probability for a given distance; for example, if g(10 nm) = 50, it is 50 times more probable that two monomers will be 10 nm apart than expected in a random distribution.

The total number of observed dimers $n_{d}$ is given by the number of pairs (across all analyzed proteins) where $r<100$ Å. For each observed dimer, the closest subunit pair was determined, and the enrichment calculated as $F_{x,y}$:

[TABLE]

where $n_{x,y}$ represents the number of dimers in which subunits $x$ and $y$ form the closest pair, and $1/25$ is the expectation for a random distribution.

We visualized and imaged all simulation results using Visual Molecular Dynamics (VMD) (Humphrey et al., 1996).

3 Results

3.1 nAChR boundary lipid preferences

In order to determine the significance of domain formation on enrichment of polyunsaturated acyl chains among boundary lipids, we calculated the fraction of embedded and annular boundary lipids in two types of membranes with equal amounts of cholesterol and saturated and polyunsaturated acyl chains. In the domain-forming membranes, all phospholipids had either two polyunsaturated acyl chains or two saturated chains (homoacids), while in the non-domain forming membranes, all phospholipids had one polyunsaturated acyl chain and one saturated chain (heteroacids).

Distributions for $f_{emb}$ and $f_{ann}$ are shown in Figure 3 in heteroacids compositions without domains and in Figure 4 for homoacids compositions with domains, while mean values are given in Table 1. For heteroacids, PUFA chains are slightly enriched among embedded chains and slightly depleted among annular chains, while the reverse trend is observed for saturated lipids. Saturated chains are significantly depleted from the embedded lipids, but due to the topology of these lipids, a fully embedded polyunsaturated acyl chain must also contribute a saturated chain to the protein annulus. Distributions were remarkably consistent across simulations containing 1,2,3, or 4 proteins.

Figure 4 shows the results of the same analysis for homoacidic membranes. Replacing heteroacids with homoacids substantially increases the fraction of polyunsaturated chains among both embedded and annular lipids, with the peak of the distribution shifting further to the right as more proteins are included (Figure 4).

Together, these results are consistent with our previous (Sharp et al., 2019) observation that polyunsaturated chains can displace cholesterol from embedded sites. For heteroacids, embedding a PUFA constrained the linked saturated chain to the nAChR annulus. In this case, cholesterol (which introduces no such constraint) was enriched among embedded lipids. For homoacids, embedded PUFAs constrained a corresponding PUFA chain to the nAChR annulus, and here cholesterol was not enriched among embedded lipids.

The fraction of embedded chains that were polyunsaturated increased as more proteins were added to the homoacidic membrane. This could be consistent with each nAChR monomer in an oligomer blocking access to embedded lipid binding sites in the other monomers: flexible PUFA chains have multiple routes to access an embedded site, and rigid lipids have only one or two.

The two dimensional density distribution of each lipid species from the center of a single nAChR is shown in Figure 5. For all three lipids, there was a five-fold symmetry of densities surrounding nAChR, with lipid preferences determined by transmembrane helix, rather than by subunit. In heteroacid mixtures, cholesterol is enriched near the subunit interfaces, and even buried more deeply within the $\gamma/\alpha$ interface, with saturated acyl chains packed just outside cholesterol. The outermost M4 helices, however, were packed with PUFAs; PUFA chains also diffused throughout the TMD bundle. In homoacid membranes, saturated lipids and cholesterol were depleted among embedded lipids, but they maintained the same helix association observed in heteroacid membranes. PUFAs were especially enriched in domain-forming membranes, with high densities around all four transmembrane helices.

3.2 nAChR clustering in the presence and absence of lipid domains

Here, we repeatedly observed spontaneous formation of receptor dimers in the simulations containing multiple proteins (Figure 2). To investigate the role of lipid domain formation in driving nAChR oligomerization, we calculated the radial distribution function for the pairwise distances between centers of mass, as shown in Figure 6. For both heteroacid and homoacid membranes, we observed a peak at $r\sim 7.5$ nm, corresponding to dimerization. The differences in profiles between domain and non-domain forming membranes were primarily quantitative: most significantly, the peak corresponding to dimers was substantially higher for 2 and 3 proteins in homoacids than in heteroacids. For two proteins, the distribution for homoacids was shifted to the left, relative to the distribution for heteroacids, indicating that the peak for dimerization was at a shorter distance in domain-forming membranes. This difference is consistent with results indicating a role for domain formation in aggregation and clustering of nAChRs (Barrantes, 2007; Oshikawa et al., 2003; Pato et al., 2008). The simplest explanation of this difference is the higher effective concentration of proteins when domains are present: all nAChRs are corralled in a single liquid-disordered domain, with approximately half the area of the overall membrane.

3.3 Closest subunits across dimerizing proteins

In order to determine whether nAChR dimers were more likely to form with specific subunit interactions, we determined the closest interacting subunit pair for nAChR dimers within each frame. To ensure that dimers, rather than larger oligomers, were being analyzed, we excluded trimers and tetramers from the analysis and only considered two protein systems, with increased sampling. Figure  7 shows the amount of enrichment for each possible subunit pairing, relative to an expected random distribution. Results were very sensitive to the use of domain-forming compositions. In heteroacidic membranes, the $\alpha_{\delta}$ subunit formed a monomer-monomer interaction with the $\beta$ subunit, while the $\alpha_{\gamma}$ most favorably interacted with the $\delta$ subunit. In domain-forming membranes, the $\delta-\alpha_{\delta}-\gamma$ interface was two to five times as likely to pair with a $\alpha_{\delta}$ subunit compared to a random distribution, but substantially more dispersion was detected than in non-domain forming membranes. This difference suggests that nAChRs have preferred dimerization orientations, but membrane organization can reorient receptors. Although $\delta$ subunits are linked by a disulfide bond at the NMJ (Chang and Bock, 1977), our simulations suggest that without this link, $\delta$ subunits do not form the closest pair among dimerizing proteins.

4 Discussion

The nAChR is one of the most well-studied, fundamental pLGICs for understanding human cognition, memory, and muscle contraction (Gotti et al., 1997). As an integral membrane protein, the function and organization of nAChR is strongly dictated by its surrounding lipid environment. DHA is an $\omega$-3 polyunsaturated fatty acid abundant in synaptic membranes and the Torpedo electric organs, which are both native nAChR membranes. With its six double bonds, DHA is considered highly disordered and can induce domain formation in membranes (Wassall and Stillwell, 2008). We previously observed (Sharp et al., 2019) favorable interactions between a single nAChR and DHA-rich, cholesterol-poor domains using coarse-grained simulations. Here, we have extended this approach to investigate the role of lipid topology and domain-formation on boundary lipids of individual nAChRs. We have also conducted the first simulations containing multiple nAChRs, which has allowed us to observe spontaneous dimerization.

While DHA is implicated in human health and disease, (Lavandera et al., 2017) experimental studies considering its interactions with nAChR have exclusively considered its free-fatty acid form (FFA), (Antollini and Barrantes, 2016) rather than as an acyl chain component of a phospholipid. Application of $\omega-3$ FFAs causes a significant reduction in open times observed through single-channel recordings (Bouzat and Barrantes, 1993). Here, we observe a substantial effect of lipid topology on both embedded and annular lipids: DHA chains in homoacidic phospholipids are far more likely to be found as either annular or embedded boundary lipids. We previously observed only quantitative effects of swapping PE with PC on partitioning, but the headgroup does serve to anchor the lipid at the membrane/protein interface (Sharp et al., 2019). DHA in its FFA form (without a headgroup) may diffuse into an open nAChR pore, blocking the channel.

We find that, consistent with coarse-grained MD simulations using one nAChR (Sharp et al., 2019), multiple nAChRs continue to prefer the liquid-disordered phase containing long-chain $\omega$-3 fatty acids. While the number of nAChRs in the system did not affect the partitioning profile in these simulations, it did affect the composition of embedded lipids. Our results are consistent with each nAChR monomer in an oligomer blocking access to embedded lipid binding sites in the other monomers, suggesting an intriguing coupling between specific binding and membrane organization.

Interestingly, upon removing membrane organization, embedded lipids cluster around specific transmembrane helices in a five-fold symmetry around nAChR. This finding suggests that intrinsic lipid preferences are primarily helix dependent, rather than subunit dependent. In homoacid membranes, a lipid preference for DHA was observed across all transmembrane helices, with shells of PUFAs found even at the border of the liquid-ordered phase. Although saturated acyl chains and cholesterol were generally depleted around nAChR, the highest densities for both lipids were found around the M1 and M3 helices, as seen in heteroacid membranes.

In native membranes, nAChR dimers can be stabilized by a disulfide bond between $\delta$ subunits (Chang and Bock, 1977). An early controversy (Anholt et al., 1980; Rüchel et al., 1981; Zingsheim et al., 1982; Schindler et al., 1984) concerned whether the disulfide bond was necessary for dimer formation. There is no mechanism for covalent bonds between monomers in these coarse-grained simulations. All dimers were stabilized by non-covalent interactions, consistent with the results of (Rüchel et al., 1981; Schindler et al., 1984). We observed far more stable dimers in the homoacid mixtures than the heteroacid mixtures, which would be consistent with high sensitivity to experimental conditions. We did not observe a consistent $\delta-\delta$ preference for interfacing subunits in either homoacid or heteroacid mixtures. Schindler et al (1984) observed an apparent gain-of-function for single-channels within dimers relative to monomers, regardless of disulfide linking (Schindler et al., 1984), suggesting that lipid modulation of oligomerization could also provide a pathway for modulating single channel function.

Acknowledgements.

GB was supported by research grants NSF MCB1330728 and NIH P01GM55876-14A1. GB and LM were also supported through a grant from the Research Corporation for Scientific Advancement. This project was supported with computational resources from the National Science Foundation XSEDE program through allocation NSF-MCB110149, a local cluster funded by NSF-DBI1126052, the Rutgers University Office of Advanced Research Computing (OARC) and the Rutgers Discovery Informatics Institute (RDI2), which is supported by Rutgers and the State of New Jersey. We are grateful to Dr. Jérôme Hénin for his helpful suggestions throughout this study.

5 Compliance with Ethical Standards

Conflict of interest: The authors declare that they have no conflict of interest. This research was supported in part by the National Science Foundation, the National Institutes of Health, and the Research Corporation for Scientific Advancement.