Spontaneous Lipid Binding to the Nicotinic Acetylcholine Receptor in a Native Membrane

The nicotinic acetylcholine receptor (nAChR) and other pentameric ligand-gated ion channels (pLGICs) 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 only a few lipid species. Here we use coarse-grained molecular dynamics (MD) simulation to probe specific binding of lipids in a complex quasi-neuronal membrane. We ran a total of 50 microseconds 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 intersubunit sites and PUFAs select for convex M4 sites, while monounsaturated and saturated lipids are unenriched in the nAChR boundary. In order to characterize binding to specific sites, we present a novel approach for calculating a “density-threshold affinity” from continuous density distributions. We find that affinity for M4 weakens with chain rigidity, which suggests 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. Combined, these insights may reconcile several apparently contradictory experiments on the role of anionic phospholipids in modulating nAChR.


I. INTRODUCTION
The nicotinic acetylcholine receptor (nAChR) is a well studied excitatory pentameric ligand gated ion channel (pLGICs). 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 memory 1 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][27][28] . Functional studies using Xenopus oocytes [29][30][31][32][33][34][35][36] require lipid additives such as asolectin 20,[29][30][31][33][34][35][36] or lipids from synaptic membranes 37 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.
Functional experiments have focused on the role of anionic lipids and cholesterol as modulators of pLGICs [18][19][20][21][22][23][24][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.
As for most membrane proteins, it is experimentally challenging to capture the boundary lipid composition of pLGIC 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 TMD and the ECD. In the TMD, each subunit has four membrane helices (M1-M4) with the five subunits forming a "star" shape around a central pore ( Figure 1A). The M2 helix lines the pore, the M1 and M3 helices form a middle ring that includes the intersubunit cavities, and the M4 helices form the tip of the star. Structural methods have resolved potential cholesterol molecules 65,66 and phospholipids [67][68][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 intersubunit sites 65,72 . Unfortunately, fully atomistic simulations suffer from slow diffusion of lipids within the membrane, which prevents spontaneous lipid sorting by proteins over accessible simulation time scales.
Coarse-grained MD simulations use simplified molecular models that can reveal spontaneous lipid sorting, domain formation, and protein partitioning over simulation timescales 73-76 .
Coarse-grained MD simulations have been used previously to probe interactions of pLGICs with propofol 77 as well as spontaneous lipid binding in model membranes 70,78,79 . In previous work, we found that nAChR embedded in multiple domain-forming model membranes partitioned to the PUFA-rich liquid disordered domain 78 , 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 five-fold symmetric sites emerged: an intersubunit site and the M4 site ( Figure  1B). Cholesterol and saturated chains were enriched at the inter-subunit interfaces and n-3 PUFA acyl-chains were enriched around the M4 helices 79 . These results were consistent with binding to minimize packing defects: the rigid lipids could fill in the concave regions at the intersubunit 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.
With collaborators in the Cheng lab, we recently 70 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 1palmitoyl-2-oleoyl phosphatidylglycerol POPG in the intersubunit sites (inner leaflet); these sites contained basic amino acids that were also implicated through mutagenesis 70 . In nAChR the highdensity of basic amino acids are in the M4 site (inner leaflet) rather than the intersubunit 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 growing more feasible. In 2014, Ingólfsson et al 46 simulated an "average mammalian" membrane containing 63 lipids species, followed in 2017 by a coarse-grained neuronal membrane 41 . Multiple accessible and realistic membranes have been developed for comparison of protein-lipid interactions between model and quasi-native membranes 75,76,[80][81][82] . 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 like DHA vs n-6 PUFAs like α-linolenic acid) or like-charged head groups (PC vs PE, or phosphoserine (PS) vs phosphoinositol (PI)). For this work, we embed the neuromuscular nAChR 83 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 results and discussion of site selectivity of neutral lipids, followed by a reoriented discussion of the same data that is 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.2 84 topology and forcefield. nAChR coordinates were based on a cryo-EM structure of the αβ γδ muscle-type receptor in native torpedo membrane (PDB 2BG9 83 ). This is a medium resolution structure (4Å) and was further coarsegrained 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 structure in 2BG9 was 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 Ingólfsson et al 41 . The neuronal membrane from described by Ingólfsson contains phospholipids, sterols, diacylglycerol, and ceramide. Membranes presented in this paper only consider phospholipids and cholesterol, for a total of 36 unique lipid species, see Table SI 1. Coarse-grained membranes were built using the MARTINI script insane.py 85 , 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 40x40x35 nm 3 with ∼ 4, 500 − 5, 000 lipids and total ∼ 450, 000 beads.

B. Simulations
Molecular dynamics simulations run using the MARTINI 2.2 84 forcefield and GROMACS 86,87 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 1000000 steps, but energy minimization tended to concluded after ∼ 5000 − 10000 steps.
Volume and pressure equilibrations were run with isothermal-isochoric (NVT) and isothermalisobaric (NPT) ensembles respectively. NVT and NPT simulations used a time step of 15 fs and run for 0.3 ns using 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.
Molecular dynamics simulations were run using a time step of 20 fs for 5 µs for 10 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 pressure coupling constant set to 3.0 ps.
Secondary structures restraints with MARTINI recommendations were constructed by the martinize.py 84 script and imposed by GROMACS 86,87 . The nAChR conformation was preserved by harmonic bonds between backbone beads separated by less than 0.5 nm and calculated using the ElNeDyn algorithm 88 associated with MARTINI 84 with a coefficient of 900 kJ·mol −1 . Boundaries of the pseudo-symmetric intersubunit (black) and M4 (white) sites. The angular components are determined by the location of the M1 and M3 alpha-helices for either two adjacent subunits (intersubunit sites) or a single subunit (M4 sites), and are listed in Table S5. Circles correspond to the helices shown in panel A. (c) The distributions P site (n) (blue) and P bulk (n) (dashed red) represent the probability distributions for number of beads of a certain lipid species in the site or in an analogously-sized area of the bulk, respectively. The value n peak maximizes P bulk . P < (pink) is the area under P site to the right of n peak while P > (light blue) is the area under P site to the right of n peak .

C. Calculation of Polar Density Distributions
As in our previous work 70,78,79 , 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 r i = 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 ∆θ = π 15 radians are the bin widths in the radial and angular direction respectively, and n B (r i , θ j ) is the time-averaged number of beads of lipid species B found within the bin centered around radius r i and polar angle θ j . In order to determine enrichment or depletion, the normalized densityρ B (r i , θ j ) is calculated by dividing by the approximate expected density of beads of lipid type B in a random mixture, where s B is the number of beads in one lipid of species B, N L is the total number of lipids in the system, and L 2 is the average projected box area: where the expected density is derived at the first frame of the simulation. Python software for these calculations are under active development and are located at 89 .
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 5 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 to 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 P site (n) of finding n beads within the site and the probability P bulk (n) of finding n beads within a region of equivalent area in the bulk, respectively. For a lipid that has no affinity for this binding site, we expect P site (n) = P bulk (n), while P site (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 n peak that corresponds to the peak of the density distribution in the bulk. As illustrated in Figure 1 C, we then integrate P site on both the left and right side of the threshold n peak to yield P < and P > respectively: Note that this step breaks the distribution into two macrostates on either side of the threshold, allowing clear analogy with 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 "densitythreshold affinity". In the special case of binary occupancy, where K D 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, so P bulk (n) = 1 for n = 0 and vanishes for all n > 0, so n peak = 0. Consequently, for this special case, P

E. Binding Site Definition and Occupancy Calculations
We consider two classes of site: intersubunit sites and M4 sites. Each pLGIC has ten of each site (five in the outer leaflet and five in the lower leaflet) for a total of twenty sites ( Figure 1B). The boundaries for each site were drawn to correspond to the localized binding hot spots observed for heteroacidic membranes 79 , and are non-overlapping. Inter-subunit sites include bins with angular components between the M1 and M3 alpha-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 alpha-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 P site (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 ∆θ = π 25 . 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.

F. Calculation of Accessible Area
Calculation of P bulk 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 A tot and the area A ex excluded by the protein: A = A tot − A ex A tot is straightforward to calculate by summing over the areas of the bins i within the site boundaries: A tot = ∑ i r i ∆r i ∆θ i .
Calculating A ex 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, P site (n) = P bulk (n) as long as P bulk (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 ∼ 370ns and determining the value of A for each site such that P site (n) and P bulk (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 Figure 2A, with representative poses shown in Figure 2B. 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 around the protein.
This enrichment is shown in Figure 3A  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: Figure 4 provides the "lipid's perspective" and is organized to identify the preferred site for a given lipid type (the lipids' "site selectivity"), while  Figure 5 has an alternate representation of the same data. Figure 5 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 Figure 4A, 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 Figure 4A, saturated lipids have similar density-threshold affinities across all sites, which is consistent with the generally flat distribution observed in Figure 3. Yet saturated lipids do yield a slightly stronger affinity for intersubunit 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 intersubunit 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 sites. As shown in Figure 4A, neutral PUFAs have significantly stronger affinities for M4 sites than for innersubunit sites in the same leaflet. Such selectivity is consistent with the PUFA enrichment density in Figure 3A, 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 ( Figure 3C). 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 intersubunit sites. Conversely, cholesterol has significantly stronger affinities for innersubunit sites than for M4 sites, which is also consistent with the enrichment density in Figure 3A and our expectations from model membranes ( Figure 3C). 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 intersubunit and M4 sites
We now switch perspectives to considering which neutral lipids are most favorable for particular sites. As shown in Figure 5 A Figure   S1.
PUFA chains yield affinities for the intersubunit site that are approximately >0.5 kcal/mol stronger than saturated lipid affinities ( Figure 5 A and B), which was unexpected based on results from model membranes but is consistent with the corresponding enrichment density in Figure   3A. More generally, neutral phospholipid affinities for intersubunit sites obey the following trend, from strongest to weakest: n-3 > n-6 > saturated > monounsaturated. Thus, even though PUFA chains prefer M4 sites to intersubunit sites, and saturated chains prefer intersubunit sites to M4 sites, PUFAs have a stronger affinity for either site type than do saturated lipids.
For intersubunit 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.
As shown in Figure 5C and 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. [90][91][92][93][94][95] Surprisingly, cholesterol had a stronger affinity for M4 sites than any acyl chains other than n-3 PUFAs. Cholesterol is rigid, small, and has asymmetric sides (rough and smooth) which potentially allows it to embed between alpha-helices and compete with n-3 PUFAs for binding.
Any cholesterol bound within the grooves of the subunit interface (as hypothesized based on atomistic simulations 71 and observed in β subunits of nAChR (using coarse-grained simulations 78 ), will also get counted within the M4 site. 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 subunit ( Figure ??) As shown in Figure   3b, anionic lipids are generally enriched around the M3/M4 helices for the α γ , γ, δ , and β subunits.
Anionic lipids are enriched at intersubunit sites and around M4 sites for all subunits but the α subunits. Non-α nAChR subunits have basic amino acids closer to M4 alpha-helices, as shown in Figure SI 2a. We incorporate data from all five pseudo-symmetric sites to obtain the densitythreshold affinities reported in Figure 4B, which suggest that anionic lipids have significantly stronger affinities for M4 sites on average. The average anionic affinity difference between intersubunit and M4 sites is ∼ −1.0 kcal/mol, as shown in Figures 4, 5d 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 Figure 5C In summary, we observe that binding sites have clear preferences for particular head group charge and acyl-chain saturation, suggesting nAChR lipid occupancy to be driven in two steps, a "coarse-sorting" by head groups, and then "fine-sorting" by acyl-chains. A neutral lipid will occupy nAChR's boundary region but acyl chains dictate where specific lipids occupy nAChR.
Anionic lipids diverge from this pattern at the inner M4 site which has the strongest affinity for anionic lipids independent of saturation.

IV. CONCLUSIONS
Using coarse-grained simulations of the nAChR within a quasi-neuronal membrane containing over thirty 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 (intersubunit site) and convex (M4 site) portions of the star-shaped nAChR and were initially observed as "hot spots" in our previous simulations 70,79 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 lab has developed promising alchemical approaches 96 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 straighforward to extend to any other transmembrane proteins with detectable regions of density enrichment.
Our results are summarized graphically in Figure 6. Based on our results from model mem-  71 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 symmetriclyarranged, while in the heteromeric nAChR they vary by subunit (Figure ??a), 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, neutral lipids with smaller head groups (PE) have stronger affinity than the larger PC headgroup. It is unclear why PE is more favorable than other neutral lipids at this time, though this is consistent with previous work 70,78 , and the most straightforward explanation is that the smaller headgroup introduces fewer clashes with the protein TMD.
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 forcefield 84 carries a charge of -2, and PIP1, PIP2, and PIP3 have charges of -3,-5, and -7. 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.
The present results highlight the utility of model membranes for developing hypotheses of specific lipid-protein interactions, and the need to test those hypothesis within more complex native membranes. The present results could 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 intersubunit 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 lipid competition, which is one of the primarily challenges of interpreting experiments in complex membranes.

ACKNOWLEDGMENTS
GB and LS were supported by the Busch Biomedical Foundation. This project was supported by generous allocation through the Rutgers University Office of Advanced Research Computing (OARC), which is supported by Rutgers University and the state of New Jersey. We are grateful to Dr. Jérôme Hénin for helpful input and suggestions.

DATA AVAILABILITY STATEMENT
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.