Free energy surface and molecular characterization of slow structural transitions in lipid bilayers

The need to incorporate specific molecular-scale features for largescale structural changes in biological membranes necessitate use of a multi-scale computational approach. Here, this comprises of Langevin dynamics in a normal mode space determined from an elastic network model (ENM) representation for lipid-water Hamiltonian. All atom (AA) MD simulations are used to determine model parameters, and Langevin dynamics predictions for an extensive set of bilayer properties, such as, undulation spectra, undulation relaxation rates, dynamic structure factor, and mechanical properties are validated against the data from MD simulations and experiments. The transferability of model parameters to describe dynamics of a larger lipid bilayer and a heterogeneous membrane-protein system is assessed. The developed model is coupled to the energy landscape for membrane deformations to obtain a set of generic reaction coordinates (RCs) for pore formation in two tensionless, single lipid-type bilayers, namely, 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) and 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC). Structure evolution is carried in an AA MD simulation wherein the generic RCs are used in a path metadynamics or an umbrella sampling simulation to investigate thermodynamics of pore formation and its molecular determinants. The transition state is characterized extensively to bring out the interplay between various bilayer motions (undulations, lateral density fluctuations, thinning, lipid tilt), lipid solvation, and lipid packing.

brane undulations(18), membrane fusion(19, 20), and pore 23 formation(21). However, they often require strong assump-24 tions on the nature of transition intermediates (19-21) and 25 predicted energetics were not quantitative for highly curved 26 states (22). This will lead to an inaccurate free energy land-27 scape and kinetics even for processes whose end-states are de-28 scribed well by the continuum description. At a more general    for S(q) of DMPC bilayer over the entire range ( Fig. S5(C)).

206
The eigen-decomposition of H gives the eigenvectors, Um 207 (Cartesian equivalent, obtained as per Eq. S8), also re-208 ferred to as the normal modes, and the associated vibra- along these modes has the characteristic parabolic shape as 222 obtained from a 1 µs MD simulation, with corresponding ENM 223 predictions in quantitative agreement ( Fig. S6(B)). The bi-224 layer motion in next set of higher frequency modes can involve 225 a single fluctuation type at multiple wave-vectors (modes 21-226 24 capture density fluctuations at q = 0.35, and 0.50 nm −1 ) 227 and, also, a combination of height, thickness, and/or lateral 228 density fluctuations at multiple wave-vectors ( Fig. S6(E)). 229 The role of these distinct set of modes in membrane remodel-230 ing is elaborated for the case of pore formation in a tensionless 231 bilayer.

232
The friction coefficients, Υ, are determined from the ve-233 locity auto-correlation function of the lipid beads, calculated 234 from a 2 ns simulation of a 256 lipid bilayer, as per Eq. S9 235 (66). These lie in the range 7 ps −1 to 40 ps −1 , with the beads 236 at or connected to the headgroup at the upper-end of the 237 spectrum and the tail beads at the lower-end (Fig. S7(A-B)). 238 This correlates well with difference in packing along the bi-239 layer normal (Fig. S7(B)), wherein the density is lowest at the 240 bilayer center. The obtained values are close to the Stoke's law 241 estimate of Υii = 62 ps −1 ( Υii = m −1 6πηa: water viscosity 242 at 298 K, η = 0.8 mPa s, lipid CG bead hydrodynamic radius, 243 a = 0.30 nm), implying a behavior close to a Newtonian fluid. 244 Fig. 1(D) shows that the diagonal elements of the friction 245 coefficients matrix in the normal mode space, Γ = U T ΥU 246 (66), lie in a narrow range of 20 ps −1 to 25 ps −1 , while the 247 non-diagonal elements, which corresponds to inter-mode hydro-248 dynamics coupling, are much smaller (Fig. S7(D)). Therefore, 249 these non-diagonal elements are set to zero, giving 3N − 3 250 uncoupled Langevin equations (Eq. S10) in the normal modes 251 space(66, 78), referred here as the uL-NMS model. In this case, 252 most lipid bilayer properties of interest, such as undulations 253 spectra and relaxation dynamics, static and dynamic structure 254 factors, can be derived analytically (see Supplemental Material 255 Sections C -D).

