Abbreviations

---

1. Introduction

Elucidating the interactions of proteins with their native membrane environment is central to our understanding of the function of cell membrane proteins and receptors ^[1,2]. Atomistic force fields enable detailed simulations of membrane dynamics ^[3,4,5,6], but the lateral diffusion of lipids in the bilayer is slow to converge and requires large system sizes ^[7]. Coarse-grained (CG) model representations group atoms into effective interaction centers to strike a balance between chemical detail and computational efficiency, expanding the range of problems that can be explored ^[8]. Coarse graining lipid bilayers has become a powerful tool for studying the dynamics of complex multicomponent membrane systems ^[9]. CG models of polypeptides, initially developed to study protein folding ^[10], can now be applied to conformational transitions and protein functional mechanisms ^{[11,12,13,14]}. To understand how the local membrane environment modulates protein function, care must be taken to incorporate a protein's conformational flexibility in the context of the lipid bilayer environment ^{[15,16,17,18,19,20]}. CG simulations have identified direct and indirect effects on protein function arising from the presence of cholesterol as well as anionic ^[21,22,23], zwitterionic and double phospholipids in the bilayer ^[20,24,25].

MARTINI is a widely-used CG force field rooted in a thermodynamics-based parameterization approach ^[26]. The resulting protein and lipid models have been characterized rather extensively and parameter sets are now available for dozens of membrane lipid types ^{[27,28,29,30,31,32,33,34,35,36]}. While simulations with the model have yielded numerous insights into multicomponent systems of biological interest ^[9,13,34,37], there are limitations inherent in the model ^[29,38]. The nonbond interaction potentials rely on the Lennard-Jones 12-6 functional form, which for CG representations has been shown to have a core repulsion that is steeper than the corresponding atomistic molecular dynamics (MD) distributions ^[39,40]. Another potential issue is that MARTINI solvent particles, each of which represents four water molecules, have a tendency to freeze at room temperature ^[38].

Explorations of protein conformational dynamics will be limited if the CG pseudodihedral potential energy function fixes the secondary structure assigned in the starting conformation ^[37]. Flexible bonded and nonbonded potentials are useful in protein folding studies but structural collapse is a potential problem when implementing them in CG force fields not intended for folding ^[40,41]. Other ongoing developments in the field include incorporating polarizable CG water beads that are still consistent with the protein parameters and capture the context dependence of interactions ^{[41,42,43,44,45]}. Martini version 2.2P incorporates polarized polar sidechains in order to obtain the expected attractive behavior at the membrane interface ^[46]. Dry Martini avoids the issues associated with large CG solvent particles by removing the explicit waters but retaining all bilayer lipids ^[33,47,48]. To better balance the treatment of protein versus solvent water interactions, Han et al. developed the PACE united-atom protein force field for simulations in the presence of CG water ^[49,50]. Unlike the Cα-only model employed for the backbone in many CG models, incorporating explicit backbone amide atoms is ideal for studies of protein folding that depend on backbone hydrogen bonding. A promising new CG approach that allows for secondary structure formation is to introduce a polarizable dipole into a Cα-only backbone bead ^[51].

The present work explores the utility of the CgProt force field for simulations of peptides and proteins in an explicit bilayer. As opposed to parameterizing particles individually using thermodynamic solvation data, CgProt is rooted in the multiscale coarse graining (MS-CG) method ^[40,52]. MS-CG takes a large set of reference configurations in one or more atomistic MD simulations and uses a variational force-matching scheme to simultaneously develop a self-consistent set of CG interaction potentials. Tabulated potentials can be constructed for the nonbond interactions with no a priori assumptions of their functional form ^[53]. Nonbond tables with 0.005 nm or smaller bins can then be read in by the Gromacs simulation package for efficient CG-MD simulations ^[54]. A technical limitation is that GPU-acceleration is not yet supported (as of Gromacs 2016.3), because the group cutoff scheme is needed for tabulated potentials. To avoid the pitfalls associated with bulky CG water particles, the water degrees of freedom are integrated out with their effects on solvation captured implicitly in the nonbond potentials ^[47]. Note that because of the additive nature of many-body interactions, each pairwise MS-CG potential in the set is not equivalent to a potential of mean force (PMF) calculation obtained from on an individual particle pair ^[53].

