Computing Curvature Sensitivity of Biomolecules in Membranes by

Subscriber access provided by University of Leicester

Article

Computing curvature sensitivity of biomolecules in membranes by simulated buckling Federico Elias-Wolff, Martin Lindén, Alexander P. Lyubartsev, and Erik G Brandt J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00878 • Publication Date (Web): 19 Jan 2018 Downloaded from http://pubs.acs.org on January 22, 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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 47 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

Computing curvature sensitivity of biomolecules in membranes by simulated buckling Federico Elías-Wolff,†,‡ Martin Lindén,¶ Alexander P. Lyubartsev,‡ and Erik G. Brandt∗,‡ †Department of Biochemistry and Biophysics, Stockholm University, Stockholm, Sweden ‡Department of Materials and Environmental Chemistry, Stockholm University, Stockholm, Sweden ¶Department of Cell and Molecular Biology, Uppsala University, Uppsala, Sweden E-mail: [email protected] Abstract Membrane curvature sensing, where the binding free energies of membrane-associated molecules depend on the local membrane curvature, is a key factor to modulate and maintain the shape and organization of cell membranes. However, the microscopic mechanisms are not well understood, partly due to absence of efficient simulation methods. Here, we describe a method to compute the curvature dependence of the binding free energy of a membrane-associated probe molecule that interacts with a buckled membrane, which has been created by lateral compression of a flat bilayer patch. This buckling approach samples a wide range of curvatures in a single simulation, and anisotropic effects can be extracted from the orientation statistics. We develop an efficient and robust algorithm to extract the motion of the probe along the buckled membrane surface, and evaluate its numerical properties by extensive sampling

1

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

of three coarse-grained model systems: local lipid density in a curved environment for single-component bilayers, curvature preferences of individual lipids in two-component membranes, and curvature sensing by a homo-trimeric transmembrane protein. The method can be used to complement experimental data from curvature partition assays, and provides additional insight into mesoscopic theories and molecular mechanisms for curvature sensing.

1

Introduction

How cells generate and maintain their shape and spatial organization, and how the cellular shape determines function, are fundamental questions in biology. 1 Membrane curvature has increasingly been recognized as an important factor for cellular organization and membrane protein function. 1–6 It was recently shown that the curvature sensitivity of E. coli chemoreceptor clusters is strong enough to entirely drive their localization to the cell poles. 7 The bacterial proteins SpoVM and DivIVA localize to organelle surfaces or to division septa in a manner directly mediated by membrane shape. 8 Other examples of proteins that are known to both sense and generate curvature include membrane associated proteins such as BAR domains, 9–13 amphipathic and antimicrobial peptides, 14–18 and transmembrane proteins. 19–22 The orientation of curvature sensing proteins can also have important biological consequences. For example, the orientation of the bacterial actin homologue MreB is strongly correlated to cell diameter in E. coli. 23 Interactions between individual lipids and membrane curvature are weak due to the small area footprint of individual lipids, 24 but can generate membrane curvature when the monolayers are very asymmetrical. 4 However, lipids show significant curvature sensing and generating effects through cooperative effects by clustering, 25,26 and/or preferential association with proteins. 27 Apart from its fundamental scientific interest, understanding the molecular mechanisms behind curvature sensing is important for medical applications and drug design. Many diseases are associated with curvature sensing proteins like Endophilin-B1 (a BAR-domain 2


Page 2 of 47

Page 3 of 47 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


protein involved in autophagy) whose malfunction is implicated in neurodegenerative, cardiovascular and neoplastic diseases. 28 Furthermore, synthetic biology offers novel avenues for curvature sensing applications. Membrane protein-mimics constructed with DNA nanostructures (and added hydrophobic molecules so the DNA constructs insert into the membrane) can perform biological functions. 29 The 3D shapes of such mimics are easier to design than the corresponding polypeptides. Engineered curvature sensing could play an important role in the future design of membrane devices. Experimental studies of curvature sensing by lipids and proteins typically rely on fluorescent methods to monitor how putative curvature sensors partition among regions of different curvatures. 24,30–32 These are performed on membrane structures that present a range of different curvatures, such as spherical vesicles with varying sizes, 15 membrane tubes pulled from giant vesicles, 2,10 or non-flat supported bilayers. 11 For understanding curvature sensing, these approaches have limitations, since they provide no details beyond the position or local concentration of fluorescent labels. A molecular understanding of the contributions from different molecular sensing mechanisms requires more detailed information. In addition, curvature sensing by asymmetric molecules such as proteins also influences the molecular orientation, not only its position. Models and data restricted to position or density are therefore fundamentally incomplete. 33 Molecular dynamics simulations could provide an important complement to understand curvature sensing mechanisms by providing structural information with atomistic resolution, including molecular orientations, and also avoid possible complications from dye/stabilizer molecules, and interactions with bilayer supports. 34 However, the standard approach of simulating single membrane proteins that interact with a flat membrane patch is less helpful for this purpose, since curvature sensing is only indirectly visible through induced perturbations of the membrane shape, which are difficult to observe and interpret. A promising way to circumvent this difficulty is to mimic the experimental setups above and enforce curved membrane shapes that the curvature probe can interact with. One way to do this is to

3



compress a rectangular membrane patch in the longitudinal direction, which produces a stable standing wave pattern, and track the positions of the curvature sensor along the curved surface. 33,35,36 Simulated buckling has previously been used to study mechanical properties of membranes. 37–41 Compared to simulating tubes or vesicles, 42,43 the buckling approach offers a continuous distribution of curvatures, symmetry between the two monolayers of the membrane, and a single solvent compartment, eliminating the need for long equilibrations between solvent molecules in each compartment and the lipids. On the other hand, in-plane equilibration of curvature-sensitive lipid species can be computationally demanding. For curvature sensing applications, the analysis presents a particular challenge since the motion of the curvature sensors consists of two components: motion along the curved surface, and large-scale movements of the standing wave pattern itself. 33 Here, we present a systematic study of the computational aspects of this analysis using coarse-grained models of lipid bilayers and trans-membrane proteins. 44 It relies on fitting each simulated frame to a theoretical buckled shape, in order to uncouple the two motion components mentioned above. We start by simulating and analyzing a small and simple system where we can achieve perfect sampling to establish a gold standard for curvature sensing simulations. We then proceed by applying the method to two-component lipid membranes, and a physically interesting transmembrane trimer model. In section 2 we describe how to simulate a buckled membrane, and how to extract positional and orientational data from such simulations. Section 3 contains our results. We begin by evaluating the performance and robustness of the method. We then present the statistical analysis of three distinct simulated systems, and discuss our results in the framework of the biochemical literature. Finally, section 4 contains our conclusions.

4


Page 4 of 47

Page 5 of 47 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


2

Methods

In this section we describe simulations of buckled membranes with and without probe molecules. First, we detail how to setup buckled bilayer simulations. Second, we present the elastic theory of buckled membranes. Finally, we describe how to analyze the motions of probe molecules. In the following, simulation parameters are given in Lennard-Jones units, with σ as the length unit, and ε as the energy unit.

2.1

Simulating buckled membranes

All molecular dynamics simulations were performed in GROMACS 5.1 at constant volume, using the leap-frog stochastic dynamics integrator. 45 The frame alignment procedure described below, as well as the methods to detect the position and orientation of probes, are implemented as a module to the software package PLUMED, 46 written in C++ and using the dlib library for linear algebra and non-linear optimizations. 47 PLUMED was chosen as it provides enhanced sampling functionality, which is a natural extension to our method that we plan to implement in the future. The statistical analysis was performed in MATLAB, 48 and VMD was used for graphical visualization of simulated trajectories. 49 Throughout this paper, simulations are performed using the solvent-free Cooke model which represents lipids with three beads. 44 To stabilize the bilayer without solvent particles, the effect of the solvent is replaced by a long-range attractive potential between tail beads, resulting in a qualitatively correct hydrophobic balance. The model adequately reproduces bilayer self-assembly, phospholipid phase diagrams and elastic properties. 50 It also introduces some peculiarities, including an artificially high lipid flip-flop rate and, in particular, a tendency for lipids to detach from the membrane and go into gas phase before returning to the bilayer. The fraction of lipids residing in gas phase is typically around 1%. 44 The Cooke model uses units σ, ε, and τ for length, energy and time, respectively. The