256
Quantitative assessment of the Langevin model. The predic-257 tions for long timescale dynamics and equilibrium properties 258 from uL-NMS are compared with AA MD simulations and ex-259 perimental data. A 1 µs AA MD simulation of a 1024 DMPC 260 bilayer is used for property calculation, however, uL-NMS 261 parameters are taken from shorter (2 ns to 100 ns) simulations 262 of a smaller (256 lipid) bilayer. The static undulation spectra 263 Iu(q) (Eq. S11) and undulation relaxation rates τu(q) (de-264 cay rate of undulations auto-correlation function, defined in 265 Eq. S12) capture out-of-plane motions, and, the lateral static 266 structure factor S(q) (Eq. S28) and the dynamic structure 267 factor S(q, ω) (Eq. S29) capture in-plane density fluctuations. 268 Fig. 2(A) shows that uL-NMS predictions for Iu(q), deter-269 mined from only the ENM parameters (Eqs. S11, S23), are in 270 excellent agreement with MD calculations for the whole range 271 of wave-vectors. A q −4 scaling at mesoscopic lengthscales 272 (q ≲ 1.5 nm −1 ) matches the prediction from a zero thickness 273 flat-sheet model of bilayer (12) (Fig. S8(A)). A more accurate 274 continuum model that includes the effects of bilayer bending, 275 lipid molecular tilt, and protrusions in membrane energetics 276 has been used to describe undulation spectra over a wide 277 q range, and estimate membrane bending (kc) and tilt (k θ ) 278 moduli as per Eq. S25 (79). Bilayer mechanical properties 279 calculated by a fit of Eq. S25 to uL-NMS undulations spectra 280 are in an excellent agreement with those determined from 281 the fit to MD data (this study, refs. (80, 81)) and neutron 282 spin-echo (NSE) spectroscopy (18, 82) ( Table 1). These elastic 283 properties of a bilayer are known to significantly influence its 284 function, stability, and behavior, for example, bending rigidity   The dynamic structure factor S(q, ω), a direct experi-305 mental observable that has been used extensively for lipid  Table 1. Bilayer properties from uL-NMS are compared with the MD simulation and the experimental data. The bending (kc) and tilt moduli (k θ ) are estimated from Iu(q) (Eq. S25), and thermal diffusivity (D T ) from Γ h (q) (Fig. 2(D)). space-domain Fourier transform of density correlation func-316 tion. Figs. S9 (E,F) shows that uL-NMS fails to reproduce 317 the short-time relaxation (∼ ps) that appears as a sharp decay 318 in both MD simulations and experiments (7), but accurately 319 captures the long-time relaxation behaviour (∼ ns-µs). This 320 deviation in short-time dynamics has been attributed to the 321 overly-smoothed harmonic approximation for a hierarchical, 322 rugged energy landscape of proteins (68,92), and such an en-323 ergy landscape has also been seen in lipid bilayers (30, 93, 94). 324 The long timescale applicability of uL-NMS is further vali-325 dated through the Rayleigh-Brillouin triplet model for S(q, ω) 326 and F (q, t) (Eq. S39), which is applicable in the hydrodynamic 327 limit and is used to analyze the data from NSE spectroscopy 328 (32, 91). A fit of Eq. S39 to S(q = 0.70 nm −1 , ω) from MD 329 of DMPC bilayer ( Fig. S8(C)) gives the width of Brillouin 330 and Rayleigh peaks, Γs = 0.3 ps −1 and Γ h = 0.44 ns −1 respec-331 tively, with Γs ≫ Γ h at all other q also (data not shown). 332 This implies that the Rayleigh peak corresponds to the long 333 time (∼ ns-µs), highly overdamped relaxation of the density 334 fluctuations while the Brillouin side-peaks correspond to the 335 short time (∼ ps), underdamped relaxations shown in Figs. S9 336 (E,F). Fig. 2(D) shows that Γ h obtained from fits to uL-NMS 337 predicted S(q, ω) are in an excellent agreement with the corre-338 sponding values from MD simulations and NSE spectroscopy 339 (88). Moreover, its q-dependence predicted from the hydro-340 dynamic theory,(32, 88, 91) viz. Γ h (q) = D T q 2 (D T : thermal 341 diffusivity), is found applicable for values from uL-NMS, MD, 342 and NSE, with estimates of D T in close agreement (Table 1). 343 Overall, the extensive comparison with AA and experimental 344 data establishes that uL-NMS can be used to study largescale 345 functional conformational changes in lipid bilayers, which are 346 expected to occur over long timescales.