The CgProt protein model is derived from folding simulations in water ^[40], while the lipid-protein interactions were developed from MS-CG calculations on unfolded peptide-phospholipid ensembles ^[20]. The lipid parameter set was recently improved upon to include lipids abundant in eukaryotic and prokaryotic membranes: POPC, POPE, POPG, cardiolipin and cholesterol ^[55]. To obtain better agreement with atomistic simulations, a more repulsive potential was incorporated for phosphatidyl ethanolamine (PE) head groups, along with a harmonic torsion to capture the configurational influence of a cis double bond on the acyl chain.

In the present work, the mechanism of insertion was examined for peptide helixes WALP27, LS3 and L24 in a DOPC bilayer. Multiple, independent unbiased CgProt simulations were performed starting from either inserted or desorbed peptide conformational states. Next, CgProt was used to demonstrate the conformational transition of the 1,188-residue, 12-transmembrane helix ABC transporter P-glycoprotein (P-gp). The conformational change involved a hinge bending motion and reorientation between the protein's two pseudosymmetric halves and did not require changes in secondary structure. To examine the potential role of lipid-protein interactions, molecular dynamics were performed in a 17 nm patch of model eukaryotic membrane containing a symmetric 2:2:1 ratio of POPC: POPE: cholesterol.

---

2. Materials and Methods

In CgProt 2.3 the polypeptide is represented as a Cα chain with up to 4 interaction sites defined for each sidechain (Figure 1). To match fluctuations observed in atomistic MD, harmonic Cα-Cα bonds are assigned an effective force constant of 80,000 kJ mol^–1 nm^–2 with 3.86 Å equilibrium length. By comparison, variants of the Martini model limit all bond constants to a maximum of 7500 kJ mol^–1 nm^–2 to allow for a larger integration time step ^{[11,56,57,58]}. CgProt employs double-well bonded potentials for 3-body backbone Cα angles and 4-body pseudotorsions that allow for both α-helix and β-sheet secondary structure ^[11,40,59]. To maintain amino acid chirality ^[60,61], an improper torsion term is applied to residues with a Cβ site in the sidechain ^[11]. All 1, 2 and 1, 3 bonded interaction pairs are excluded from nonbonded interactions, as are 1, 4 backbone Cα pairs. Flexible angle terms for CG sidechain sites were developed as 4th order polynomials via Boltzmann inversion of flexible peptide simulations ^[40]. Only the Trp and Tyr sidechains have dihedral terms to maintain ring planarity. Harmonic sidechain pseudobond potentials are representative of their corresponding atomistic fluctuations in solution with a maximum force constant of 100,000 kJ mol^–1 nm^–2 to allow for up to an 8 fs integration time step ^[11]. The liquid bilayer phase is maintained at kT = 2.24 kJ/mol in Gromacs 4.6.7 using stochastic Langevin dynamics with a damping friction coefficient of γ = 0.5 ps^{–1 [54,55]}. The neighbor list is updated every step for the tabulated nonbond potentials (1.2 nm cutoff distances), which combine van der Waals and electrostatic contributions into a single PMF. As γ = 50 ps^–1 in real water ^[62], the CG bilayer and implicit solvent represent a significant computational speedup. Lipids in CgProt exhibit lateral diffusion coefficients ≥ 250 μm² s^–1, a 30-fold enhancement relative to experimental values for the liquid bilayer phase ^[7,20].

---

2.1. Peptide insertion simulations

The mechanism of membrane binding and insertion of synthetic peptides is explored in a DOPC bilayer ^{[30,31,63,64,65,66]}. The sequence GWW(LA)₁₀LWWA (WALP27) of the WALP n series of single-pass transmembrane peptides was simulated since its 31.5 Å hydrophobic length is the closest match for the hydrophobic thickness of the DOPC bilayer: 29 Å as measured by the mean distance between acyl ester groups ^[67]. The highly charged transmembrane peptide K₂L₂₄K₂ (L24) was also assessed for its bilayer insertion stability ^[68]. Third, the designed amphipathic peptide (LSSLLSL)₃, referred to as LS3, was assessed for its binding behavior at the membrane-water interface ^[69]. Peptide N-and C-terminal alpha carbons were assigned the default positive and negative site types, respectively, as there is no Coulomb potential in CgProt.

