Subscriber access provided by Stockholm University Library
Biomolecular Systems
Analysis of Lipid Order States and Domains in Lipid Bilayer Simulations Soohyung Park, and Wonpil Im J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00828 • Publication Date (Web): 23 Nov 2018 Downloaded from http://pubs.acs.org on November 27, 2018
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Analysis of Lipid Order States and Domains in Lipid Bilayer Simulations
Soohyung Park* and Wonpil Im* †Departments
of Biological Sciences and Bioengineering Program, Lehigh University, Bethlehem, Pennsylvanian
Corresponding Authors *
[email protected];
[email protected] ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
ABSTRACT We propose a general procedure to analyze lipid order states and domains in lipid bilayer simulations using surface areas and hydrophobic thicknesses of lipids. In our approach, the observable order states of individual lipids are inferred by a hidden Markov model analysis of their time series and by considering the deformation of a lipid at different packing environments. The assigned lipid order states are mapped onto the Voronoi tessellation of lipids, from which the ordered/disordered lipids are robustly clustered by the Getis-Ord local spatial autocorrelation statistics. The usefulness of this method is illustrated by its application to the quinary mixed bilayers consisting of cholesterol (Chol), 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), 1,2-dimyristoyl-sn-glycero-3phosphoethanolamin (DMPE), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE), where any type of phospholipids does not show strong preference over the other types to be enriched in lipid domains. The independent order state analysis for each lipid type allows straightforward applications of our method to arbitrarily complex bilayer simulations. KEYWORDS Lipid domain; Order parameter; Hidden Markov model; Spatial autocorrelation statistics
ACS Paragon Plus Environment
Page 2 of 22
Page 3 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Introduction Lipid rafts are nano- or micro-sized dynamic membrane domains, and they are involved in many important biological processes such as membrane trafficking, signaling, protein sequestration, and so on.1,2 Cholesterol, glycosphingolipids (and other saturated lipids), and specific membrane proteins are known to be enriched in lipid rafts.1,2 While the existence of these domains has been supported by studies using various experimental techniques3–11 it is still challenging to experimentally monitor their dynamics (size, composition, and lifetime) due to various limitations (detection limits such as resolution in optical microscopy and the intrinsic limitations in the experiments such as ambiguity in lateral motion of the species8 and their compositions9 in vesicle experiments and slower dynamics due to the interactions with the support in supported bilayers10,11). Molecular dynamics (MD) simulation can be a useful approach that complements experiments and provides insight into the properties and dynamics of domains by monitoring and analyzing membrane structures at the molecular level. So far, in most computational studies, the lipid domains were analyzed for model membranes composed of a few types of lipids (binary or ternary mixtures with or without cholesterol) with clearly different properties (such as tail saturation and melting temperature).12–17 In these model membranes, the domains were assigned based on an observable such as the surface area,12 order parameters of lipid acyl chains13–15 and their lateral packing,13 or the local lipid compositions,16,17 which showed clear contrast between the low and high ordered regions. While these observables have been effective, there are potential issues that could complicate the domain analysis of complex multi-component bilayers as described below. It has been reported that the liquid-ordered phase (Lo) of ternary mixed bilayers with cholesterol (Chol), low melting, and high melting lipids has substructures of hexagonally packed lipid arrays bordered by Chol.12,17 In this case, order parameters from local chain packing would result in mis-assigned lipid domains due to low order parameter values around Chols, as Chol is known to disrupt surrounding hexagonal lipid packing in Lo. In another case of mixed bilayers consisting of Chol and lipids with subtly different physical properties (e.g., phospholipids of the same acyl chains with choline and ethanolamine head groups), a specific type of lipid would not be strongly preferred in domains over the other types.18,19 In such cases, neither the chain packing order parameter nor the local composition approach (looking for aggregation of a certain type of lipids) can assign lipid domains correctly. The chain orientation order parameters could handle this situation, but it needs careful verification whether a single criterion can be applied to different lipid types. Thus, it could be an issue whether the observables for order parameters is properly chosen for the bilayer systems of interest. In this context, it is desired to have a generally applicable analysis protocol to analyze lipid order states and domains in lipid bilayer simulations, which does not suffer from the aforementioned issues. For such an analysis protocol, it is required to have a physically clear descriptor of lipid order states that reflects the local environment of lipids robustly. Among many such observables, the deformation of a cylinder representation of a lipid can be a simple and intuitive descriptor of lipid order states as elaborated below. Let us consider a lipid in a membrane environment. Depending on its interactions with (heterogeneously distributed) different lipid types, the lipid can be either loosely or
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
tightly packed with its neighboring ones. When the lipid is tightly packed, its surface area (A) becomes smaller and chain thickness (d) gets larger than those in loosely packed environments. Each observable reflects the local environment, so that a descriptor using both A and d should be robust. Therefore, the deformation of lipids can describe the lipid order states on a clear physical basis. In this approach, because the characteristics of lipids differ from one lipid type to the others, it is natural to determine lipid order states of a given lipid type independently from the other types. Therefore, the lipid order states for each lipid can be determined, from whose spatial distribution the domains can be assigned. In this work, we follow the above approach and propose a two-step analysis protocol to characterize lipid order states and domains based on the observable order states of A and d (see Figure 1). In the first step of the protocol, the lipid order states of individual lipids at each time are determined from the distributions of A and d of each lipid type, which are robustly inferred by a hidden Markov model (HMM) analysis.20,21 A similar approach was applied for local composition to assign lipid domain.16 In the next step, the lipid domains are assigned as clusters of high order state lipids determined by the Getis-Ord local spatial autocorrelation (Gi*) statistics22,23 of the order states. Our analysis protocol is conceptually simple yet robust due to the employed statistical methods (HMM and Gi* statistics). The analysis protocol was then applied to analysis of quinary mixed bilayers of CHOL, DMPC, DMPE, POPC, and POPE, which are challenging systems (see below). Using our analysis protocol, the size distributions and compositions of leaflet ordered lipid clusters were obtained, from which the local bilayer properties were analyzed.
ACS Paragon Plus Environment
Page 4 of 22
Page 5 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Figure 1. The two-step workflow for analysis of lipid order states and domains in lipid bilayer simulations. In the first step, the time series of lipid order state {O(t)} are obtained from the (hidden) order state of A and d, {SA(t)} and {Sd(t)}, which are inferred by a hidden Markov model (HMM) analysis of {A(t)} and {d(t)}. In the second step, ordered-lipid clusters are assigned by the Getis-Ord local spatial autocorrelation (Gi*) statistics of the Voronoi tessellation of lipids, on which lipid order states are mapped. Analysis protocol for lipid order states and domains The first step of our analysis procedure is the determination of the lipid order states using A and d, which are inferred from the HMM analysis of the time series of A and d. An HMM is well suited for a robust inference of the statistically most probable sequences of order states for a set of observables {X} that are generally noisy time series. As the terminology implies, (i) the hidden order states are not directly observable and the probabilities of the observable emission states from these hidden states are modeled by an emission probability matrix and (ii) the evolution of the hidden state is assumed to be Markovian (independent of history other than just the previous hidden state) modeled by a transition matrix between hidden states. Though there is a variety of available HMM models (continuous HMM,24 multi-variate Gaussian HMM,25 etc.), in this work, we will focus on the simplest HMM, i.e., two univariate discrete HMMs for A and d, from which the lipid order states are determined (see below). The underlying assumption in the HMM analysis is that P({X}) for each lipid type can be decomposed into those from low and high order states, which is approximated by a bi-Gaussian distribution, PG(X) (Figure 2): 𝑃G(𝑋) = 𝑃L(𝑋) + 𝑃H(𝑋) = 𝛼L𝒩(𝑋│𝜇𝑋,L,𝜎𝑋,L) + 𝛼H𝒩(𝑋│𝜇𝑋,H,𝜎𝑋,H)
(1)
where PL(X) and PH(X) are the distributions for low and high order-states, 𝛼L and 𝛼H are coefficients, and 𝒩(𝑥│𝜇,𝜎) is a normal distribution of a variable x with the mean 𝜇 and the standard deviation σ. In the current analysis scheme, the mean values are used as the boundary values of the low and high order states, which are utilized in the setup of (discrete) HMM along with their standard deviations (see below). When P(X) for a given observable is indecomposable (though it is rare for lipids), we assign the observable into the high order state.
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Figure 2. Bi-Gaussian approximation of (A) P(A) and (B) P(d) for DMPC (in black lines) from an all-atom simulation of a quinary mixed bilayer consisting of Chol (10 mol %) and an equimolar (22.5 mol %) mixture of DMPC, POPC, DMPE, and POPE. In general, the simulation distributions are well approximated by the sum of the low (PL(X)) and high (PH(X)) order-state Gaussian distributions. In this work, the discrete HMMs for A and d consist of nine emission states and three hidden states: low (L), intermediate (I), and high (H). The (discrete) emission states are assigned by partitioning the observable space into nine sub-spaces (defined as 0,1, 2, ⋯, 8), where a larger value of emission state implies a higher probability of its hidden state being H state. For d, the lowest and highest emission states (0 and 8) are assigned when d < dL (= d,L) and d ≥ dH (= d,H). The intermediate emission state i (i = 1, 2, ⋯, 7) are assigned when 𝑑L + (𝑖 ― 1)𝛥𝑑 ≤ 𝑑 < 𝑑L +𝑖𝛥𝑑, where 𝛥𝑑 = (𝑑H ― 𝑑L)/7. Similarly, for A, the lowest and highest emission states are assigned when A > AL (= A,L) and A ≤ AH (= A,H), while the intermediate emission state i are assigned when 𝐴L ― (𝑖 + 1)𝛥𝐴 ≥ 𝐴 > 𝐴L ―𝑖𝛥𝐴, where 𝛥𝐴 = (𝐴L ― 𝐴H)/7. It is known that the result of an HMM is not so sensitive to its initial transition probability matrix and thus it has been typically assigned as a uniform matrix.20 However, the HMM is sensitive to its initial emission probability matrix.20 Therefore, to obtain reliable results from the HMM, the initial emission probability matrix should be carefully set up. We make use of the integrals of 𝒩(𝑋│𝜇𝑋,L,𝜎𝑋,L) and 𝒩(𝑋│𝜇𝑋,H,𝜎𝑋,H) over the subspaces of the observables to assign the initial emission probabilities for the low and high hidden ordered states, respectively. For the intermediate ordered state, we make use of the integral of 𝒩(𝑋│𝜇𝑋,I,𝜎𝑋,I), where we define 𝜇𝑋,I = (𝜎𝑋,H𝜇𝑋,L + 𝜎𝑋,L𝜇𝑋,H)/(𝜎𝑋,L + 𝜎𝑋,H) and 𝜎𝑋,I = min {|𝜇𝑋,H ― 𝜇𝑋,I|,|𝜇𝑋,L ― 𝜇𝑋,I|}/3 to make its initial emission probabilities be localized within 3X,I around X,I in between the upper (X,H) and lower (X,L) bounds. With these initial transition and emission probability matrices, the HMM for observable X is trained by the Baum-Welch algorithm20,21 using the input time series of the emission states for each lipid type in a given leaflet. Then, the times series of the X order state (SX) for each lipid is decoded by Viterbi algorithm26 to obtain its most likely order state sequence.
ACS Paragon Plus Environment
Page 6 of 22
Page 7 of 22
From the obtained observable order state pairs (a total of nine possible combinations of SA-Sd pairs), the lipid order states can be assigned by considering the deformation and tilt of a lipid (Figure 3). We start by defining three canonical order states OH, OI, and OL for an upright cylinder, whose SA-Sd pairs are (H, H), (I, I), and (L, L), respectively. These canonical order states undergo transitions by concerted changes in SA and Sd (elastic deformation). The neighboring four combinations around these canonical order states (along the axial direction) can be accessed by changes in one of SA and Sd (partial deformation), which we define as intermediate order states (OHI between OH and OI and OLI between OI and OL). Before considering remaining combinations, we note that OHI and OLI remain to be the same by the changes of SA and Sd in opposite directions (e.g., (I, L) (L, I)). Physically, such SA-Sd pair changes can be interpreted as tilt of a cylinder. Then, by considering the tilt of canonical OI state, the order state of (H, L) can be assigned to OI. While physically unrealistic, for completeness, we define the last combination of (L, H) as OI state. By this scheme, all nine combinations of SA-Sd pairs are mapped to five lipid order states (OH, OHI, OI, OLI, and OL). For more robust inference of lipid order state time series, one can assign it from the (exponential) running averages of observable order states over a given interval (e.g., the past 1 ns in this work).
H
Thickness order, Sd
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Tilt
Tilt
Tilt
Tilt
I
L
L
I Area order, SA
H
Figure 3. A pictorial representation of the mapping scheme of SA-Sd pairs to lipid order states. A lipid is shown as a cylinder, and its order states are indicated in different colors: high (OH in blue), high-intermediate (OHI in green), intermediate (OI in yellow), lowintermediate (OLI in orange), and low (OL in red). Two diagonal directions of SA-Sd pair changes correspond to elastic deformation and tilt of a cylinder. The changes in only one of SA or Sd can be interpreted as partial deformation. The second step in our analysis protocol in Figure 1 is assignment of domains by clustering the ordered lipids from the spatial distribution of their order states obtained in the previous step. First, the spatial distribution of the order states is prepared by mapping these onto the Voronoi tessellation of lipids. The order state map is then analyzed by the
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 22
Getis-Ord local spatial autocorrelation (Gi*) statistics22,23 whose standardized form (zscore) for a lipid i, Gi*, can be written as 𝐺𝑖∗ =
∑𝑗𝑤𝑖𝑗𝑂𝑗 ― 𝑊𝑖∗ 𝑂
[
]
2
𝑠 (𝑛𝑆1𝑖∗ ― 𝑊𝑖∗ )/(𝑛 ― 1)
1/2
(2)
where n is the number lipids in the leaflet, 𝑤𝑖𝑗 = 𝑏𝑖𝑗/𝑑𝑖𝑗 is a symmetric spatial weight matrix with bij and dij being the shared border and distance between a lipid pair, i and j, Oi represents the order state of lipid i, 𝑊𝑖∗ = ∑𝑗𝑤𝑖𝑗, 𝑂 and s denote usual sample mean and variance, and 𝑆1𝑖∗ = ∑𝑗𝑤2𝑖𝑗. In the calculation of Gi*, to consider the self-influence, we assign wii as an average weight to its neighbors, 𝑤𝑖 = ∑𝑗𝑤𝑖𝑗/𝑛𝑖, where j belongs to the contacting neighbors of i and ni is the number of these neighbors. Gi* is a measure of clustering of the high and low values of O around a lipid i, whose interpretation is analogous to a standard z-score (for a normal distribution). As Gi* deviates more from 0, the clustering of high or low values of O is more intense (i.e., hot/cold spots). In the original Gi* statistics, the hot and cold spots are determined relative to the mean values of a given spatial distribution. Thus, the original Gi* statistics may introduce systematic errors in domain analysis, which becomes more severe as 𝑂 deviates more from OI. For example, for a disordered membrane, 𝑂 is smaller than OI, so that the hot spots (i.e., ordered regions) assigned by Eq. (2) may contain lipids with lower order states than OI, resulting in overestimated hot spots. In the opposite case, where the membrane is highly ordered, the hot spots may be underestimated. To remedy such systematic errors, let us consider the deviation of 𝑂 from OI around lipid i by expressing in an analogous form to Gi*: 𝐺𝑖∗
=
𝑊𝑖∗ (𝑂 ― 𝑂I)
[
]
2
𝑠 (𝑛𝑆1𝑖∗ ― 𝑊𝑖∗ )/(𝑛 ― 1)
1/2
(3)
Eq. (3) is a z-core that describes how significant the system’s O globally deviated from the reference state OI. As a simple remedy, we define the modified Gi* by adding Eq. (3) to Eq. (2) as ∗
Γ𝑖∗ = 𝐺𝑖∗ + 𝐺𝑖
(4)
where the z-score of 𝑂 for lipid i (in a bilayer whose global order state is OI) is added as a correction term. In other words, i* is Gi* calculated relative to an intermediate order state OI instead of 𝑂. Using the modified Gi* (i.e., i*) statistics in Eq. (4), the lipids can be clustered as follows. First, the lipids with i* > 1.645 (90% confidence interval) are assigned as members of cluster core, and those with i* < -1 (likely to be low ordered) are not considered further in the member assignment. Among remaining unassigned lipids, those with O > OI and contacting sufficiently to the core members are assigned to additional cluster members.
ACS Paragon Plus Environment
Page 9 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Here, we define a sufficient contact of lipid i to core members when ∑𝑘 ∈ Core𝑤𝑖𝑘 > 𝑤𝑖, where lipid k belongs to core. These tentatively assigned ordered lipids are then clustered hierarchically as follows. Initially, all members are assigned to different clusters. Analogous to a contact between a lipid and core members, a pair of clusters CM and CN are merged when these are in a sufficient contact, i.e., if there exists i in CM with ∑𝑘 ∈ 𝐶 𝑤𝑖𝑘 > 𝑤𝑖 or j in CN with ∑𝑘 ∈ 𝐶 𝑤𝑗𝑘 > 𝑤𝑗. The clustering iterates until there remains N
M
no more contacting cluster. We note that the size of an assigned cluster is bounded by the total number of lipids in a leaflet, while the cluster may be connected to its own images and extend across the primary system due to a nature of the periodic boundary conditions. A similar method was applied to bilayer simulations of ganglioside GM1 and POPC.27 Illustrative example To illustrate the usefulness of our analysis protocol, we used quinary mixed bilayers composed of equimolar mixture of DMPC, DMPE, POPC, and POPE with three different Chol concentrations: 10 mol % (Chol10%), 20 mol % (Chol20%), and 30 mol % (Chol30%). Because any type of the phospholipids in these mixed bilayers would not show strong preference over the others due to subtle difference in their physical properties (chain saturation and head groups), these systems would be challenging to analyze by either the local composition or order parameter approach. A
B
C
Figure 4. The last snapshots of quinary mixed all-atom bilayers composed of equimolar mixture of DMPC (blue), DMPE (cyan), POPC (pink), and POPE (red) with three different Chol (green) concentrations: (A) 10 mol % (Chol10%), (B) 20 mol % (Chol20%), and (C) 30 mol % (Chol30%). The phospholipids are shown in sticks and cholesterols are shown in spheres. For clarity, the head groups, ions, and water were omitted. We first considered the all-atom (AA) simulations (Figure 4). Each system was composed of a total of 200 lipid and Chol molecules in each leaflet with bulk water and 0.15M KCl with initial sizes of ~110×110×80 Å3. At each Chol concentration, we also considered coarse-grained (CG) models with three different system sizes: 1X (same number of lipids to the corresponding AA systems), 4X (4 times larger), and 16X (16 times larger). These systems are denoted by adding a superscript to that for the corresponding AA system (e.g., Chol10%-4X for four times larger CG simulation system at 10 mol % Chol). The system information is summarized in Table S1 in Supporting
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 22
information (SI). The systems were generated using Membrane Builder28–30 (AA) and Martini Maker31,32 (CG) in CHARMM-GUI.33 We then preformed 6-s simulations for AA models using ANTON34 and 10-s (1X and 4X) and 7-s (16X) simulations for CG models using GROMACS35 under constant pressure and temperature (NPT) conditions at p = 1 atm and T = 303.15 K (see S1.1 Simulation details in SI for details). The last 4-s trajectories were analyzed for AA simulations and the last 7-s (1X and 4X) and 4-s (16X) trajectories were analyzed for CG simulations. Our analysis protocol begins with obtaining the time series of A and d of individual lipids and Chol, as well as their distributions P(A) and P(d) for each lipid type and Chol as inputs for the HMM (Figure 2 and see S1.2. Surface area and hydrophobic thickness in SI for details). A and d were moderately anti-correlated (Figure S1 in SI and see Table S2 for the Pearson correlation coefficients). This confirms that A and d can jointly describe the lipid order state better than a single observable (either A or d) as the uncertainty from one observable can be corrected by the other. Before proceeding our analysis, we checked the applicability of the local composition approach to our systems, for which we analyzed the nearest neighbors (NN) for each component (Tables 1 and S3-S6). The difference between each lipid type was most clearly shown in the AA model at the lowest concentration (Chol10%, Table 1), where the number of NNs around Chol was decreased slightly in the order of DMPC > DMPE POPC POPE, which is consistent with both the template36 and umbrella mechanism37 of Chol-lipid interactions.38 This trend became diminished at higher Chol concentrations (Table S3), which we attributed to the smaller number of available lipids around Chol that can preferentially interact with. In the CG models, there was no preference of a specific lipid type over the others (Tables S4-S6). Among the phospholipids, we did not observe strong aggregation of a specific lipid type, as the number of NNs were similar for any pair of lipid types. Overall, there was no strongly preferred lipid type around Chol and no strong aggregation of a specific lipid type. Thus, the local composition approach based on a specific lipid type (other than Chol) cannot be applied to the domain analysis for these systems. Table 1. Number of the nearest neighbors around each component in Chol10%.† NN type Chol DMPC DMPE POPC POPE †Average
Chol 0.3 1.3 1.1 1.1 1.0
Lipid type DMPC DMPE POPC 0.6 0.5 0.5 1.4 1.4 1.4 1.4 1.4 1.5 1.4 1.5 1.4 1.4 1.4 1.5
POPE 0.4 1.4 1.4 1.5 1.4
over 1-s blocks with standard errors < 0.1.
The next check was for the order parameters, for which we calculated two types of commonly used order parameters, the liquid crystal order parameter (P2) and bondorientational (hexatic) order parameter (6).39 Briefly, P2 is based on the chain
ACS Paragon Plus Environment
Page 11 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
orientations and 6 represents hexatic chain packing with neighbors (see S1.3 Order parameters in SI for details). As shown in Figure 5, there seem to be little correlation between them, which was quantified by the Pearson correlation coefficient in Table 2. Such little correlation again raises a critical issue in choosing an appropriate order parameter in domain analysis. These order parameters were further examined by calculating the Pearson coefficient between the lipid order state (O) obtained in our analysis. As shown in Table 2, the P2 and O were moderately correlated at all Chol concentrations, whereas 6 and O were very weakly correlated even in Chol10%, implying that 6 is not an appropriate order parameter for our system. Consistent correlation patterns were observed for CG models (Tables S7-S9). Considering that P2 only considers the orientation of lipid tails, the lipid order states in our approach (from both A and d) would be more generally applicable in the analysis of complex bilayer simulations. Below, we present the results from our analysis
Figure 5. The last snapshots of |P2|, |6|, and the assigned ordered-lipid states for (A) Chol10%, (B) Chol20%, and (C) Chol30% (AA). |P2| and |6| were calculated for each chain and their values at each point on the xy-plane were estimated using a Gaussian filter. The assigned lipid order states (O) were mapped on the Voronoi tessellation with the same color code as in Figure 3, where the assigned cluster members were indicated as black borders. The primary cell (simulation box) was indicated by a red box. Table 2. Pearson correlation coefficients between |P2|, |6|, and O. Modela
System
Component
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Chol |P2|-|6|
|P2|-O
|6|-O
aThe
Page 12 of 22
DMPC
DMPE
POPC
POPE
Alla
Chol10%
0.01
0.03
0.03
0.03
0.04
0.03
Chol20%
0.01
0.02
0.03
0.02
0.04
0.02
Chol30%
0.01
0.02
0.03
0.01
0.03
0.02
Chol10%
–
0.55
0.57
0.55
0.56
0.55
Chol20%
–
0.50
0.55
0.49
0.55
0.51
Chol30%
–
0.35
0.43
0.36
0.44
0.38
Chol10%
–
0.02
0.02
0.01
0.02
0.02
Chol20%
–
0.00
0.01
0.01
0.02
0.01
Chol30%
–
-0.01
0.01
-0.01
0.00
0.00
Pearson correlation coefficients were calculated for all lipids.
In the second step of our analysis, the spatial distributions of lipid order states were analyzed using the i* statistics in Eq. (4) to obtain the ordered-lipid clusters. For the lipid types considered in the present work (DMPC, DMPE, POPC, and POPE), the main driving force for lipid order state and domain formation are Chol-lipid interactions. These interactions arise from hydrophobic interactions between sterol rings and lipid tails,36 whose strengths would be similar among these lipid types40 (also shown as no preference of any lipid type over the others, Table 1). Thus, a Chol could easily change its interaction partner from one lipid type to another, so that the leaflet domains were observed to be highly dynamic (i.e., they kept forming and breaking depending on the local Chol distributions (and/or diffusion), see also Movies S1-S3). Though the domains are highly dynamic, as shown in Figure 5, even at the lowest Chol concentration, ordered lipids appeared to be aggregated. Such aggregations formed bigger clusters at higher Chol concentrations, as also shown in the probability (P(f)) of a lipid belonging to f-size cluster where f is the fractional cluster size normalized by the total number of lipids in the leaflet (Figure 6A). The Chol content in the clusters was increased compared to its bulk concentration, while the relative amounts of DMPC, DMPE, POPC, and POPE were similar (close to equimolar) (Figure 6B). This result again shows that any type of phospholipid is not preferred to be enriched in these clusters over the other types. As expected, the Chol content is decreased outside the ordered-lipid cluster with comparable relative amounts among DMPC, DMPE, POPC, and POPE.
ACS Paragon Plus Environment
Page 13 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Figure 6. (A) Probability of a lipid (phospholipid or Chol) belonging to a f-size cluster, P(f) with standard errors (gray area) for AA model. (B) Lipid compositions of orderedlipid clusters and disordered regions, which are indicated as O and D on top of each bar. The fraction of each component was shown in different colors: green (Chol), blue (DMPC), cyan (DMPE), pink (POPC), and red (POPE). The standard errors were shown as error bars. Considering rather weak and similar interactions between different lipid types and Chol, the ordered-lipid clusters in one leaflet are likely to be anti-correlated with those in the other leaflet at low Chol concentration to minimize the stress on the bilayer (Figure 7A). At higher Chol concentration, however, domains in opposing leaflet become more positively correlated because there would be less available for disordered regions to avoid the overlap between the ordered-leaflet domains (Figures 7B and 7C). The correlation between these domains were analyzed by using the results of leaflet domain analysis as follows. For a molecule (lipid or Chol) in a given leaflet, its nearest molecule in the other leaflet was assigned by using the xy-COM of the representative atoms (see S1.2. Surface area and hydrophobic thickness in SI).27,41 Then, for the chosen lipid pair of molecules (between the leaflets), we categorized the local bilayer order states into three classes: ordered (O, both lipids from ordered-lipid clusters), intermediate (I, only one lipid from ordered-lipid cluster), and disordered (D, none of them from the ordered-lipid clusters).27 The fraction of lipids in O-region (fO) and the fraction of lipids in leaflet ordered-lipid clusters belonging to O-region (fOO) are consistent with the aforementioned picture. As shown in Figure 8, fO fluctuated around ~ 0.2 and fOO around ~ 0.4 at 10 mol % Chol concentration throughout simulation, where the low fOO reflects the anti-correlation between domains in the opposing leaflets. At higher Chol concentration, fO fluctuates around ~ 0.3 and fOO around ~ 0.6 at 20 mol % Chol, and fO around ~ 0.5 and fOO around ~ 0.7 at 30 mol % Chol. For the assigned bilayer regions, the hydrophobic thickness (dH) were estimated (Table 3), where one can see the ordered regions are thicker than the disordered regions. At 30 mol % Chol, dH of the D-region is comparable to that of the Oregion in Chol10%. Thus, it appears that Chol30% system is globally ordered, where our analysis shows the sub-structures inside the lipid domain. Table 3. Hydrophobic thickness of bilayer regions, dH (Å)† System
Bilayer regions O I D
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Chol10% Chol20% Chol30% †Average
31.7 32.8 33.5
29.9 31.5 32.8
28.0 30.0 32.2
over 1-s blocks. The standard errors were < 0.1.
Figure 7. The last snapshots of the assigned ordered-lipid states for (A) Chol10%, (B) Chol20%, and (C) Chol30% (AA). The left and right panels are the snapshots of upper and lower leaflets, respectively. The assigned lipid order states (O) were mapped on the Voronoi tessellation with the same color code as in Figure 3, where the assigned cluster members were indicated as black borders. The primary cell (simulation box) was indicated by a red box.
ACS Paragon Plus Environment
Page 14 of 22
Page 15 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Figure 8. The fraction of lipids in O-region (fO, black) and the fraction of lipids in leaflet ordered-lipid clusters belonging to O-region (fOO, red) for (A) Chol10%, (B) Chol20%, and (C) Chol30% (AA). The same analyses for the CG simulations show consistent results with the AA simulations in terms of the correlation between lipid order state (O) and other order parameters (|P2| and |6|) (Figure S2-S4), the spatial distribution and dynamics of the ordered-lipid domains (Movies S4-S12), and their lipid compositions (Figure S5). The size distribution P(f) (Figure S5) were also consistent with those of AA model except those for the lowest Chol concentration, where the peak of P(f) (i.e., f ~ 0.6 in Chol10%-1X) became weaker and shifted to f ~ 0.5 in Chol10%-4X and disappeared in Chol10%-16X. The system size dependence of the mean cluster size is shown in Figure 9, where the mean cluster size at 10 mol % Chol grows sub-linearly as a function of the total number of lipids and Chol (N) compared to those at 20 and 30 mol %. These results imply that the clusters at Chol 10 mol % would not aggregate into a big cluster. For systems at 20 and 30 mol % Chol, there is one notable change in size distributions, which became sharper and their peaks were shifted to smaller f-size for larger systems. It appears that the behavior is related to the increased global order for larger systems as reported in a previous CG MD study of lipid domains for ternary mixed bilayer composed of Chol, DIPC, and DOPC.13
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Figure 9. System size (N) dependence of the average size of ordered-lipid clusters from CG simulations with their fit to N, where N is the number of lipids and Chols in a leaflet and is an exponent. The error bars are standard errors over 1-s blocks. The correlation between the domains in opposing leaflets in the CG model were consistent with the AA model at all sizes (Figure S6), where fluctuation of fO and fOO were decreasing with increasing system size (due to a greater number of components in the system). Note that, while fO and fOO are similar for different system sizes, the actual domain sizes are increasing with the system size (Figures S7-S9). However, the behavior of estimated dH was inconsistent with those obtained from the AA model and the concept of lipid domains. As shown in Table S10, dH values were comparable in all bilayer regions in Chol10% systems and slightly decreased in ordered bilayer regions in Chol20% and in Chol30% systems. To check the validity of our analyses, we also analyzed A and d of the bilayer regions, which behaved consistently with those from the AA model. Therefore, the inconsistent dH statistics appears to originate from the bead selection for dH calculations in CG models (ROH for Chol and GL1 and GL2 for phospholipids). The GL1 and GL2 beads are positioned above the ROH bead, so the local bilayer thickness in a region of high Chol concentration tended to be smaller than that in a region of low Chol concentration. Although dH depended on the bead selection (which needs a separate care, so that the reference beads for each types of CG lipids closely matches those of corresponding AA lipids), the lipid domains can be analyzed without difficulty in our analysis protocol. Conclusions In this study, we propose a robust method for analysis of lipid order states and domains in lipid bilayer simulations using A and d of lipids, where lipid order states are assigned by considering deformation of a lipid at different packing environments. To illustrate the
ACS Paragon Plus Environment
Page 16 of 22
Page 17 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
usefulness of our approach, we analyzed simulations of quinary mixed bilayer simulations consisting of Chol, DMPC, DMPE, POPC, POPE. For these mixed bilayer systems, the nearest neighbor statistics based on the local composition approach failed to provide proper lipid order state information. In addition, the two widely used order parameters, P2 (chain orientation) and 6 (chain packing), were not consistent to each other (low Pearson correlation coefficients) for these challenging systems. The lipid order states from our method were moderately correlated to the results based on P2, whereas it was not correlated with 6. The lipid domains were assigned by the modified Gi* (i.e., i*) statistics of the lipid order states, which appeared to aggregate even at the lowest Chol concentration (10 mol %) and formed bigger clusters at higher Chol concentrations. The composition of each phospholipid type in the clusters were comparable to each other, which again showed the inapplicability of the local composition approach. The leaflet domains do not appear to be co-localized at low Chol concentration and the increase in leaflet-leaflet correlation for our model bilayer is likely due to the increased ordered domains. Our analysis protocol is conceptually simple yet robust due to the employed statistical methods (HMM and Gi* statistics), where order states of each lipid types are analyzed independently. Thus, our method is general and can be applied to various complex bilayer simulations. ASSOCIATED CONTENT Supporting Information The Supporting Information is available free of charge. Section S1, Methods; Table S1, System information; Table S2, Pearson correlation coefficients between A and d; Tables S3-S6, The number of the nearest neighbors around each component; Table S7-9, Pearson correlation coefficients between |P2| and |6|, |P2| and O, and |6| and O of coarse-grained systems; Table S10, Hydrophobic thickness of bilayer regions dH in coarse-grained systems; Table S11-13, Average surface area (A) and chain hydrophobic thickness (d) in coarse-grained systems; Figure S1, Two dimensional histograms of A and d, P(A,d); Figures S2-4. The last snapshots of |P2|, |6|, and the assigned ordered-lipid states of coarse-grained systems; Figure S5, Probability of a lipid belonging to a f-size cluster, P(f) and lipid compositions of ordered-lipid clusters and disordered regions; Figure S6, The fraction of lipids in O-region (fO) and the fraction of lipids in leaflet ordered-lipid clusters belonging to O-region (fOO) of coarse-grained systems; Figures S7-S9, The last snapshots of the assigned ordered-lipid states of coarse grained-systems; Movies S1-S12, Dynamics of ordered-lipid clusters. AUTHOR INFORMATION Corresponding Author *(S. P.) Tel: +1-610-758-4707. Fax: +1-610-758-4004. E-mail:
[email protected]. *(W. I.) Tel: +1-610-758-4524. Fax: +1-610-758-4004. E-mail:
[email protected].
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Notes The authors declare no competing financial interests. ACKNOWLEDGEMENT This work was supported by NSF MCB-1727508 and XSEDE MCB070009. Anton 2 computer time was provided by the Pittsburgh Supercomputing Center (PSC) through Grant R01GM116961 from the National Institutes of Health. The Anton 2 machine at PSC was generously made available by D.E. Shaw Research.
ACS Paragon Plus Environment
Page 18 of 22
Page 19 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
References (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11)
(12) (13)
(14) (15) (16)
Alberts, B.; Johnson, A.; Lewis, J.; Raff, M.; Roberts, K.; Walter, P. Molecular Biology of the Cell, 6th ed.; Garland Science: New York, 2014. Lingwood, D.; Simons, K. Lipid Rafts As a Membrane- Organizing Principle. Science (80-. ). 2010, 327 (5961), 46–50. Simons, K.; Van Meer, G. Lipid Sorting in Epithelial Cells. Biochemistry 1988, 27 (17), 6197–6202. Stöckl, M. T.; Herrmann, A. Detection of Lipid Domains in Model and Cell Membranes by Fluorescence Lifetime Imaging Microscopy. Biochim. Biophys. Acta - Biomembr. 2010, 1798 (7), 1444–1456. Pan, J.; Heberle, F. A.; Petruzielo, R. S.; Katsaras, J. Using Small-Angle Neutron Scattering to Detect Nanoscopic Lipid Domains. Chem. Phys. Lipids 2013, 170– 171, 19–32. Sanchez, J.; Badia, A. Atomic Force Microscopy Studies of Lateral Phase Separation in Mixed Monolayers of Dipalmitoylphosphatidylcholine and Dilauroylphosphatidylcholine. Thin Solid Films 2003, 440 (1–2), 223–239. Richter, R. P.; Brisson, A. R. Following the Formation of Supported Lipid Bilayers on Mica: A Study Combining AFM, QCM-D, and Ellipsometry. Biophys. J. 2005, 88 (5), 3422–3433. Aliaskarisohi, S.; Tierno, P.; Dhar, P.; Khattari, Z.; Blaszczynski, M.; Fischer, T. H. M. On the Diffusion of Circular Domains on a Spherical Vesicle. J. Fluid Mech. 2010, 654, 417–451. Larsen, J.; Hatzakis, N. S.; Stamou, D. Observation of Inhomogeneity in the Lipid Composition of Individual Nanoscale Liposomes. J. Am. Chem. Soc. 2011, 133 (28), 10685–10687. Sonnleitner, A.; Schütz, G. J.; Schmidt, T. Free Brownian Motion of Individual Lipid Molecules in Biomembranes. Biophys. J. 1999, 77 (5), 2638–2642. Przybylo, M.; Sýkora, J.; Humpolíčová, J.; Benda, A.; Zan, A.; Hof, M. Lipid Diffusion in Giant Unilamellar Vesicles Is More than 2 Times Faster than in Supported Phospholipid Bilayers under Identical Conditions. Langmuir 2006, 22 (6), 9096–9099. Pandit, S. A.; Jakobsson, E.; Scott, H. L. Simulation of the Early of Nano-Domain Formation in Mixed Bilayers of Sphingomyelin, Cholesterol, and Dioleylphosphatidylcholine. Biophys. J. 2004, 87 (5), 3312–3322. Pantelopulos, G. A.; Nagai, T.; Bandara, A.; Panahi, A.; Straub, J. E. Critical Size Dependence of Domain Formation Observed in Coarse-Grained Simulations of Bilayers Composed of Ternary Lipid Mixtures. J. Chem. Phys. 2017, 147 (9), 095101. Baoukina, S.; Mendez-Villuendas, E.; Tieleman, D. P. Molecular View of Phase Coexistence in Lipid Monolayers. J. Am. Chem. Soc. 2012, 134 (42), 17534–17553. Baoukina, S.; Mendez-Villuendas, E.; Bennett, W. F. D.; Tieleman, D. P. Computer Simulations of the Phase Separation in Model Membranes. Faraday Discuss. 2013, 163, 63–75. Sodt, A. J.; Sandar, M. L.; Gawrisch, K.; Pastor, R. W.; Lyman, E. The Molecular Structure of the Liquid-Ordered Phase of Lipid Bilayers. J. Am. Chem. Soc. 2014,
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
(17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30)
(31) (32)
136 (2), 725–732. Sodt, A. J.; Pastor, R. W.; Lyman, E. Hexagonal Substructure and Hydrogen Bonding in Liquid-Ordered Phases Containing Palmitoyl Sphingomyelin. Biophys. J. 2015, 109 (5), 948–955. Engberg, O.; Hautala, V.; Yasuda, T.; Dehio, H.; Murata, M.; Slotte, J. P.; Nyholm, T. K. M. The Affinity of Cholesterol for Different Phospholipids Affects Lateral Segregation in Bilayers. Biophys. J. 2016, 111 (3), 545–556. Williams, J. A.; Wassall, C. D.; Kemple, M. D.; Wassall, S. R. An Electron Paramagnetic Resonance Method for Measuring the Affinity of a Spin-Labeled Analog of Cholesterol for Phospholipids. J. Membr. Biol. 2013, 246 (9), 689–696. Baum, L. E.; Petrie, T.; Soules, G.; Weiss, N. A Maximization Technique Occuring in Statistical Analysis of Probabilistic Functions of Markov Chains. Ann. Math. Stat. 1970, 41 (1), 164–171. Rabiner, L. R. A Tutorial on Hidden Markov-Models and Selected Applications in Speech Recognition. Proc. Ieee 1989, 77 (2), 257–286. Getis, A.; Ord, J. K. The Analysis of Spatial Association by Use of Distance Statistics. Geogr. Anal. 1992, 24 (3), 189–206. Ord, J. K.; Getis, A. Local Spatial Autocorrelation Statistics: Distributional Issues and an Application. Geogr. Anal. 1995, 27 (4), 286–306. Leggetter, C. J.; Woodland, P. C. Maximum Likelihood Linear Regression for Speaker Adaptation of Continuous Density Hidden Markov Models. Comput. Speech Lang. 1995, 9 (2), 171–185. Gauvain, J.-L. Maximum a Posteriori Estimation for Multivariate Gaussian Mixture Observations of Markov Chains. IEEE Trans. Speech Audio Process. 1994, 2 (2), 291–298. Viterbi, A. J. Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm. IEEE Trans. Inf. Theory 1967, 13 (2), 260–269. Patel, D. S.; Park, S.; Wu, E. L.; Yeom, M. S.; Widmalm, G.; Klauda, J. B.; Im, W. Influence of Ganglioside GM1 Concentration on Lipid Clustering and Membrane Properties and Curvature. Biophys. J. 2016, 111 (9), 1987–1999. Jo, S.; Kim, T.; Im, W. Automated Builder and Database of Protein/Membrane Complexes for Molecular Dynamics Simulations. PLoS One 2007, 2 (9), e880. Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI Membrane Builder for Mixed Bilayers and Its Application to Yeast Membranes. Biophys. J. 2009, 97 (1), 50–58. Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; Davila-Contreras, E. M.; Qi, Y. F.; Lee, J. M.; Monje-Galvan, V.; Venable, R. M.; Klauda, J. B.; Im, W. CHARMM-GUI Membrane Builder Toward Realistic Biological Membrane Simulations. J. Comp. Chem. 2014, 35 (27), 1997–2004. Qi, Y.; Ingólfsson, H. I.; Cheng, X.; Lee, J.; Marrink, S. J.; Im, W. CHARMMGUI Martini Maker for Coarse-Grained Simulations with the Martini Force Field. J. Chem. Theory Comput. 2015, 11 (9), 4486–4494. Hsu, P. C.; Bruininks, B. M. H.; Jefferies, D.; Cesar Telles de Souza, P.; Lee, J.; Patel, D. S.; Marrink, S. J.; Qi, Y.; Khalid, S.; Im, W. CHARMM-GUI Martini Maker for Modeling and Simulation of Complex Bacterial Membranes with Lipopolysaccharides. J. Comput. Chem. 2017, 38 (27), 2354–2363.
ACS Paragon Plus Environment
Page 20 of 22
Page 21 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
(33) (34)
(35) (36) (37) (38) (39) (40) (41)
Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A Web-Based Graphical User Interface for CHARMM. J. Comput. Chem. 2008, 29 (11), 1859–1865. Shaw, D. E.; Grossman, J. P.; Bank, J. A.; Batson, B.; Butts, J. A.; Chao, J. C.; Deneroff, M. M.; Dror, R. O.; Even, A.; Fenton, C. H.; Forte, A.; Gagliardo, J.; Gill, G.; Greskamp, B.; Ho, C. R.; Ierardi, D. J.; Iserovich, L.; Kuskin, J. S.; Larson, R. H.; Layman, T.; Lee, L. S.; Lerer, A. K.; Li, C.; Killebrew, D.; Mackenzie, K. M.; Mok, S. Y. H.; Moraes, M. A.; Mueller, R.; Nociolo, L. J.; Peticolas, J. L.; Quan, T.; Ramot, D.; Salmon, J. K.; Scarpazza, D. P.; Ben Schafer, U.; Siddique, N.; Snyder, C. W.; Spengler, J.; Tang, P. T. P.; Theobald, M.; Toma, H.; Towles, B.; Vitale, B.; Wang, S. C.; Young, C. Anton 2: Raising the Bar for Performance and Programmability in a Special-Purpose Molecular Dynamics Supercomputer. IEEE 2014, 41–53. Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindah, E. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25. Sugahara, M.; Uragami, M.; Yan, X.; Regen, S. L. The Structural Role of Cholesterol in Biological Membranes. J. Am. Chem. Soc. 2001, 123 (32), 7939– 7940. Huang, J.; Feigenson, G. W. A Microscopic Interaction Model of Maximum Solubility of Cholesterol in Lipid Bilayers. Biophys. J. 1999, 76 (4), 2142–2157. Krause, M. R.; Regen, S. L. The Structural Role of Cholesterol in Cell Membranes: From Condensed Bilayers to Lipid Rafts. Acc. Chem. Res. 2014, 47 (12), 3512– 3521. Nelson, D. R.; Halperin, B. I. Dislocation-Mediated Melting in Two Dimensions. Phys. Rev. B 1979, 19 (5), 2457–2484. Park, S.; Im, W. Quantitative Characterization of Cholesterol Partitioning between Binary Bilayers. J. Chem. Theory Comput. 2018, 14 (6), 2829–2833. Pandit, S. A.; Vasudevan, S.; Chiu, S. W.; Mashl, R. J.; Jakobsson, E.; Scott, H. L. Sphingomyelin-Cholesterol Domains in Phospholipid Membranes: Atomistic Simulation. Biophys. J. 2004, 87 (2), 1092–1100.
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
TOC GRAPHICS
ACS Paragon Plus Environment
Page 22 of 22