5



Page 6 of 47

sizes of the lipid beads are determined by the parameter b in a Weeks-Chandler-Andersen potential 51 VWCA (r) =

 [ ]   4ε (b/r)12 − (b/r)6 + 1 , r ≤ 21/6 b 4   0 ,

r>2

1/6

(1)

b.

The values are set to bhead,head = bhead,tail = 0.95σ and btail,tail = σ, which gives the lipids an effective cylindrical shape. 50 The three lipid beads are linked by two finite extensible nonlinear elastic (FENE) bonds with stiffness 30ε/σ 2 and divergence length 1.5σ. Lipids are kept straight with a harmonic potential between the head bead and the second tail bead with rest length 4σ and bending stiffness 10ε/σ 2 . We constructed a model protein consisting of 216 Cooke beads with the same size as the lipid beads. It is formed by 162 tail beads layered between two layers of 27 head beads each. As shown in Fig. 1B, the protein is a trimer, constructed from rotated copies of an asymmetric monomer with a three-fold rotational symmetry around the membrane normal. The monomers have a wedge-like shape, with curvature radii 5σ and 25σ in the transverse and longitudinal direction, which induces a preference for non-flat membranes. The trimer structure in Fig. 1 is stabilized by an elastic network of harmonic bonds of strength 100ε/σ 2 between all pairs of particles within a distance of 4σ. The long-range attractive potential between tail beads, intended to mimic hydrophobicity, has a decay range parameter set to wc = 1.6σ. In all simulations, the temperature is set to kB T = 1.08ε. This combination of wc and T ensures a stable fluid bilayer phase. 50 The time step is δt = 0.01τ , and the friction constant is Γ = τ −1 , which defines the time scale. All beads are electrically neutral. Based on typical values of bilayer thickness and lipid diffusion, the length and time scales are approximately given by σ ≈ 1 nm, τ ≈ 10 ns. 44 A buckled membrane shape is obtained by compressing a flat bilayer patch along its longest axis (Fig. 1). We denote the length of the uncompressed bilayer (or equivalently the arc length of the bilayer’s midplane) as L, and the projected length across the α-direction as Lα . We define the compression factor as γ = 1 − Lα /L, so that γ = 0 denotes a flat bilayer. 6


Page 7 of 47 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


To simulate a buckled lipid bilayer we proceed as follows: We first assembled and equilibrated a flat rectangular bilayer patch with a Lx /Ly ratio of 2. The system is then compressed in the x-direction by a factor 1 − γ, scaling all x-coordinates and the simulation box. For high γ values (≳ 0.3) this can cause collisions between atoms, in which cases the problematic atoms are manually adjusted (for a discussion of different methods to buckle a membrane computationally see Ref. 37). Following compression, we performed energy minimization and equilibration. First, we simulated a pure single-component system with 1024 lipids, compressed to four different γ values (γ =0.2, 0.3, 0.4 and 0.5). Each simulation was run for 107 τ . The purpose of these systems is to evaluate the alignment procedure and study curvature dependent properties of the bilayer itself. Second, we simulated two-component lipid mixtures, each consisting of 512 small-headed, and 512 big-headed lipids, with different curvature preferences. 52 These lipids are obtained by modifying the value of bhead , while keeping the tail beads size constant (btail = σ). A first system, with head bead sizes bhead = 0.95σ and bhead = 1.05σ (denoted 0.95:1.05), was simulated for two values of γ = 0.3 and 0.5. The second system, with bhead = 0.9σ and bhead = 1.1σ (0.9:1.1), was simulated for γ = 0.3. In all three cases simulations were run for 6 × 106 τ . Third, we simulated a larger patch with 5966 lipids and one added transmembrane protein model (described above) acting as the curvature sensing probe molecule (Fig. 1). The bilayertrimer system was assembled with custom Matlab scripts to build a flat bilayer, insert the trimer into the bilayer, and remove overlapping lipids. This flat bilayer was energy minimized and subjected to equilibration with pressure control in the long (x) direction only, from which the average rest length in the x-direction was estimated. After equilibration, the system was compressed by 20% in the x-direction compared to the estimated rest length, again energy minimized, and then subjected to a warmup run of length 105 τ at constant box size, during which the buckled membrane shape formed and equilibrated. One thousand replicas with

7



Page 8 of 47

different initial positions for the probe were run, each for time 106 τ . The initial time 5×105 τ in each simulation were discarded as equilibration. Each replica required about 800 corehours on Intel Haswell processors.

2.2

The elastic theory of buckled membranes

The shape of a compressed lipid bilayer is the two-dimensional analogue to a compressed elastic rod as described by the Euler elastica. 53 To provide a mathematical description, we follow Ref. 37. The energy of the lipid membrane is expressed as a Helfrich-type functional of the surface shape S, 54 (

) 1 2 dA κ(K(S) − K0 ) + κ F [S] = ˜ KG (S) , 2 S ∫

(2)

where κ is the mean curvature modulus, κ ˜ is the Gaussian curvature modulus, K denotes the total curvature (twice the mean curvature), K0 the spontaneous curvature, and KG the Gaussian curvature. The contribution from Gaussian curvature is constant due to the GaussBonnet theorem for periodic surfaces with constant topology, 55 but the ground state is flat in the y direction, so KG = 0. The total curvature is given by the derivative of the tangent angle w.r.t. the arc length along the curve. We introduce the arc length parameter 0 ≤ s ≤ 1 normalized by the length L, so that K = ψ ′ (s)/L, and Ly κ F [ψ(s)] = 2L

∫

1

ds(ψ ′ (s) − LK0 )2 .

(3)

0

Minimizing F with the projected length Lx = L

∫1 0

ds cos(ψ(s)) fixed, leads to

∂ 2 ψ(s) = −λL sin(ψ(s)) , ∂s2

(4)

where λL is a Lagrange multiplier. All material constants have either dropped out or been absorbed into the Lagrange multiplier. Hence, the buckled shape is determined by the ge8


Page 9 of 47 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 1: A: Schematic representation of the buckled bilayer with a bound probe molecule (green arrow). The position s of the curvature sensing probe is defined by the projection its center of mass to the midplane of the buckled bilayer (blue surface). The orientation of the peptide is described by the angle θ between the projection of the peptide to the tangential plane and the tangent vector t. The angle between the x axis and the local tangential vector t is denoted by ψ. B: Representations of a Cooke lipid, of one of the monomers forming the curvature sensing trimer, consisting of 54 tail beads between two layers of 9 head beads each, and of the trimer (in green in panel C) from the top and lateral perspectives. C: Snapshot of simulation containing a buckled bilayer with 5966 lipids and one trimer molecule, in a configuration with ψ < 0.

9



Page 10 of 47

ometric constraints set by L and Lx . Equation (4) can be solved exactly with periodic boundary conditions. The solution is given in terms of Jacobi elliptic functions. 37 By integrating the cosine and sine of ψ(s) we can obtain the Cartesian coordinates of the buckled shape X(s), Z(s) in terms of special Jacobi functions. To avoid special functions and get a numerically tractable representation of the buckled shape, we proceed via an alternative route. We solve numerically for ψ(s) and the buckled shape, ∫

∫

s

s

dσ cos(ψ(σ)), and Z(s) =

X(s) = 0

(5)

dσ sin(ψ(σ)) , 0

with the boundary conditions

ψ(0) = ψ(L) = 0,

X(0) = Z(0) = Z(L) = 0,