While previous simulations relied an elastic network model of native contacts to stabilize the backbone in tertiary protein domains ^[20], distance restraints were not sufficient to stabilize a single ideal α-helix. In the case of LS3, polar residues did not remain aligned on the same side of the helix. Backbone alpha carbons were therefore restrained using GROMOS-96 angles (θ₀ = 89°, k_θ = 2000 kJ mol^–1) and harmonic dihedrals (φ₀ = 50°, k_φ = 2000 kJ mol^–1 rad^–2) in all three peptides. For MD analysis, translational motion of the bilayer center of mass in the z-dimension was removed from each trajectory by root-mean-square alignment of the ester sites with the starting configuration. The angle of the peptide insertion with respect to the bilayer normal was calculated using the vector spanning alpha carbons 5 and 16 for LS3, 6 and 24 for WALP27, and 8 and 26 for L24. Peptides were simulated in multiple independent 500 ns simulation runs with 1,148 DOPC molecules in a 20 nm square box, each requiring 80 hours on 48 CPUs.

---

2.2. P-gp membrane simulation

The refined crystal structure (PDB ID: 4M1M) for mouse P-gp was simulated to investigate the conformational dynamics of the ABC transporter in a eukaryotic model membrane ^[70]. The missing flexible linker in the structure between its N-and C-terminal halves was omitted; these two terminal residues employed the Cα nonbond site type (uncharged). The seven cis peptide bonds in P-gp were assigned a 3.1 Å harmonic Cα-Cα bond. While complex algorithms have been developed for inserting protein structures into a lipid bilayer ^[71], the following procedure was used to expedite creating a pore around the protein. System preparation began with deleting the L24 peptide coordinates from a 17.1 nm square bilayer patch previously equilibrated with 472 palmitoyl-oleoyl phosphatidyl choline (POPC), 446 palmitoyl-oleoyl phosphatidyl ethanolamine (POPE), and 230 cholesterol lipids (50.7 Å area/lipid). A 1 ns simulation was performed with applied surface tension (140 mN/m) to expand the pore left behind in the bilayer using Berendsen pressure coupling with a 2 ps coupling constant every 10 steps. The principal axis of P-gp was aligned normal to the bilayer for insertion into the pore. Lipids lining the pore were observed to collapse onto the protein surface within 750 ps of simulation at constant zero tension ^[72]. To maintain headgroup positions during opening and closing of the pore, phosphate sites were restrained in the z-dimension (k_z = 1000 kJ mol^–1 nm^–2).

The resulting P-gp system was first equilibrated for 50 ns with the protein backbone fixed. A production run of 0.5 μs was then performed with free translational motion of the protein and a stabilizing homogeneous elastic network applied for the backbone starting crystal structure. A total of 11,365 virtual harmonic bonds (k = 150 kJ mol^–1 nm^–2) were applied to Cα carbon pairs separated by two or more bonds that fell within a 10.5 Å cutoff distance in the crystal structure. Choosing this weak spring constant but large cutoff was previously shown to stabilize helical bundle and beta sheet domains while allowing the closing transition between two symmetric transmembrane (TM) domains in P-gp homolog, MsbA, via a flexible hinge region ^[20]. CgProt lipid bonded and nonbonded parameters were employed as reported previously ^[55]. The cholesterol bonded terms were adapted from Daily et al. with a harmonic force constant of 100,000 kJ mol^–1 nm^{–2 [73]}.

---

3. Results and Discussion

---

3.1. Peptide insertion simulations

Independent simulations were performed for the designed helical peptides LS3, WALP27, and L24. The first batch of five simulations started with each peptide displaced 10 nm into the aqueous phase, while the second batch of five simulations for LS3 and WALP/L24 started from interfacial and transmembrane conformations, respectively. The mean length of the simulated LS3 peptide from the Cα of residue 1 to that of residue 21 was 29.5 ± 0.3 Å. This is in good agreement with the 30 Å length observed in atomistic MD ^[74], confirming that an ideal α-helix was maintained by the angle and dihedral restraints. The LS3 interfacial conformation exhibited an average tilt angle of 74.0° relative to the bilayer normal, consistent with the N-terminus penetrating further into the bilayer than the C-terminus (Figure 2). The helix center exhibits a mean displacement of 11.6 Å along the bilayer normal (Table 1), placing it near the carbonyl region of the bilayer centered at 1.5 nm. LS3 rotates about the helical axis so that polar sidechains point away from the nonpolar bilayer interior.