347
ENM for a protein embedded in a lipid bilayer. γ-secretase em-348 bedded in a DPPC bilayer, a key therapeutic target for the 349 Alzheimer's disease (95), is taken as a model system for a 350 membrane protein (Fig. 3(A)). A heterogeneous environ-351 ment inherent to a membrane-embedded protein makes in-352 adequate the typical development of ENM without explicit 353 consideration of the environment degrees of freedom. This 354 approach, referred to here as ENM-P, has been used to suc-355 cessfully model dynamics of proteins and protein assemblies 356  S11(B and C)). This results from some select 395 differences in the set of slow modes in ENM-P and ENM-PL, 396 even when 6 out 10 slowest modes in two sets have > 75 % 397 overlap ( Fig. S11(A)). For example, there is a noticeable 398 underprediction of the conformational variance along the 3 rd 399 ENM-P mode, implying an inaccurate estimate of mode vi-400 brational frequency. This is also the first ENM-P mode that 401 has a low match with the corresponding ENM-PL mode. This 402 analysis of conformational variance in the normal mode space 403 clearly explains and establishes the importance of including 404 protein-lipid interactions in ENM development for membrane 405 protein.  Table 2). The next set of 430 contributing modes capture shear deformation (no associated 431 wave-vector) and undulations spanning the lateral size (q = 432 0.71 nm −1 ). The overalldpore projections shows a concerted, 433 synchronized motion of lipid beads spanning the entire bilayer 434 leading to a density reduction around the cylinder axis, i.e., an 435 increase in ξ ( Fig. S12(C)). Lateral diffusion in at-equilibrium 436 lipid bilayers, characterized from µs-long MD simulations, 437 was found to be governed by a predominantly continuous 438 motion with correlations spanning over tens of nanometers (6). 439 These flow patterns were conjectured to be important for the 440 molecular mechanisms of several processes in a bilayer and are 441 observed here for displacement in a plausible pore formation 442 direction (increasing ξ) in a free-standing bilayer.

443
Conformational sampling alongdpore is accelerated by 444 adding an extra velocity to the initial velocity distribution, 445 given by a weighted sum in the normal mode space, viz. 446 vextra = λ βmUm, where λ controls the overall degree of 447 excitation (excitation factor reported in K as ∆T ). In absence 448  with ∆T from 0 to 0.2 ( Fig. S13(A,C)). However, the initial (Eq. S44) and z(R) (Eq. S45), which respectively de-493 scribe the progress along and the distance from this pathway 494 (36, 104)(details in Supplemental Material Sec. F). The evolu-495 tion of the lipid bilayer in water system between a flat bilayer 496 and a prepore state, calculated from an umbrella sampling (US) 497 simulation along s(R) (details in Methods Sec. E), is reported 498 in terms of pore lipid density order parameter ξ (Eq. 2, Ref. 499 (44)) and pore water density order parameter ξ MF (Eq. S58 , 500 Ref. (45)). Fig. 4(A, B) (Fig. S15(A, B)) respectively show 501 that ξ and ξ MF profiles along the pore-opening and closing 502 pathways of DMPC (DPPC) bilayer are essentially coincident, 503 with maximum deviation observed near the transition state 504 (s(R) ≈ 0.6). The potential of mean force (PMF) profiles, 505 ∆Gpore, are also coincident along the two pathways for both 506 bilayer systems (Fig. 4(D), Fig. S15(D)). The activation bar-507 riers, Ea, along the two pathways are within 5 k B T for both 508 bilayer systems, and overall root-mean square deviation of 509 ∆Gpore is 3.2 k B T and 5.4 k B T for DMPC and DPPC bilayers, 510 respectively. In an earlier study, direct use of ξ to obtain a set 511  erties calculated using all atom MD simulations (lipid bead 572 RMSF, Fig. 1(A-B); undulation spectra, Fig. 1(C) and veloc-573 ity auto-correlation, Fig. 1(D)). Several independent (not used 574 in parameterization) properties estimated from uL-NMS, viz. 575 mechanical parameters, and equilibrium and long-timescale 576 dynamic properties, were shown to be in excellent agreement 577 with MD simulations from this and earlier studies, and experi-578 mental data (Fig. 2, Fig. S9, Table 1). Importantly, uL-NMS 579 parameterized from short simulations (2 ns to 100 ns) simu-580 lations of a small (256 lipid) bilayer was found in-agreement 581 with properties calculated from a long simulation (1 µs) of a 582 larger (1024 lipids) bilayer. While the long-time dynamics of 583 the DMPC bilayer were in overdamped regime, the short-time 584 (∼ ps) relaxation of density fluctuations was in underdamped 585 regime and were not captured by uL-NMS (Fig. S9 (E,F)). 586 This is a typical deficiency resulting from local smoothing of 587 the energy landscape in ENM representations (30, 68, 92-94). 588 The use of AA simulations for affecting structural change 589 along the direction predicted by the coarse-grained uL-NMS 590 representation allowed for relaxation of these fast degrees of 591 freedom during membrane remodeling simulations.

592
The as-parameterized lipid bilayer ENM was shown to carry 593 over to the development of ENM for a membrane-protein 594 system, with γ-secretase in a DPPC bilayer taken as model 595 system. The neglect of membrane environment in a protein-596 only ENM (ENM-P) gave a poor accuracy for the RMSF of 597 membrane-embedded residues (Fig. 3) and residue-residue 598 fluctuation correlations (Fig. S11), both of which improved 599 significantly on integration of the protein and lipid ENMs for 600 representation of overall Hamiltonian of the system (ENM-601 PL). This could serve as a valuable tool for characterization 602 of vibrational dynamics in membrane proteins, wherein the 603 surrounding membrane environment plays an integral role in 604 governing allosteric communication and signalling pathways 605 involving embedded proteins.