(6)

X(L) = Lx .

The buckle shape is symmetric: Z(s) has period 1, is even around s = 0.5 and odd around s = 0.25, while X(s) − sLx has period 0.5, is even around s = 0, and odd around s = 0.25. This suggests that Z(s) and X(s)−sLx are well approximated by a Fourier series, with many coefficients being zero due to symmetry. For our application, where the simulation box size is constant, it convenient to use Lx as the unit length. We therefore choose an approximate representation of the buckled shape parameterized by Lx and the compression factor γ, [ XM (s, γ) = Lx s +

M ∑

] a(x) n (γ) sin(4πns) ,

n=1

[ (z)

ZM (s, γ) = Lx a0 (γ) +

M ∑

] ) a(z) . n (γ) cos 2π(2n − 1)s (

(7)

n=1

(x)

(z)

The coefficients an (γ), an (γ) are obtained by a least squares fit to the numerical solution (x)

(z)

Eq. (5). For fast evaluation of X(s) and Z(s), we construct lookup tables for an (γ), an (γ) in the range 0 ≤ γ ≤ 0.85 (at higher values of γ the buckled bilayer intersects itself), and use spline interpolation to compute the membrane shape for arbitrary γ. The splines

10


Page 11 of 47 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


can be analytically differentiated, thus derivatives w.r.t. both s and γ can be computed quickly, speeding up the alignment procedure described in the next section. The constant (z)

term a0 does not influence the shape and can be dropped with this parameterization; this corresponds to the constraints Z(0.25) = Z(0.75) = 0. In general M = 7 Fourier components are required for numerical precision of at least 5 significant digits, while M = 3 is a sufficient approximation at moderate compressions (see section S.1 in the SI.)

2.3

Analysis of simulated buckled membranes

The phase and arc length of the buckled bilayer are subject to strong fluctuations. In simulations it is difficult to filter out this motion to collect position data of the probe molecule relative to the buckled profile. We have therefore developed an alignment procedure which relies on fitting the bilayer mid-plane of each simulated frame to a reference shape which corresponds to the theoretical buckled shape of a compressed membrane. We have employed this alignment method in Refs. 33,35. Here we describe this method in detail. In order to align each frame we fit the curve in Eqs. (7), including offsets in the x and z directions, to the positions of the innermost lipid tail beads. The optimization of this fit is performed by minimizing the sum of square distances from the innermost lipid tail beads (xj , zj ) to the buckled shape X(s, γ), Z(s, γ): 1∑ (x0 + X(sj ; γ) − xj )2 + (z0 + Z(sj ; γ) − zj )2 . min χ = min x0 ,z0 ,γ,sj x0 ,z0 ,γ,sj 2 j=1 N

2

(8)

The minimization parameters are the offsets x0 and z0 , the compression factor γ (to prevent area fluctuations from biasing the fit), and the sj -values (j = 1, . . . , N where N is the number of innermost tail beads), which are the projected positions of the innermost tail-beads in terms of the arc length-like parameter s, where 0 ≤ s ≤ 1 for one period. The number of fit parameters can be reduced by using only a subset of the sj , thus decreasing N . This decreases the resolution of the fit, which we define as N divided by the total number of

11



Page 12 of 47

lipids. The differences in Eq. (8) must adhere to the minimum image convention because of the periodic boundaries. The fit is strongly nonlinear, necessitating good initial conditions. We generate them by first performing a preliminary fit in the least-squares sense to the function zj = z1 + b sin(2π(xj − x0 )/Lx − π/2). The second step is to perform the fit in Eq. (8) with a fixed γ-value (the actual value used when the simulation box was compressed), with initial sj = (xj − x0 )/Lx and the values of x0 and z0 = z1 − |b| obtained in the previous step (if b is negative then x0 = mod(x0 + Lx /2, Lx ), where mod is the modulo operation). The final step is the actual fit in Eq. (8) with γ as a free parameter (the difference between these two γ-values is around 1–2%). The first step is only necessary when (x0 , z0 ) are not available from a previous fit. Occasionally the fit will fail, particularly when aligned frames are far apart in time and the fit is performed in environments with a lot of noise, or the number of beads used for the fit is too small. In such cases a new fit will be attempted, starting with the preliminary fit and random values for x0 and z0 . Further attempts are made until the fit is successful. Once the theoretical surface shape has been obtained (i.e., after fitting for γ, x0 , and z0 ), another point of interest is to find the instant position sp ∈ (0, 1) of a membrane-embedded probe molecule. This position is given by the probe center-of-mass coordinates, (¯ xp , z¯p ), projected on the buckled surface. sp is the s-value which minimizes the squared-distance to the theoretical surface:

sp = arg min(x0 + X(s, γ) − x¯p )2 + (z0 + Z(s, γ) − z¯p )2 .

(9)

s

This nonlinear equation on the domain 0 ≤ s ≤ 1 is solved by the standard golden-section search method with the fitted parameters obtained from Eq. (8).

12


Page 13 of 47 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


The local curvature of the buckled shape is computed using the standard formula

K(s) =

Z ′ (s)X ′′ (s) − X ′ (s)Z ′′ (s) . (X ′ (s)2 + Z ′ (s)2 )3/2

(10)

Equation (10) determines the curvature sign as positive if the membrane bends away from the probe, that is, for a probe in the upper monolayer, curvature is positive if the bilayer is locally concave. Finally, to compute the orientation angle θ of a probe (Fig. 1), we consider the probe to be rigid, and compute the optimal rotation, in the least-squares sense, for a reference configuration (where θ = 0 is defined as the orientation shown in the top view of Fig. 1B) to fit the simulated probe. This procedure is known as the orthogonal Procrustes problem, 56,57 and we refer the numerical details to Section S.2 in the SI.

3

Results and discussion

In this section we evaluate the method described above. We discuss the performance of the method in terms of the number of lipids used in the analysis and robustness to noise. Then we discuss the results obtained from the single-component, the two-component mixture, and the trimer simulations described in section 2.1.

3.1

Evaluation of method performance

The method relies heavily on the alignment of each frame in the simulated trajectory, so this process must be robust to noise, as the system is subject to strong thermal fluctuations. Additionally, because of the simplified Cooke model, a small fraction of the lipids are found in a coexisting gas phase around the membrane, potentially introducing a significant amount of noise. We therefore tested the robustness to noise by adding a given number of beads in uniformly distributed random positions across the simulation box for each frame, before performing the frame alignment. 13



Figure 2 illustrates the results of the fitting procedure. The lower panel gives the rootmean-squared deviation RMSD = (χ2 /N )1/2 for each fitted frame in a 106 τ -long trajectory. χ2 is defined in Eq. (8), and N is the number of lipid beads used in the fit, excluding noise beads (similarly, though the presence of noise beads affects the value of χ2 they are not used in the evaluation of Eq. (8)). Two trajectories are shown, one with no added noise, and a second where 24 noise beads have been artificially added to each frame. The upper panels are snapshots of four individual fits, representative of the best and worst fits in both trajectories. The first panel illustrates the fitted frame with the lowest RMSD value of around 0.6. In the second panel we can observe a single lipid gone into the gas phase, contributing to a poorer fit, with an RMSD value close to 1. Still, from this second panel it is clear that the fit is good. The third and fourth panels correspond to fits in noisy environments. We observe that the worst fits, like that depicted in the fourth panel of Fig. 2 with RMSD value slightly above 1, result from the noise beads being distributed inhomogeneously. In this case there are more than twice as many noise beads below the bilayer than above it, thus the value z0 is slightly lower than ideal, yet overall the resulting fit is very reasonable. Figure 3 shows the accuracy of the frame alignment for a range of noise and resolution levels. In addition to evaluating a fit based on its RMSD, we also show the results in terms of the average deviation from the projected positions of the tail beads. We define serror = |s−s∗ | where s∗ corresponds to the zero-noise, 100% resolution simulation, which is the baseline. Panel A in Fig. 3 indicates that for a noise level up to 9.4% we can expect the RMSD to be below 0.9σ, and thus very good fits as shown in Fig. 2, with 98% confidence levels. Even the worst outliers have RMSD-values below 1.2σ, which still represent reasonable fits. In terms of serror we observe (Fig. 3C) that with the same level of confidence the worst fitted frames show an average deviation below 0.015, while the worst outliers present 5% errors in s. Noise levels higher than 10% can introduce larger deviations while also affecting performance rather significantly (Fig. S.2). The computational performance of the method depends cubically on the number of lipids

14


Page 14 of 47

Page 15 of 47 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


used in the alignment. Thus, up to a point, lowering the resolution of the fit can significantly increase time efficiency at the cost of fit accuracy, as shown in detail in Section S.3 in the SI. Panels B and D in Fig. 3 show the same indicators as above, but as functions of the fit resolution. Fig. 3B shows that the median RMSD is almost constant as resolution is decreased, though the proportions of poorer individual fits increase. When simulating a 1024lipid system and using only 50% of the innermost tail beads, almost the same distribution of RMSD values is obtained compared to 100% resolution. We find that lowering the resolution to about 25% still gives acceptable results, with the worst outliers having RMSD values below 1.2σ. In terms of serror at 25 percent level of resolution we observe that the error in the worst outliers lies around 1%. For lipid models that are less coarse-grained than the Cooke model, or even all-atom models, a scheme for defining the equivalent to the position of the inner-most tail bead must be implemented. A convenient way to treat such models would be to randomly select a subset of lipids, and define their tail position as the center of mass of their innermost carbon atoms. In conclusion, we find the minimum acceptable resolution in our one lipid component simulations to be 25%, or 256 lipids. Going any lower introduces larger deviations (Fig. 3B,D) without improving time performance (Fig. S.3). The performance penalty is due to the fact that for too low a resolution, fits will fail with high frequency (Fig. S.2), and repeated attempts (with different initial conditions) impact negatively on overall performance. We note that when fits are performed more frequently (i.e., the time lapse between fits is smaller) the probability of fitting failure decreases because of the improved initial conditions, obtained from the previous fit. Neglecting thermal fluctuations, the buckle shape depends only on the dimensionless compression parameter γ. Thus, the number of lipids needed to fit the shape is independent of system size. As the system size increases, the molecular dynamics dominates the use of computer resources over the fitting procedure. Could thermal fluctuations change the conclusion about a size-independent minimum number of fitting points? We argue that this is not the case. For a flat patch, the Helfrich

15



Page 16 of 47

model predicts that the variance of thermal undulations in a lipid membrane is 58

⟨u2 ⟩ =

AkB T , 8.3π 3 kc

(11)

where u is the membrane height deviation from the flat surface, A is the membrane patch area, and κc the bending modulus with a typical value of 20kB T . 59 Though Eq. (11) was derived for flat membranes we expect the result to be a reasonable approximation for curved membranes on short length scales. Both Rmin and ⟨u2 ⟩1/2 scale with A1/2 . For kc = 20kB T , a fixed aspect ratio L = 2Ly and γ ≤ 0.5, the value of Rmin is larger than ⟨u2 ⟩1/2 by at least one order of magnitude. Thus, regardless of the membrane patch size, thermal fluctuations enter at a higher curvature scale and are effectively averaged out during the fit to the membrane midplane. By the above arguments, the fitting procedure can be applied to buckled membranes with the patch sizes that are routinely used in molecular simulations. A possible complication for large systems is that the system may be found in higher buckling modes, especially if the aspect ratio Ly /Lx is small.

3.2

Single-component bilayer under buckling

Here we establish a gold standard in terms of sampling by simulating a single-component bilayer system, which can be run for as long as needed. The main question addressed here is how long simulation times are required to adequately sample the curvature probe at all s-values. We also determine the lipid density distribution of a membrane (for Cooke lipids) along the buckled surface. We consider simulations of a system of 1024 Cooke lipids, compressed to four increasing values of γ, and treat the lipids themselves as curvature probes. The s-value for each lipid is taken directly from the frame alignment procedure. Since the lipids are constrained to maintain a continuous bilayer, a simple interpretation in terms of curvature dependent free energy is less useful in this case.

16


Page 17 of 47 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: The RMSD of the fit to Eq. (8) as a function of time for a 1024-lipid simulation with γ = 0.2, and N = 256, corresponding to resolution of 25%. Colors correspond to noise levels: blue (no noise added) and red (24 noise beads added). The snapshots correspond to the best and worst fits for each of the noise levels; here the black curve corresponds to the fitted curve, while the colored circles represent the beads used for the fit. The simulation time is 106 τ .

17



Figure 3: A: Box plot of the fitted RMSD as function of noise level (percentage of added noise beads with respect to the number N = 256 of beads used for the fit). B: RMSD as a function of resolution level indicated by the percentage of innermost tail beads used for the fit. C–D: the same as upper panels for serror . The simulations consist of 1024 lipids, the buckling factor is γ = 0.2, simulation time is 106 τ , each box was calculated from 2 × 104 points. The whisker length is 1.5 times the interquartile range.

18


Page 18 of 47

Page 19 of 47 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


The convergence time for the curvature sensing statistics is determined by the slowest motion of interest, which is lipid lateral diffusion. To measure convergence times, we calculate the bilayer midplane density and consider the ergodic assumption that the time average of n probes, simulated for a time tsim equals the time average of one probe, simulated for a time ntsim . The midplane density is approximately uniform as the monolayers are symmetric. The computed distribution, shown in Fig. 4A, shows small systematic deviations, with maxima at s = 0.25 and s = 0.75, where the bilayer is flat. These higher order deviations from constant density are perhaps due to a difference between the bilayer midplane and the neutral plane, 40 and might be specific for the coarse-grained model. Thus, while we do not place much physical significance to this quantity, it is useful as a slowly converging measure, and an example of the ability of our method to reveal such small features. Fig. 4A shows the converged distribution, ρ∞ , in gray bars, which takes the positions of all 1024 innermost tail beads in a simulation of tsim = 107 τ . We observe that after this time the distribution of a single probe (blue line), is far from converged. Only with 100 times as many simulated time steps (orange line) does the distribution agree with the converged distribution. 2 2 Figure 4B shows the root mean squared error, RMSE = (⟨(ρi − ρ∞ i ) ⟩) , of the midplane 1

distribution as a function of the number Nlip of lipid probes (here i = 1 . . . M denotes the index of the histogram bin). As the distribution is approximately uniform, all bins ρi are equivalent and the occupancy of each is binomially distributed with probability 1/M . Accordingly, the RMSE scales as ( RMSE =

1( 1) 1 1− M M neff

) 12 ,

(12)

with the effective number of uncorrelated samples being neff = Nlip (tsim /tac ). In Section S.5 of the SI we show that the autocorrelation time for the lipids in the single-component system is tac ≈ 1090τ . The relation between RMSE and neff calculated from our simulation data is in excellent agreement with Eq. (12). Consider the case Nlip = 100, where the distribution

19



is close to convergence. For a simulation with a single probe, this would roughly correspond to a simulation time of 10 seconds (using the Cooke model mapping of units discussed in section 2.1), vastly exceeding current capabilities for all but the most coarse-grained models. For all-atom MD simulations, the longest published simulations (to our knowledge) of small proteins in explicit water are in the 1 ms range. To access such timescales in physically realistic systems requires both the use of a coarse grained model like Martini (which can provide about 3 orders of magnitude faster phase space sampling), 60 and enhanced sampling. The application of a metadynamics scheme is a natural extension to our method, using the values of s and θ from the probe molecule to coarse-grain the free energy landscape.

Figure 4: A: Probability density ρ(s) of the lipids arc length position s, projected to the bilayer midplane. The histogram is constructed from the s-values of the 1024 tail beads used in the fit to the buckled shape. Colored lines correspond to using the innermost tail beads of only 1,10, and 100 lipids. B: Symbols correspond to the average root mean squared error as a function of the number Nlip of lipids used to compute ρ(s). The red line corresponds to the theoretical RMSE, Eq. (12), with tac = 1090τ . Simulation time is tsim = 107 τ , the compression factor is γ = 0.2. The curvatures of the buckled bilayers in the present work are very high. The maximum curvature corresponds to a radius of 3.5 σ (roughly 3.5 nm in the Cooke model), which is comparable or even smaller than the bilayer thickness. We therefore expect the lipid density difference between monolayers in such structures to be large. In regions of positive

20


Page 20 of 47

Page 21 of 47 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


curvature lipids have an increased surface area available, and as Fig. 5B shows, they are more concentrated there. This effect is of the order of γ. The different degrees of buckling in the simulations imply that a wide range of curvatures is sampled, as shown in Fig. 5A. The density is to an excellent approximation a function only of curvature for all simulations (Fig. 5B.) The simplest geometric model is that the local density is linear in the curvature,

ρ(s) = ρ0 + a1 K(s) .

(13)

Fig. 5B shows that the linear model is an excellent approximation. Since all the lipids have the same curvature preferences, the density in the buckled membrane is mostly a geometric effect, reflecting the increased area available in the monolayer under positive curvature. However, higher order deviations from this model are discernible. Fig. 5C depicts the relationship between midplane density and curvature. The linear model implies that the upper and lower monolayers balance out, yielding a flat density (dashed black line), but the coarse-grained model shows a local density maximum at K = 0, minima near |K| ∼ 0.2, and increasing density at higher curvatures which are only accessible for the most strongly buckled systems. Given that the model is highly coarse-grained, and the possibility that small deviations from harmonic elasticity could potentially bias the fit, it is not clear how relevant this observation is to real lipids, but it does illustrate the potential of the buckling analysis to reveal new and interesting phenomena.

3.3

Two-component bilayer

To validate our method quantitatively, we study curvature dependent lipid sorting in a binary mixture. Following Ref. 52, we construct 1:1 mixtures of lipids with different head bead sizes, and therefore different curvature preferences. The two monolayers of a curved bilayer are subject to opposite mean curvatures, which leads to different compositions. In

21



Figure 5: A: Curvature K(s) for γ = 0.2, 0.3, 0.4 and 0.5, for the single-component system, with L = 35 σ. The curves correspond to the curvature of the bilayer midplane, from the perspective of the upper leaflet, where curvature is defined as positive if the membrane bends away from the probe. B: Curvature dependent probability density of lipid tails from one monolayer for different values of γ. The dashed line is the linear model ρ(s) = ρ0 + a1 K(s), where ρ0 = 1 (a flat monolayer) and a1 is obtained from a least-squares fit. C: Taking ρ(s) + ρ(s + 0.5) corresponds to taking the midplane density including both monolayers [ρ(s) = ρupper (s), ρ(s + 0.5) = ρlower (s) and K(s) = −K(s + 0.5)]. Ref. 52, lipid sorting in a range of curvatures is investigated by simulating spherical vesicles of varying sizes. The result is shown as the relationship between total curvature K and the so-called enhancement ratio, defined as ln(ϕ+ (s)/ϕ− (s)), where ϕ+ , ϕ− are the molar fractions of large-headed lipids, with local positive and negative curvature, respectively (in vesicles, the curvature of the outer leaflet is positive). Using the buckling method, we can access a range of curvatures in a single run. We simulated two mixtures, the first with head bead size ratio 0.95:1.05, and the second with 0.9:1.1, in a range of curvatures similar to those of Ref. 52. As shown in Fig. 6, larger size differences produce stronger sorting using both methods, and different amounts of buckling produce consistent results, just as for the one-component case. Moreover, the results for the 1.1:0.9 mixture show excellent quantitative agreement, except for two points at high curvature. These may be outliers, since they deviate from the predicted linear relationship 52 shown by the rest of the data. In theory, they could reflect sorting due to Gaussian curvature, which is zero in the buckled system but 1/R2 = (K/2)2 on a spherical vesicle surface with radius R. However, this option is unlikely given that the 22


Page 22 of 47

Page 23 of 47 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


monolayer Gaussian modulus of the lipids depend only weakly on head bead size. 61 The two lipid species are therefore expected to respond similarly to Gaussian curvature, and produce a weak sorting effect. On the other hand, the total curvature is coupled to the lipid spontaneous curvature, which does vary strongly with head bead size. 52,61 Also, the two monolayers have equal Gaussian curvatures when a single surface (such as the bilayer mid-plane) is used as reference. Hence, detecting sorting by Gaussian curvature between opposing monolayers requires a theory that accounts for internal bilayer structure. An indepth analysis is needed to sort out these interesting questions, but this is beyond the scope of the present work.

Figure 6: Enhancement ratio as function of absolute total curvature |K| (twice the mean curvature), for different head-size ratios. The data marked C&D in the legend was taken from Fig. 4 in Ref. 52.

3.4

Trimer-bilayer system

We now turn to simulating a model trimer protein embedded in a buckled bilayer, as described in section 2.1. In addition to tracking the position of the trimer across the buckled membrane, its orientation also becomes relevant. Some of the largest membrane curvatures in cells are found in the tubules of endoplasmic reticulum in Saccharomyces cerevisiae, with curvature radii of 15 nm. 4 This is comparable 23



to the maximum curvature of the system presented here. We expect the trimer to localize to regions in the bilayer where the surrounding bilayer shape matches its conical cross-section, which minimizes elastic bilayer deformations. The minimal radius of curvature in the buckled membrane analyzed here is approximately 15σ. As the monomers have curvature radii 5σ and 25σ, it is reasonable to see the trimer prefer a membrane region with intermediate curvature. In addition, the higher bead density at the bottom of the trimer may lead to stronger interactions with surrounding lipids, which could also influence the curvature preference. To obtain an estimate of the simulation time required to obtain adequate sampling levels, we consider the single-component lipid system, as both systems are diffusion limited (typical values for the diffusion constant of membrane proteins are about 5–10 µm2 /s). 62 Reported times for simulations of peptides in buckled bilayers using the Martini force field are of the order of 104 τ for system sizes of 1024 lipids. 33,35 In section S.5 of the SI we show that the autocorrelation time for the trimer’s s-position in a flat bilayer is two orders of magnitude larger than for individual lipids. From Fig. 4, we would thus expect a simulation time of order 1010 τ . On the other hand, the membrane phase fluctuations can be considered a motion that speeds up the phase space sampling in terms of s for a single probe. Proper sampling in this system also requires the trimer to fully explore all six peaks in phase space (Fig. 7B) corresponding to the symmetries of the system. We proceed by simulating the system for 109 τ . We achieve this by running a thousand replicas, each simulated for 106 τ . We estimated the trimer autocorrelation time to about 2.3 × 105 τ (Fig. S.6 in the SI). Thus, the initial 5 × 105 τ of each replica simulation were discarded as equilibration. To get the initial conditions for each replica, we first ran a simulation of 2.6 × 107 τ , and obtained a time series of (s(t), θ(t)). We then created a 35 × 35 lattice (s ∈ (0, 1), θ ∈ (0, 2π)), and labeled the trajectory frames close to each lattice point (∆s < 1/100 and ∆θ < 2π/100). From the lattice points populated by more than one trajectory frame we randomly chose one to be the initial frame of each replica. From the 1225 points, 155 were unpopulated, and 70 more

24


Page 24 of 47

Page 25 of 47 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


were randomly thrown out, resulting in 1000 initial frames. To extract a binding free energy profile for the trimer in this lipid environment we model e the distribution of the probe position as depending on local curvature, so G(s) = G(K(s)), e where G(s) denotes the orientation averaged binding free energy in units of kB T . The dise

tribution density should be proportional to the Boltzmann factor e−G(s) , thus we calculate the free energy G from curvature histograms by weighting each bin with its Jacobian transformation factor: e

ρe(s)ds ∝ e−G(s) ds ∝ e−G(K) |dK/ds|−1 dK ∝ ρ(K)dK ,

(14)

where normalization constants were dropped. The curvature-dependent free energy is then ( ) G(K) = − ln ρ(K)|dK/ds| + const.

(15)

The Jacobian |dK/ds| compensates for intermediate values of curvature corresponding to shorter ranges in arc length while extremal values of curvature have longer arc length foote prints. Fig. 7D shows that the free energies G(s) and G(K), calculated from position and curvature histograms respectively, are identical. The distribution of the trimer position in Fig. 7A and the corresponding free energy profile in Fig. 7D both show the protein preferring a region of intermediate curvature, as expected from its structure. From a quadratic fit to the free energy (dashed in Fig. 7D), we estimate the preferred curvature to 0.035σ −1 . The binding free energy of single proteins is usually modeled with a Taylor expansion in the local curvature. 2 If second order terms are included, a free energy minimum is expected at some preferred curvature, as in Fig. 7D. This prediction has been borne out in a few recent experiments. 63,64 In other cases, the protein sensors have localized to the maximal available curvature, 9,15,30,65 either because the assays did not probe a wide enough range of curvatures, or due to effects beyond the phenomenological quadratic model. 25



Page 26 of 47

We now discuss the symmetries of the system with respect to the joint distribution of the trimer’s position and orientation ρ(s, θ), depicted in Fig. 7B. Given the 3-fold structural symmetry of the trimer, we expect the distribution to be symmetric under rotations by ±120 degrees, (s, θ) → (s, θ + 120◦ ). Similarly, the distribution is symmetric under a rotation of the whole system by 180 degrees around the z-axis, (s, θ) → (1 − s, θ + 180◦ ). These two symmetries can be clearly observed in Fig. 7B. However, a common approach to modeling direction-dependent curvature sensing is to take the binding energy to be a function of the local curvature tensor seen from a coordinate system moving with the protein, 33,66–70 which implies an additional local symmetry. To see this, we note that in a coordinate system rotated by θ with respect to one of the principal curvature directions, the local curvature tensor is given by 33   D sin 2θ  H + D cos 2θ Cij =  , D sin 2θ H = D cos 2θ

(16)

where H = (c1 + c2 )/2 and D = (c1 − c2 )/2 are the mean and deviatoric curvatures, and c1 and c2 the principal curvatures. Thus, if the binding energy is a function only of Cij , we would expect the Boltzmann distribution to be invariant under 180-degree rotations, (s, θ) → (s, θ + 180), but this is clearly not the case in Fig. 7B. We believe this to reflect slight differences in local curvature at the individual monomers, meaning that the trimer effectively senses a curvature gradient, which is not invariant under a 180 degree rotation.

4

Conclusions

Membrane curvature sensing is an important driving force for shaping and organizing cell membranes, but is incompletely understood both at the atomistic level of single molecules and sensing motifs, as well as the mesoscopic level of cellular structures. Molecular simulations offer an important complement to biophysical experiments, but standard approaches

26


Page 27 of 47 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 7: A: Normalized distribution of the arc length position s from the trimer simulation. B: Joint distribution ρ(s, θ) of the position and orientation of the trimer. Density is indicated by colors from blue (low probability) to red (high probability). The data includes 1000 replicas, each simulated for t = 106 τ , excluding the first 5 × 105 τ as equilibration. C: Trajectory snapshots from top and side perspectives, the s and θ values correspond ) ( to the star in panel B. D: Free energy G(K(s)) versus curvature, where G(K(s)) = − ln ρ(K(s)) dK ds e and G(s) = − ln(e ρ(s)). The red dashed line is a quadratic fit to the data shown in circles, with fitted parameters K0 = 0.035 σ −1 , b = 1005 kB T σ 2 , and G0 = −1.26 kB T .

27



rely on indirect signals like induced perturbations on flat membranes or pressure profiles which are often weak and difficult to interpret. We study the computational aspects of an alternative approach based on simulated buckling. Here, the main signal is the position and orientation of the sensor molecule interacting with a curved membrane, which can be directly related to the relevant curvature-dependent free energy. The buckling analysis requires alignment of different frames to differentiate the relevant probe motion from the large-scale motion of the buckled profile. We solve this problem by a nonlinear least squares fit to a theoretical shape using an efficient parameterization, and show that the procedure is robust against different noise sources. Compression-induced buckling also produces complex in-plane tensions in the bilayer, 37,38 and consistency is therefore a potential concern: how equivalent are local curvatures from different compression ratios? To address this, we show that lipid density of a single component bilayer, as well as the molar ratio in two-component mixtures, are to very good approximation functions of only local curvature over a range of compression ratios (Figs. 5 and 6). Our results also highlight the potential of the buckling method to reveal features not easily obtained by other methods, such as the nonlinear curvature dependence of the lipid density in the one-component system (Fig. 5B,C), and the in-plane orientation of the trimer (Fig. 7B). The main challenge for applying the method to more detailed models is the simulation time needed for the configuration distribution of the probe molecule to convergence. Sampling is limited by the lipid diffusion rate when the curvature profile is nearly flat. In a more realistic scenario, such as determining the curvature sensing of a symmetry-constrained protein, the embedded probe is forced to cross membrane regions with large curvatures. Long times are required to overcome the associated free energy barriers (Fig. 7D). As an illustrative example, we determined the curvature sensing properties of a model trimer protein by brute-force simulations of 1000 replicas, with the total simulation time reaching seconds. A more efficient way to handle fine-detailed simulation models would require enhanced sampling of regions in phase space with unfavorable curvatures. Our method is easily extended

28


Page 28 of 47

Page 29 of 47 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


to enhanced sampling algorithms such as metadynamics, which bias selected “collective” degrees of freedom to speed up calculations. 71 The position of the protein inclusion along the curved membrane, defined by us in Eq. (9), could favorably be used as a collective variable. Our method applies to a simple curved geometry, amenable to exact analysis. A buckled membrane is flat in the transverse direction, and can only access a subset of all possible curvature states. In particular, the Gaussian curvature is always zero. While Cooke lipids appear to be insensitive to Gaussian curvature (Fig. 6), this is a significant restriction, as some proteins like MreB have been shown to sense local Gaussian curvature. 23 Our method has the potential to be generalized to other curved geometries by discarding the Euler elastica parameterization. Using local fitting to quadratic shapes provides a general way to fit a variety of curved structures. 72–74 This opens up the possibility to investigate highly curved and asymmetric shapes where Gaussian curvature plays a key role, such as stalk formation and vesicle fusion. 75

Acknowledgement This work has been supported by the Swedish Research Council (Vetenskapsrådet), the EU Horizon 2020 “SmartNanoTox” project, and the Swedish Foundation for Strategic Research via the Center for Biomembrane Research. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing at the KTH Royal Institute of Technology. The authors would like to thank Jon Kapla for sharing a Matlab library to read and write Gromacs trajectories, Arnold Maliniak and an anonymous reviewer for constructive criticism of the manuscript.

Supporting Information Available The Supporting Information contains mathematical details on the Fourier representation 29



of the theoretical buckled shape, a detailed description of the solution of the orthogonal Procrustes problem to detect the orientation angle of the trimer, an analysis of the time performance of our method, additional details on the geometric effects of one-component buckled bilayers, and the autocorrelation analysis for the one-component and trimer simulations. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Zimmerberg, J.; Kozlov, M. M. How proteins produce cellular membrane curvature. Nat. Rev. Mol. Cell Bio. 2006, 7, 9–19. (2) T. Baumgart et al, Thermodynamics and Mechanics of Membrane Curvature Generation and Sensing by Proteins and Lipids. Annu. Rev. Phys. Chem. 2011, 62, 483–506. (3) McMahon, H. T.; Gallop, J. L. Membrane curvature and mechanisms of dynamic cell membrane remodelling. Nature 2005, 438, 590–596. (4) Shibata, Y.; Hu, J.; Kozlov, M. M.; Rapoport, T. A. Mechanisms Shaping the Membranes of Cellular Organelles. Annu. Rev. Cell Dev. Biol. 2009, 25, 329–354. (5) Graham, T. R.; Kozlov, M. M. Interplay of proteins and lipids in generating membrane curvature. Curr. Opin. Cell Biol. 2010, 22, 430–436. (6) Phillips, R.; Ursell, T.; Wiggins, P.; Sens, P. Emerging roles for lipids in shaping membrane-protein function. Nature 2009, 459, 379–385. (7) Draper, W.; Liphardt, J. Origins of chemoreceptor curvature sorting in Escherichia coli. Nat. Commun. 2017, 8, 14838. (8) Huang, K. C.; Ramamurthi, K. S. Macromolecules that prefer their membranes curvy. Mol. Microbiol. 2010, 76, 822–832.

30


Page 30 of 47

Page 31 of 47 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


(9) Peter, B. J.; Kent, H. M.; Mills, I. G.; Vallis, Y.; Butler, P. J. G.; Evans, P. R.; McMahon, H. T. BAR domains as sensors of membrane curvature: the amphiphysin BAR structure. Science 2004, 303, 495–499. (10) Sorre, B.; Callan-Jones, A.; Manzi, J.; Goud, B.; Prost, J.; Bassereau, P.; Roux, A. Nature of curvature coupling of amphiphysin with membranes depends on its bound density. Proc. Natl. Acad. Sci. USA 2012, 109, 173–178. (11) Hsieh, W.-T.; Hsu, C.-J.; Capraro, B. R.; Wu, T.; Chen, C.-M.; Yang, S.; Baumgart, T. Curvature Sorting of Peripheral Proteins on Solid-Supported Wavy Membranes. Langmuir 2012, 28, 12838–12843. (12) Eriksson, H. M.; Wessman, P.; Ge, C.; Edwards, K.; Wieslander, Å. Massive formation of intracellular membrane vesicles in Escherichia coli by a monotopic membrane-bound lipid glycosyltransferase. J. Biol. Chem. 2009, 284, 33904–33914. (13) Stachowiak, J. C.; Schmid, E. M.; Ryan, C. J.; Ann, H. S.; Sasaki, D. Y.; Sherman, M. B.; Geissler, P. L.; Fletcher, D. A.; Hayden, C. C. Membrane bending by protein–protein crowding. Nat. Cell Biol. 2012, 14, 944–949. (14) Drin, G.; Casella, J.-F.; Gautier, R.; Boehmer, T.; Schwartz, T. U.; Antonny, B. A general amphipathic α-helical motif for sensing membrane curvature. Nat. Struct. Mol. Biol. 2007, 14, 138–146. (15) Hatzakis, N. S.; Bhatia, V. K.; Larsen, J.; Madsen, K. L.; Bolinger, P.-Y.; Kunding, A. H.; Castillo, J.; Gether, U.; Hedegård, P.; Stamou, D. How curved membranes recruit amphipathic helices and protein anchoring motifs. Nat. Chem. Biol. 2009, 5, 835–841. (16) Bhatia, V. K.; Madsen, K. L.; Bolinger, P.-Y.; Kunding, A.; Hedegård, P.; Gether, U.; Stamou, D. Amphipathic motifs in BAR domains are essential for membrane curvature sensing. EMBO J. 2009, 28, 3303–3314. 31



(17) Mizuno, N.; Varkey, J.; Kegulian, N. C.; Hegde, B. G.; Cheng, N.; Langen, R.; Steven, A. C. Remodeling of lipid vesicles into cylindrical micelles by α-synuclein in an extended α-helical conformation. J. Biol. Chem. 2012, 287, 29301–29311. (18) Braun, A. R.; Sevcsik, E.; Chin, P.; Rhoades, E.; Tristram-Nagle, S.; Sachs, J. N. αSynuclein induces both positive mean curvature and negative Gaussian curvature in membranes. J. Am. Chem. Soc. 2012, 134, 2613–2620. (19) Tonnesen, A.; Christensen, S. M.; Tkach, V.; Stamou, D. Geometrical Membrane Curvature as an Allosteric Regulator of Membrane Protein Structure and Function. Biophys. J. 2014, 106, 201–209. (20) Schmidt, N. W.; Mishra, A.; Wang, J.; DeGrado, W. F.; Wong, G. C. Influenza virus A M2 protein generates negative Gaussian membrane curvature necessary for budding and scission. J. Am. Chem. Soc. 2013, 135, 13710–13719. (21) Arechaga, I.; Miroux, B.; Karrasch, S.; Huijbregts, R.; de Kruijff, B.; Runswick, M. J.; Walker, J. E. Characterisation of new intracellular membranes in Escherichia coli accompanying large scale over-production of the b subunit of F(1)F(o) ATP synthase. FEBS Lett. 2000, 482, 215–219. (22) Davies, K. M.; Anselmi, C.; Wittig, I.; Faraldo-Gómez, J. D.; Kühlbrandt, W. Structure of the yeast F1Fo-ATP synthase dimer and its role in shaping the mitochondrial cristae. Proc. Natl. Acad. Sci. USA 2012, 109, 13602–13607. (23) Ouzounov, N.; Nguyen, J. P.; Bratton, B. P.; Jacobowitz, D.; Gitai, Z.; Shaevitz, J. W. MreB orientation correlates with cell diameter in Escherichia coli. Biophys. J. 2016, 111, 1035–1043. (24) Kamal, M. M.; Mills, D.; Grzybek, M.; Howard, J. Measurement of the membrane curvature preference of phospholipids reveals only weak coupling between lipid shape and leaflet curvature. Proc. Natl. Acad. Sci. USA 2009, 106, 22245–22250. 32


Page 32 of 47

Page 33 of 47 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


(25) Mukhopadhyay, R.; Huang, K. C.; Wingreen, N. S. Lipid localization in bacterial cells through curvature-mediated microphase separation. Biophys. J. 2008, 95, 1034–1049. (26) Renner, L. D.; Weibel, D. B. Cardiolipin microdomains localize to negatively curved regions of Escherichia coli membranes. Proc. Natl. Acad. Sci. USA 2011, 108, 6264– 6269. (27) Sorre, B.; Callan-Jones, A.; Manneville, J.-B.; Nassoy, P.; Joanny, J.-F.; Prost, J.; Goud, B.; Bassereau, P. Curvature-driven lipid sorting needs proximity to a demixing point and is aided by proteins. Proc. Natl. Acad. Sci. USA 2009, 106, 5622–5626. (28) Frost, A.; Unger, V. M.; De Camilli, P. The BAR Domain Superfamily. Cell 2009, 137, 191–196. (29) Howorka, S. Changing of the guard. Science 2016, 352, 890–891. (30) Ramesh, P.; Baroji, Y. F.; Reihani, S. N. S.; Stamou, D.; Oddershede, L. B.; Bendix, P. M. FBAR Syndapin 1 recognizes and stabilizes highly curved tubular membranes in a concentration dependent manner. Sci. Rep. 2013, 3, 1565. (31) Zhu, C.; Das, S. L.; Baumgart, T. Nonlinear Sorting, Curvature Generation, and Crowding of Endophilin N-BAR on Tubular Membranes. Biophys. J. 2012, 102, 1837–1845. (32) Aimon, S.; Callan-Jones, A.; Berthaud, A.; Pinot, M.; Toombes, G. E.; Bassereau, P. Membrane shape modulates transmembrane protein distribution. Dev. Cell 2014, 28, 212–218. (33) Gómez-Llobregat, J.; Elías-Wolff, F.; Lindén, M. Anisotropic membrane curvature sensing by amphipathic peptides. Biophys. J. 2016, 110, 197–204. (34) Alejo, J. L.; Blanchard, S. C.; Andersen, O. S. Small-Molecule Photostabilizing Agents are Modifiers of Lipid Bilayer Properties. Biophys. J. 2013, 104, 2410–2418.

33



(35) Martyna, A.; Gomez-Llobregat, J.; Lindén, M.; Rossman, J. S. Curvature sensing by a viral scission protein. Biochemistry 2016, 55, 3493–3496. (36) Wang, H.; de Joannis, J.; Jiang, Y.; Gaulding, J. C.; Albrecht, B.; Yin, F.; Khanna, K.; Kindt, J. T. Bilayer Edge and Curvature Effects on Partitioning of Lipids by Tail Length: Atomistic Simulations. Biophys. J. 2008, 95, 2647–2657. (37) Hu, M.; Diggins IV, P.; Deserno, M. Determining the bending modulus of a lipid membrane by simulating buckling. J. Chem. Phys 2013, 138, 214110. (38) Noguchi, H. Anisotropic surface tension of buckled fluid membranes. Phys. Rev. E 2011, 83, 061919. (39) Otter, W. K. d. Area compressibility and buckling of amphiphilic bilayers in molecular dynamics simulations. J. Chem. Phys. 2005, 123, 214906. (40) Wang, X.; Deserno, M. Determining the pivotal plane of fluid lipid membranes in simulations. J. Chem. Phys. 2015, 143, 164109. (41) Wang, X.; Deserno, M. Determining the Lipid Tilt Modulus by Simulating Membrane Buckles. J. Phys. Chem. B 2016, 120, 6061–6073. (42) Harmandaris, V. A.; Deserno, M. A novel method for measuring the bending rigidity of model lipid membranes by simulating tethers. J. Chem. Phys. 2006, 125, 204905. (43) Risselada, H. J.; Marrink, S. J. Curvature effects on lipid packing and dynamics in liposomes revealed by coarse grained molecular dynamics simulations. Phys. Chem. Chem. Phys. 2009, 11, 2056–2067. (44) Cooke, I. R.; Deserno, M. Solvent-free model for self-assembling fluid bilayer membranes: stabilization of the fluid phase based on broad attractive tail potentials. J. Chem. Phys. 2005, 123, 224710.

34


Page 34 of 47

Page 35 of 47 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


(45) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19–25. (46) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604–613. (47) King, D. E. Dlib-ml: A Machine Learning Toolkit. J. Mach. Learn. Res. 2009, 10, 1755–1758. (48) MATLAB 2014a, The MathWorks Inc.: Natick, Massachusetts, United States. (49) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33–38. (50) Cooke, I.; Kremer, K.; Deserno, M. Tunable generic model for fluid bilayer membranes. Phys. Rev. E 2005, 72, 011506. (51) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237–5247. (52) Cooke, I. R.; Deserno, M. Coupling between lipid shape and membrane curvature. Biophys. J. 2006, 91, 487–495. (53) Antman, S. Nonlinear Problems of Elasticity; Applied Mathematical Sciences; Springer New York, 2005. (54) Helfrich, W. Elastic properties of lipid bilayers: theory and possible experiments. Z. Naturforsch. C 1973, 28, 693–703. (55) Pressley, A. N. Elementary differential geometry; Springer Science & Business Media, 2010.

35



(56) Brandt, E. G.; Lyubartsev, A. P. Molecular dynamics simulations of adsorption of amino acid side chain analogues and a titanium binding peptide on the TiO2 (100) surface. The Journal of Physical Chemistry C 2015, 119, 18126–18139. (57) Gower, J. C.; Dijksterhuis, G. B. Procrustes problems; Oxford University Press on Demand, 2004; Vol. 30. (58) Edholm, O. Time and length scales in lipid bilayer simulations. Curr. Top. Membr. 2008, 60, 91–110. (59) Nagle, J. F.; Jablin, M. S.; Tristram-Nagle, S.; Akabori, K. What are the true values of the bending modulus of simple lipid bilayers? Chem. Phys. Lipids 2015, 185, 3–10. (60) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812–7824. (61) Hu, M.; Briguglio, J. J.; Deserno, M. Determining the Gaussian Curvature Modulus of Lipid Membranes in Simulations. Biophys. J. 2012, 102, 1403–1410. (62) Weiß, K.; Neef, A.; Van, Q.; Kramer, S.; Gregor, I.; Enderlein, J. Quantifying the diffusion of membrane proteins and peptides in black lipid membranes with 2-focus fluorescence correlation spectroscopy. Biophys. J. 2013, 105, 455–462. (63) Prévost, C.; Zhao, H.; Manzi, J.; Lemichez, E.; Lappalainen, P.; Callan-Jones, A.; Bassereau, P. IRSp53 senses negative membrane curvature and phase separates along membrane tubules. Nat. Commun. 2015, 6, 8529. (64) Alnaas, A. A.; Moon, C. L.; Alton, M.; Reed, S. M.; Knowles, M. K. Conformational Changes in C-Reactive Protein Affect Binding to Curved Membranes in a Lipid Bilayer Model of the Apoptotic Cell Surface. J. Phys. Chem. B 2017, 121, 2631–2639.

36


Page 36 of 47

Page 37 of 47 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


(65) Bigay, J.; Casella, J.-F.; Drin, G.; Mesmin, B.; Antonny, B. ArfGAP1 responds to membrane curvature through the folding of a lipid packing sensor motif. EMBO J. 2005, 24, 2244–2253. (66) Fournier, J. Nontopological saddle-splay and curvature instabilities from anisotropic membrane inclusions. Phys. Rev. Lett. 1996, 76, 4436. (67) Perutková, Š.; Kralj-Iglič, V.; Frank, M.; Iglič, A. Mechanical stability of membrane nanotubular protrusions influenced by attachment of flexible rod-like proteins. J. Biomech. 2010, 43, 1612–1617. (68) Ramakrishnan, N.; Kumar, P. S.; Ipsen, J. H. Monte Carlo simulations of fluid vesicles with in-plane orientational ordering. Phys. Rev. E 2010, 81, 041922. (69) Akabori, K.; Santangelo, C. D. Membrane morphology induced by anisotropic proteins. Phys. Rev. E 2011, 84, 061909. (70) Ramakrishnan, N.; Sunil Kumar, P. B.; Ipsen, J. H. Membrane-Mediated Aggregation of Curvature-Inducing Nematogens and Membrane Tubulation. Biophys. J. 2013, 104, 1018–1028. (71) Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601. (72) Bolle, R. M.; Vemuri, B. C. On three-dimensional surface reconstruction methods. IEEE T. Pattern Anal. 1991, 13, 1–13. (73) Ohtake, Y.; Belyaev, A.; Alexa, M.; Turk, G.; Seidel, H.-P. Multi-level partition of unity implicits. ACM Siggraph 2005 Courses. 2005; p 173. (74) Levin, D. Geometric modeling for scientific visualization; Springer, 2004; pp 37–49.

37



(75) Siegel, D. P. The Gaussian curvature elastic energy of intermediates in membrane fusion. Biophys. J. 2008, 95, 5200–5215.

38


Page 38 of 47

Page 39 of 47 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


Graphical TOC Entry

39


Journal of Chemical Theory and Computation Page 40 of 47 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


Page 41 of 47 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



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


Page 42 of 47

Page 43 of 47 1 2 3 4 5 6 7 8 9 10 11 12



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


Page 44 of 47

Page Journal 45 ofof 47Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12


Journal of Chemical Theory and ComputationPage 46 of 47 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


Page 47 Journal of 47 of Chemical Theory and Computation 1 2 3 4


Computing Curvature Sensitivity of Biomolecules in Membranes by

Recommend Documents