For the LS3 simulations starting from the aqueous phase, the C-terminus adsorbs to the bilayer surface as early as 1.5 ns (Figure 3). By 3 ns, the N-terminus also adsorbs onto the bilayer. The resulting interfacial state resembles the equilibrated conformation in simulations starting from the interfacial state. In all five desorbed simulations and in three out of five interfacial simulations, the N-terminus eventually crosses over to the other monolayer, with LS3 remaining in the transmembrane orientation for the remainder of the simulations (Table 2). Simulations from the interfacial state took an average of 175 ns to adopt the transmembrane configuration, while desorbed simulations took only 81 ns. These results indicate that the interfacial conformation is a metastable state. It is noteworthy that in both the interfacial and transmembrane conformations the positive N-terminus resides at a lower depth in the bilayer than the negative C-terminus (Table 1), indicating preferential membrane interactions with the positive charge. Observing both interfacial and transmembrane orientations is consistent with previous LS3 simulations and the known ability of LS3 to oligomerize into transmembrane pores at high concentrations ^[66,74,69].

The fact that the peptide termini are charged helps explain why the transmembrane orientation of LS3 remains anchored throughout the second simulation. A PMF was previously calculated for the insertion of uncharged LS3 peptide in DPPC using the Bond force field ^[66]. The simulation with neutral termini obtained a 4:1 ratio of interfacial: transmembrane helix. The interfacial conformation had an 80–90° tilt angle and 15–18 Å helix displacement, while the transmembrane conformer had a tilt of 5–20° and 0–2 Å helix displacement. These states correspond well with the two conformers predicted by CgProt, though the transmembrane state is more stabilized in the present model.

In contrast to LS3, only the transmembrane conformation of WALP27 is stable in the desorbed and inserted simulations. Starting from the aqueous phase, the C-terminus first docks onto the bilayer surface as early as t = 3 ns (Figure 4). The N-terminus follows suit as early as 10 ns, but the interfacial orientation with the peptide aligned with the head group bilayer region is not stable. The N-terminus continually decreases its distance to the bilayer center. As early as 12 ns, Trp 2 forms a bridge to glycerol and ester lipid sites in the other monolayer. The N-terminus then fully crosses over to the other head group region in an average of 22 ns. The resulting stable conformation is indistinguishable from the transmembrane orientation observed in simulations started from the fully inserted state.

Transmembrane WALP maintains a 22° tilt angle with its helix center < 1 Å from the bilayer center. The C-terminus is found in the head group region of the bilayer near the positions of the phosphates (Table 1), while the N-terminus sits lower in the ester region that spans 1 to 2 nm from the bilayer center, again indicating enhanced permeability of the positive N-terminus. The observed degree of tilt is in agreement with the 24° ± 5° angle found by a recent fluorescence spectroscopy study of WALP23, which was intended to address discrepancies between atomistic MD simulations and previous NMR experiments ^[75,76]. In comparison to the transmembrane stability in CgProt, the insertion of neutral WALP in the transmembrane orientation has been observed with 80%, 100%, and 100% frequency using the MARTINI, PLUM, and GROMOS force fields, respectively ^[31,66].

A membrane insertion event was not observed during the desorbed simulations of the highly charged L24 peptide. Like LS3 and WALP, the C-terminus of L24 is first to adsorb to the bilayer surface. The N-terminus soon associates with the head group region of the bilayer and adopts a lower position than the C-terminus (Table 1). The interfacial conformation of L24 was then stable throughout the remainder of the simulations, exhibiting a mean tilt angle of 82° and z-displacement of 14 Å. During the simulations starting from the transmembrane state, L24 remains anchored in a highly stable transmembrane orientation (Figure 2). The transmembrane helix exhibits a mean tilt angle of 34° and z-displacement of < 1 Å. The relative positions of the N-and C-termini in transmembrane L24 are similar to the values observed for WALP, despite the two charged lysines at either end. Despite the leucine stretch of L24 being more hydrophobic than the Leu-Ala repeat in WALP, the terminal lysines apparently provide a considerable kinetic barrier to bilayer permeation. Taken together, these simulations reveal the utility and limitations of unbiased CG MD for investigating the membrane insertion of small peptides.

---

3.2. P-glycoprotein conformational transition