606
One prominent finding of this study is development of a 607 general framework for studying membrane remodeling, applied 608 here for determining the free energy landscape and molecular-609 level characterization of pore formation in DMPC and DPPC 610 bilayers. The connection between membrane remodeling and 611 intrinsic dynamics is exploited to estimate the direction of pore 612 formation (dpore) at the coarsegrained uL-NMS description 613 and actual structural change is carried in an AA simulation. 614 Several degrees of freedom are shown to contribute todpore in 615 addition to membrane undulations (Table 2, Fig. 5: thickness 616 fluctuations, lateral displacement, shear deformation, lipid 617 tilt), as expected given the involvement of nanometer scale 618 features(106). The inherent connection between uL-NMS and 619 the underlying free energy landscape (unknown at outset) 620 allowed determination of a generic reaction coordinate (PCV 621 s(R)) by affecting small structural change along an arbitrary 622 lipid density based progress variable, ξ, whose direct use as a 623 RC led to hysteresis in pore formation and closing pathways 624 (46). In contrast, the PMF along s(R) is found to be coincident 625 on the forward and reverse pathways, with free energy profiles 626 and activation barriers in agreement with recent findings (5, 627 43, 45).

628
The metastability of the prepore state governs the transport 629 of hydrophilic molecules and ions across the lipid bilayer(42). 630 For both DMPC and DPPC bilayers, pore closure is observed 631 in only one out of the five independent, unbiased 500 ns simula-632 tions started from the pore state (Fig. S18(A-B)). In contrast, 633 DPPC prepore states identified in some other studies were 634 found to close within a ∼ 10 ns MD simulation (5, 42). These  Fig. S19(U)). It represents a net outward migration of 668 lipid atoms away from the pore centre, consequently reducing 669 lipid density in this area (Fig. S12(C)). The gradient of t XY 670 alongdpore (Figs. 5(B), S19) clearly shows that membrane 671 thinning is predicted by the uL-NMS model, and is not merely 672 a sequential after-effect of a decrease in lipid density near 673 the pore centre. Evolution of lipid orientation at any bilayer 674 surface grid (details in Supplemental Material Sec. G.2) was 675 defined by the projection of a vector from tail to headgroup, 676 n XY , on the vector connecting the grid to the pore centre, 677 ⟨n XY ·r XY ⟩ (Figs. 5(D), S19(P-T)).

678
The average orientation of lipids within the pore-spanning 679 cylindrical region changes from being along the bilayer normal 680 (w.r.t. flat bilayer) to along the normal to the developing 681 curved surface near the indentation in almost direct proportion 682 to the extent of bilayer thinning near the pore center ( Fig. 683  5(E)). This is indicative of a cooperative shielding of alkyl 684 tails in an indented membrane by lipid headgroups (39).

685
The thinning of the lipid bilayer also correlate with an 686 increase in the coordination between the lipid headgroup po-687 lar atoms (LP: N + and all PO 3− 4 group atoms, Fig. S16(A)) 688 and water oxygen atoms (Fig. 5(G), Fig. S19(V)). This was 689 accompanied by a decrease in the inter-lipid polar-polar coor-690 dination, maintaining a nearly constant value for total polar 691 atom coordination of lipid headgroup polar atoms throughout 692 the pore formation pathway (see Fig. S19(W)). The average 693 fractional contribution of water molecules to the total polar-694 polar coordination of lipid headgroup atoms within the pore 695 spanning cylinder (⟨f LP-water ⟩ cyl ) increased from 0.81 in the 696 flat-bilayer state to 0.89 in the stable pore state. The number 697 of water molecules coordinated to only lipid tails, i.e., exclud-698 ing those coordinated to both lipid tails and polar headgroups, 699 N tail w , is used to monitor presence of a hydrophobic defect 700 (Supplemental Material Sec. G.3). Fig. 5(H) shows that 701 N tail w ≈ 0 in the initial stages of pore formation (s(R) ≈ 0. 3). 702 The PMF profiles shows a large increase in the free energy 703 till this stage (∼ 10 k B T from Fig. 4(D)), essentially equal 704 to the barrier for pore formation, and therefore, these "de-   Fig. S17(B), resulting 768 in lipid-lipid LJ penalty. In case of DPPC, this penalty is 769 higher compared to DMPC bilayer due to high packing order 770 in its flat-bilayer state (Fig. S16(C)) and larger radial extent 771 of order disruption from pore centre (Fig. S17(B)), resulting 772 in net unfavourable enthalpic contribution of pore state.

D R A F T
where, X and V are the mass-weighted coordinate and velocity vec- where, r i is the radial distance (from pore center) of the i th