The 4M1M crystal structure (Figure 5) adopts an inward-facing conformation in which the two cytosolic nucleotide-binding domains (NBDs) belonging to its pseudo-symmetric halves are not contacting ^[70]. Binding of the two NBDs is required for hydrolysis of two ATP molecules at symmetrically opposed sites and subsequent substrate translocation as driven by the outward-facing conformational transition of the TM domains. In the course of the elastic network simulation, the 12 TM helices lining the accessible substrate binding chamber did not significantly change their relative orientations (3.1 Å final Cα RMSD), but motion about each TM-NBD hinge allowed the NBDs to bind each other. There was negligible intra-domain motion in each NBD (1.6 ± 0.1 Å mean RMSD).

A measure of structural openness between the NBDs is taken as the Cα distance between residues N607 and T1252 in P-gp (T561 in dimeric homolog MsbA). A range of distances is observed in crystallography and electron microscopy in the absence of substrate, suggesting flexibility in the inward-facing state ^[77]. In the CgProt simulation, the NBD-NBD distance rapidly decreases from 52 Å in the starting crystal structure to a final contacting distance of 30 Å. Analogously, MsbA has been crystallized in an inward-facing state in which the NBDs are partially contacting and an outward-facing ATP-bound state in which the NBD-NBD interface is fully formed. The NBDs are separated in these "closed-apo" and "nucleotide-bound" states by 37 and 35 Å, respectively ^[78]. The final conformation observed in the CG simulation resembles a pre-hydrolysis intermediate in which the NBD-NBD interface is fully formed and competent for ATP binding, but the opening between the TM domains still faces inward.

A homology model for P-gp was recently constructed based on the outward-facing crystal structure of Sav1866 and subjected to refinement in atomistic MD simulations ^[79]. One simulation in which the two MgATPs remained tightly bound resulted in an NBD separation of 37 Å (Table 3). In other simulation runs one or both of the ATPs became loosely bound, increasing the NBD separation to 39 and 41 Å, respectively. Runs performed in the absence of ATP had an NBD distance of up to 46 Å. The equilibrated conformation from the CgProt simulation resembles the "loose-tight occluded" conformation generated by atomistic MD. Each ATP binding site is formed by a Walker A motif on one NBD and an LSGGQ motif on the other NBD. Pan and Aller identified a tight ATP(O_γ) binding interaction by a 9.4 Å Cα distance between S528 and K1072 (site 2) and a loose ATP(O_α) site 1 interaction characterized by a 11-12 Å separation between the symmetrically equivalent residues K429 and S1173. The mean Walker A-LSGGQ distances in the CG simulation were 9.4 ± 0.6 Å and 17.6 ± 1.7 Å, respectively, with an even larger opening observed for site 1 (Figure 5). The equilibrated CG trajectory exhibits a mean Cα RMSD from the loose-tight structure of 5.3 ± 0.3 Å for the NBD-NBD interface after alignment in its own reference frame.

Several residues have been shown to participate in cysteine crosslinking across the NBDs ^[80,81], but this requires an alternate alignment between the NBDs. Contacts between the identified residues are only formed in the closed-apo crystal structure for MsbA, in which the NBDs are partly contacting but out of register ^[78]. Ward et al. posited that the closed-apo state is stabilized by bound substrate and precedes ATP binding. This would explain why ATPase activity is dramatically enhanced in the presence of substrate. The CgProt simulation adopts the fully formed NBD-NBD interface with large distances between the cross-link residues (Table 3). This is notable because the backbone potential incorporated an elastic network encoding only the open conformation. The network of contacts in the open structure thus appears predisposed for formation of the nucleotide-bound state but not the closed-apo state. In order to sample the closed-apo intermediate, an elastic network could be developed from the MsbA crystal structure based on a sequence or structural alignment with P-gp. To model transitions between the conformations, two structural states can be encoded simultaneously by the implementation of inverted Gaussian potentials in place of the unbreakable harmonic springs ^[11]. Changes in secondary structure have yet to be observed in elastic network simulations, but may be possible with attractive but dissociable Gaussian potentials.

---

3.3. Annular lipids

The P-gp simulation was performed in a physiologically representative 2 POPC: 2 POPE: cholesterol mixture to examine the role of interactions with the membrane environment. Previous simulations of MsbA in a model DOPC: DOPE bilayer revealed a solvation shell of non-specifically bound lipids bound in the grooves between TM helices ^[20]. With association energies of at most 1 kcal/mol and residence times below 100 ns, lipid tails were observe to dissociate from the protein and then be replaced by other lipids that adopted similar conformations on the surface of the protein, consistent with measurements on spin-labeled phospholipids ^[82]. To examine whether lipids preferentially associate with P-gp, the time-averaged lipid density was determined in the reference frame of the freely diffusing protein.

While POPC, POPE and cholesterol molecules do not segregate within the bulk membrane environment ^[55], sequestration was observed in the annulus of lipids that remained around P-gp (Figure 6) ^[84]. The cholesterol and phosphatidyl choline density is enriched relative to the bulk in the immediate vicinity of the protein, while no phosphatidyl ethanolamine molecules are found in the first solvation shell. A few phospholipids were also observed to diffuse into the large transmembrane chamber between the N-and C-terminal TM domains. The passage of amphiphilic lipids through the opening between the N-and C-terminal domains could potentially either mediate or compete with the uptake of xenobiotic substrates or drugs, which have been shown to have multiple distinct binding sites within the transmembrane chamber ^[85]. Interestingly, the smaller cholesterol molecules were not observed to enter the pore.

The lipid densities indicate that POPC and cholesterol have a slightly enhanced affinity for the exterior protein surface. Additional simulations would be needed to test whether there are any synergistic effects between the two lipids. An annular ring of loosely associated cholesterols could facilitate strong specific binding events as required for function ^[86] by enhancing the local sterol concentration. The preference of phosphatidyl choline for the protein was not found to be driven by direct interactions with charged sidechains. P-gp exposes predominantly basic residues to the head group region of the inner leaflet, and choline and ethanolamine experience a similar repulsive term for positive groups in the model. The bare ethanolamine does experience a stronger attraction for acidic sidechains in the model compared to choline's methylated amine, but POPE was not found in the vicinity of the several acidic sidechains exposed to the outer leaflet. Clusters of phospholipid density in the vicinity of MsbA were previously found where the protein radius decreases near the upper portion of the inward-facing conformation ^[20]. We conjecture that irregularly shaped proteins create voids in the bilayer that are better filled by the bulkier PC head groups, which have weaker interactions with phosphate than PE ^[87].

---

4. Conclusion

CgProt enables the rapid and favorable membrane insertion of small peptides that bear minimal charge. All observed insertion events involved the N-terminus crossing the bilayer. Multiple lines of evidence point to the attractive nature of positive-charged peptide moieties for the interface of zwitterionic membranes ^[88,89,90]. PMFs for sidechain insertion in DOPC, for example, suggest that protonated arginine and lysine are strongly attracted to the carbonyl region of the bilayer and experience a lower permeation barrier, compared to the strong repulsion experienced by deprotonated glutamate and aspartate ^[91]. The attraction may in part be due to water defects that have been identified in atomistic simulations, which involve the snorkeling of polar head groups to accommodate the buried charge. A similar interaction is observed in CgProt to mediate the crossing of the N-terminus of WALP between monolayers. Our observations that the free N-terminus sits deeper in the membrane than the C-terminus for three different peptides is supported by double electron-electron resonance measurements of WALP23's position in a 7:3 DPPC:PG bilayer ^[92].

Elastic network simulations of the P-gp open to closed conformational transition in conjunction with the CgProt force field lend support to the loose-tight structural model of ATP binding recently generated by homology modeling and atomistic MD simulation. The structural model of Pan and Aller contains an outward-facing arrangement of the TM domains, which though extensively characterized for bacterial protein MsbA, has remained elusive in experiments of the mammalian transporter P-gp. Future conformational studies with CgProt will explore the use of dissociable Gaussian potentials for modeling more subtle changes in supersecondary or tertiary structure.

In addition to driving the insertion of transmembrane helices, lipid-protein interactions are implicated in regulating the function of transmembrane proteins. Distinct from high affinity lipid binding at specific protein sites ^[93], molecular dynamics with CgProt reveal that the large hydrophobic TM region of P-gp and the conformational changes therein alter the local structure and diffusion of the surrounding annular lipids ^[20,84]. Altogether, these molecular dynamics simulations of LS3, WALP, L24, and P-gp illustrate the potentially important structural and functional roles of lipid-protein interactions and demonstrate the necessity of conducting dynamics simulations of proteins in the context of their native membrane environment.

---

Acknowledgments

This work was supported by National Science Foundation grant MCB-1516826 (to R.D.H.). Cody Black was supported by an American Foundation for Pharmaceutical Education Gateway to Research Award.

---

Conflict of Interest

The authors declare no conflicts of interest in this paper.

---