Stress Propagation through Biological Lipid Bilayers in Silico

Communication pubs.acs.org/JACS

Stress Propagation through Biological Lipid Bilayers in Silico Camilo Aponte-Santamaría,*,†,‡,¶ Jan Brunken,§ and Frauke Graẗ er*,†,‡ †

Molecular Biomechanics Group, Heidelberg Institute for Theoretical Studies, 69118 Heidelberg, Germany Interdisciplinary Center for Scientific Computing, Heidelberg University, 69120 Heidelberg, Germany ¶ Max Planck Tandem Group in Computational Biophysics, University of Los Andes, 111711 Bogotá, Colombia § Faculty of Biosciences, Heidelberg University, 69120 Heidelberg, Germany ‡

S Supporting Information *

Here, by using molecular dynamics (MD) simulations, we address the pressing question over which length scales stress propagates through a lipid bilayer, given that damping is expected to dominate in condensed phase. We considered a fully solvated and equilibrated dipalmitoylphosphatidylcholine (DPPC) lipid bilayer of length Lx ∼ 288 nm and cross-sectional area LyLz ∼ 8 × 11 nm2 (Figure 1A). We used the coarsegrained MARTINI model20 to describe the bilayer and simulated its dynamics using the GROMACS MD package.21 We considered two lipid slabs (1.5 nm each). Each slab was

ABSTRACT: Membrane tension plays various critical roles in the cell. We here asked how fast and how far localized pulses of mechanical stress dynamically propagate through biological lipid bilayers. In both coarse-grained and all-atom molecular dynamics simulations of a dipalmitoylphosphatidylcholine lipid bilayer, we observed nanometer-wide stress pulses, propagating very efficiently longitudinally at a velocity of approximately 1.4 ± 0.5 nm/ps (km/s), in close agreement with the expected speed of sound from experiments. Remarkably, the predicted characteristic attenuation time of the pulses was in the order of tens of picoseconds, implying longitudinal stress propagation over length scales up to several tens of nanometers before damping. Furthermore, the computed dispersion relation leading to such damping was consistent with proposed continuum viscoelastic models of propagation. We suggest this mode of stress propagation as a potential ultrafast mechanism of signaling that may quickly couple mechanosensitive elements in crowded biological membranes.

M

echanical stress defines key aspects in biological membranes such as shape, deformability, dynamics, and fusogenicity.1−3 In particular, mechanical stress serves as a stimulus for the activation of mechanosensitive membraneembedded ion4,5 and possibly water channels.6 A long-standing question is whether the lipid bilayer by itself provides a pathway to dynamically transmit the forces necessary to activate these channels, in the absence of auxiliary scaffold proteins, and if that was the case, how fast and how far membrane-mediated mechanical stress propagates. Ultrasound experiments suggested a very fast mode of propagation through solvated lipid bilayers (comparable to the speed of sound in water).7,8 In addition, computer simulations linked membrane undulations and fluctuations to key elastic and thermodynamic properties for the propagation of stress.9−13 Furthermore, sound has been proposed to exhibit soliton-like modes of propagation in lipidic systems, due to their nonlinear behavior near their melting temperatures,14 an effect observed recently in lipid monolayers at water−air interfaces.15,16 Mechanical stress upon a perturbation has also been shown to propagate through proteins for few nanometers within picoseconds,17,18 resulting in directed vibrational modes with potential functional implications.19 © 2017 American Chemical Society

Figure 1. Longitudinal stress-propagation through a DPPC lipid bilayer revealed by force-probe MD simulations. (A) Coarse-grained system highlighting the lipid acyl-chains (green) and headgroups (brown and blue). Solvent not shown for clarity. Two portion of lipids at the extremes of the simulation box (blue) were subjected to a force F = 5 nN during the first 7 ps of the simulation. (B) Punctual stress PS as a function of the time t for each position of the bilayer x. PS is presented according to the color bar, highlighting reductions (brownred) and increments (blue-black) of up to 3 standard deviations (σ) with respect to the total average background value ⟨⟩ at t = 0 (green). At the right panels, PS versus x is presented for a time window of 5 ps centered at times 10, 20, ..., 60 ps (gray, raw data; black, smoothed curve; green dot, highest peak at each side of the box). Height of each panel is 12σ. The inset shows the PS for a smaller system with Lx ≃ 72 nm (Ly and Lz are similar to panel A). Received: May 8, 2017 Published: August 30, 2017 13588

DOI: 10.1021/jacs.7b04724 J. Am. Chem. Soc. 2017, 139, 13588−13591

Communication

Journal of the American Chemical Society

increase in lipid density, is due to the weak dependency of c on [L] imposed by the experimental u values (∼10−4 m3/kg). Therefore, the speed of propagation of our simulated pulses (1.4 ± 0.5 km/s) is in remarkable good agreement with the experimental estimations of the speed of sound through DPPC bilayers. We also monitored the attenuation of the pulses by extracting their maximum amplitudes as a function of time (Figure 2B). Pulses decayed within tens of picoseconds. More quantitatively, by fitting an exponentially decaying function exp(−t/τ), also accompanied by a bootstrap analysis to estimate the error, yielded a characteristic decay time τ (for which the amplitude decayed by 1/e of its initial value) of approximately 70 ± 30 ps. Therefore, we predict that stress pulses travel distances of the order of cτ = 100 ± 50 nm before attenuation. To discard possible artifacts due to the coarse-grained level of resolution, we also induced the pulse propagation through all-atom DPPC bilayers.23 All-atom MD simulations were conducted for smaller systems (x1, Lx ∼ 36 nm; x2, Lx ∼ 72 nm) and compared to coarse-grained simulations of the same as well as of larger systems (x4, Lx ∼ 144 nm; x8, Lx ∼ 288 nm). Similar propagation patterns were observed for both levels of resolution (Figure S1A). As expected, the time until pulses encountered one another increased with the system size. For smaller systems, the pulses encountered themselves several times before vanishing and the dissipation after encounter was slightly less pronounced for the coarse-grained systems (see e.g. x2 size in Figure S1A). For x1, the pulse width was almost comparable to the box size and pulses met almost instantaneously, thereby preventing from a reliable estimation of c, τ, and cτ. For the other sizes, these parameters approached plateau values with increasing system size (Figure 3A). The estimates from coarse-grained simulations approached the allatom predictions when reducing the system size (Figure 3A). Increments of the cross-sectional area by a factor of two or four, while maintaining the applied pressure constant (see methods in Supporting Information), did not significantly affect neither the pulse patterns (Figure S1B) nor the c, τ and cτ parameters derived from them (Figure 3B). Similar results were observed in the all-atom simulations (which had fewer available water molecules per lipid molecule than the CG simulations) and in extra coarse-grained simulations with increased Lz (Figure S1C), indicating that the size of the water solvation shell did not artificially restrict dissipation. Therefore, the pulses are robust with regard to the system size (length, crosssectional area, and height) and to the level of resolution (allatom or coarse-grained). Pulses were also observed after reducing the effective applied force per lipid, Flipid (Figure S2). In this case, the initial amplitude of the pulse was reduced as Flipid was decreased, and this caused the pulses to propagate for smaller distances, but around the same speed. Note that even for the smallest Flipid ∼ 33pN, pulses propagated by at least 10 nm before attenuation. Thus, pulses are predicted to propagate of the order of tens of nanometer, provided that the initial perturbation is sufficiently larger than the background average stress (achieved here by initially pushing with Flipid > 30 pN). We also obtained similar results for a more characteristic lipid in biological membranes (DOPC) at a physiological temperature of 310 K (Figure S3A), when changing the boundary conditions (Figure S3B), or in absence of a barostat or thermostat (Figure S3C), corroborat-

Figure 2. Speed and attenuation of stress pulses propagating longitudinally through a DPPC lipid-bilayer. (A) Displacement and (B) amplitude of the PS peak as a function of time. Normalized displacement with respect to a reference position xref = 145 nm is presented in panel A. Two peaks (one at x < 0 and one at x > 0) were considered for each time (dots). Fits to the MD data are presented with the black line. Gray area displays the error of the fits obtained by bootstrap analysis (80% confidence intervals). See details of the smoothing and bootstrap analysis in Supporting Information.

located near each box periodic boundary and was oriented perpendicular to the longest bilayer axis. We pushed the slabs longitudinally in opposite directions, by exerting a force F = 5 nN on them during the first 7 ps of the simulation (Figure 1A). This force ensured a sufficiently high signal-tonoise ratio. Because F was distributed on all the lipids contained in the pulled portions (35 in average), the effective force per lipid was of the order of 140 pN, closer to the pN biological range. We released the force and monitored the spatiotemporal changes in mechanical stress by computing the punctual stress PS (the spatially and temporally resolved sum of all interlipid pairwise forces) from force distribution analysis.22 See details of simulations and punctual stress calculation in the Supporting Information. By averaging the punctual stress over 12 independent trajectories, at each position of the bilayer x and at each time t, PS(x,t) = (1/12)∑r12= 1PSr(x,t), we detected the formation of two mechanical pulses (Figure 1B). They were initially of approximately 5 nm in width, and later propagated toward each other traveling tens of nanometers in tens of picoseconds. With their amplitude decreasing while the width increased, the pulses dispersed and almost vanished at around 100 ps, shortly before they met in the middle of the box (Figure 1B). In simulations with a 4-fold reduced system size, pulses continued propagating while maintaining their shape, after repeatedly crossing with each other at the middle of the box and at the periodic boundaries (Figure 1B inset). Extraction of the time-dependent position of the pulse peaks before they encountered for the fist time, followed by a smoothing at each time window to improve the signal-to-noise ratio, allowed us to quantify the speed of propagation. A linear regression to the curve of the peak positions versus time, together with a bootstrap analysis, revealed that pulses propagate at a constant speed c of ∼1.4 ± 0.5 nm/ps (km/s) (Figure 2A). Ultrasonic velocity measurements7,8 yielded experimental estimates of the speed of sound of around 1.54−1.58 km/s (see details of the estimation in the Supporting Information). Note that in these cases the speed is (almost) equal to the speed of the surrounding medium, which is not surprising given the highly diluted conditions of these experiments ([L] of the order of 1 kg/m3). If we extrapolate to a higher lipid concentration such as that of our simulations ([L] = 576 kg/m3), we obtain a 10−15% reduction in c, yielding an estimate from experiments and simulations of 1.39 km/s. This small reduction, despite of the 100-fold 13589

DOI: 10.1021/jacs.7b04724 J. Am. Chem. Soc. 2017, 139, 13588−13591

Communication

Journal of the American Chemical Society

Figure 4. Spectral analysis of the stress pulses. (A) Spectrum of the punctual stress PS(x,t) recovered by two-dimensional Fourier analysis. Harmonics are sorted according to their amplitude A (a.u.: arbitrary units). (B) Dispersion relation (κ vs ω) recovered from the Fourier analysis of PS(x,t) recovered from the MD simulations (colored squares). The amplitude A of each harmonic mode is color-coded based on the color scale at the bottom. A weighted power fit using main harmonics (A > Amin = 1.3 a.u.) is presented with the red line (α = 0.53 ± 0.10 nm−1 ps0.76 and n = 0.76 ± 0.30). (C) Dependence of the resulting fit parameters F and n on the cutoff amplitude Amin is displayed with the colored solid lines. Average over values obtained between 1.3 ≤ Amin ≤ 6 a.u. is indicated with colored dashed lines. Gray areas indicate estimated errors from bootstrap analysis (80% confidence intervals). Black dashed lines in panels A and C indicate the cutoff amplitude Amin used to compute the dispersion relation in panel B and throughout the text. Details of the Fourier analysis, the weighted-fit and the bootstrap analysis are in the Supporting Information.

Figure 3. Stress propagation in all-atom and coarse-grained systems. (A,B) Propagation speed c, attenuation time τ, and propagation distance cτ are shown as a function of the bilayer length Lx (A) and cross-sectional area Ly (B). (C) c and τ for the solvated DPPC lipid bilayer is compared to the values for pure coarse-grained water (H2O) and hexadecane (Hexa.). Experimental c values mentioned in the text are also displayed. Lx ≃ 144 nm and Ly ≃ 8 nm were used here.

ing our finding of longitudinal stress propagation in lipid bilayers. To investigate the molecular factors defining the mode of propagation, stress was induced to propagate through two pure media: water and hexadecane (Figures 3C and S4). The MARTINI force-field predicted the speed c in very good agreement with the experimental values at 320 K.24,25 Furthermore, c of solvated DPPC was found to be closer to that of water than to that of hexadecane (Figure 3C,top). Concerning attenuation, τ for the bilayer was comparable to the values for water and hexadecane (within the uncertainty bars in Figure 3C,bottom). Compensation of the excess density of the headgroups by the acyl-chains may be the key factor defining similar stress propagation in the bilayer and in the bulk media (see below). We compared our results with a recent continuum model developed by Kappler and Netz for the propagation of pressure waves through materials coupled viscoelastically to their surrounding media. 26,27 They obtained the following generalization of the Lucassen dispersion relation 28 relating the wavelength κ and the wave-frequency ω: κ 2 = (ω 2 /K 2D) × ( iρη /ω + ρ2D ). Here, ρ is the density and η the viscosity of the surrounding bulk media. The bilayer is coupled to these media through the area elastic modulus K2D and the bilayer excess density ρ2D. Experimental and computational estimates of these quantities9,11,29,30 (see Supporting Information) indicate that the DPPC bilayer at 320 K appears to be much more similar to the surrounding bulk water, regarding the density, and thereby the dispersion relation approximates to |κ| ∼ αω3/4, with α = (ρη/K2D2)1/4 (1.6−1.8 nm−1 ps3/4). We computed the dispersion relation by applying a twodimensional Fourier transform to the time- and spatially-

resolved punctual stress PS(x,t) (Figure 4). At the most, few hundred harmonics contributed to the PS spectrum (Figure 4A). By considering only high-amplitude harmonics (A ≥ Amin, with Amin = 1.3 a.u., which is 16% of the maximum amplitude Amax = 8 a.u.), a weighted fit to the spectrum yielded |κ| = αωn, with α = 0.53 ± 0.10 nm−1 ps0.76 and n = 0.76 ± 0.30 (Figure 4B). Note that these fitting parameters did not substantially vary depending on the number of harmonics considered, as α and n stabilized around constant values for increasing Amin (Figure 4C). The estimate of the power n (0.76 ± 0.30) is encouragingly close to the factor 3/4 predicted by the theory. Our estimate of α (0.53 ± 0.10 nm−1 ps0.76) appears to slightly underestimate the range roughly predicted above from thermodynamic quantities (1.6−1.8 nm−1 ps3/4). Deviations are acceptable given the coarse level of resolution and the continuum nature of the analytic model. The good agreement with experimental propagation speeds suggests that thermodynamic quantities such as density and elastic compression are well captured by the simulations. On the contrary, the inherent difficulties of MD simulations to properly account for the diffusion (and thereby for the viscosity) of interfacial lipidic systems12,13 may affect the damping estimates. Furthermore, transverse membrane undulations (not considered here) may also play an important role smearing out the lipid excess density. Finally, the continuum model may not fully capture molecular details at the headgroup−water interface which could become relevant at the studied nanometer−picosecond scale. 13590

DOI: 10.1021/jacs.7b04724 J. Am. Chem. Soc. 2017, 139, 13588−13591

Communication

Journal of the American Chemical Society The group velocity vg(κ0) = ∂ω(κ0)/∂κ for the highest harmonic (κ0, ω0) = (0.13 nm−1, 0.17 ps−1), using the fitting parameters α and n, was found to be 1.6 km/s. This value is consistent with the propagation speed directly computed from the positions of the pulses (∼1.41 ± 0.5 km/s). Our analysis in the reciprocal space therefore suggests that pulses can be semiquantitatively explained by the continuum viscoelastic models proposed by Kappler and Netz26,27 and Lucassen28 and are in relatively good agreement with existing experimental and simulation data of DPPC bilayers. In conclusion, we here present simulation evidence for the propagation of nanometer-short pulses through biological lipid bilayers at speeds of the order of nm/ps (km/s) during tens of picoseconds implying propagation distances up to tens of nanometers. Our speed estimates are in very good agreement with estimates from ultrasound experiments. Furthermore, our data is robust with regard to the level of resolution (coarsegrained or all-atom) and the simulation details. The predicted dispersion relation, which leads to the observed speeds and attenuation properties, is consistent with existing theories of pulse propagation through viscoelastic continuum materials.26−28 Our work supports that the lipid bilayer can by itself couple mechanosensitive elements in crowded biological membranes, by quickly transmitting localized pulses of stress. We speculate the bilayer to confer an advantage by providing a narrow two-dimensional medium for directed-stress propagation. It will be highly interesting to test the dispersion of these pulses experimentally (by, e.g., terahertz spectroscopy or neutron scattering methods), to check whether they may exhibit nonlinear properties such as the ones observed for monolayers, and to elucidate their biological relevance for multicomponent complex membrane systems.

■

(3) Kozlov, M. M.; Chernomordik, L. V. Curr. Opin. Struct. Biol. 2015, 33, 61−67. (4) Kocer, A. Curr. Opin. Chem. Biol. 2015, 29, 120−127. (5) Ranade, S.; Syeda, R.; Patapoutian, A. Neuron 2015, 87, 1162− 1179. (6) Fischer, G.; Kosinska-Eriksson, U.; Aponte-Santamaría, C.; Palmgren, M.; Geijer, C.; Hedfalk, K.; Hohmann, S.; de Groot, B. L.; Neutze, R.; Lindkvist-Petersson, K. PLoS Biol. 2009, 7, e1000130. (7) Mitaku, S.; Ikegami, A.; Sakanishi, A. Biophys. Chem. 1978, 8, 295−304. (8) Schrader, W.; Ebel, H.; Grabitz, P.; Hanke, E.; Heimburg, T.; Hoeckel, M.; Kahle, M.; Wente, F.; Kaatze, U. J. Phys. Chem. B 2002, 106, 6581−6586. (9) Lindahl, E.; Edholm, O. Biophys. J. 2000, 79, 426−433. (10) Ayton, G.; Smondyrev, A. M.; Bardenhagen, S. G.; McMurtry, P.; Voth, G. A. Biophys. J. 2002, 82, 1226−1238. (11) Marrink, S. J.; de Vries, A. H.; Mark, A. E. J. Phys. Chem. B 2004, 108, 750−760. (12) den Otter, W.; Shkulipa, S. Biophys. J. 2007, 93, 423−433. (13) Venable, R. M.; Ingólfsson, H. I.; Lerner, M. G.; Perrin, B. S.; Camley, B. A.; Marrink, S. J.; Brown, F. L. H.; Pastor, R. W. J. Phys. Chem. B 2017, 121, 3443. (14) Heimburg, T.; Jackson, A. D. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 9790−9795. (15) Griesbauer, J.; Bössinger, S.; Wixforth, A.; Schneider, M. F. Phys. Rev. Lett. 2012, 108, 198103. (16) Shrivastava, S.; Schneider, M. F. J. R. Soc., Interface 2014, 11, 20140098. (17) Young, H. T.; Edwards, S. A.; Gräter, F. PLoS One 2013, 8, e64746. (18) Brinkmann, L. U. L.; Hub, J. S. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 10565−10570. (19) Niessen, K. A.; Xu, M.; Paciaroni, A.; Orecchini, A.; Snell, E. H.; Markelz, A. G. Biophys. J. 2017, 112, 933−942. (20) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. J. Phys. Chem. B 2007, 111, 7812−7824. (21) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. SoftwareX 2015, 1−2, 19−25. (22) Costescu, B. I.; Gräter, F. BMC Biophys. 2013, 6, 5. (23) Ulmschneider, J. P.; Ulmschneider, M. B. J. Chem. Theory Comput. 2009, 5, 1803−1813. (24) Bilaniuk, N.; Wong, G. S. K. J. Acoust. Soc. Am. 1993, 93, 1609− 1612. (25) Outcalt, S.; Laesecke, A.; Fortin, T. J. J. Chem. Thermodyn. 2010, 42, 700−706. (26) Kappler, J.; Netz, R. R. EPL 2015, 112, 19002. (27) Kappler, J.; Shrivastava, S.; Schneider, M. F.; Netz, R. R. ArXiv eprints; https://arxiv.org/abs/1702.08864, 2017. (28) Lucassen, J. Trans. Faraday Soc. 1968, 64, 2221−2229. (29) Kestin, J.; Sokolov, M.; Wakeham, W. A. J. Phys. Chem. Ref. Data 1978, 7, 941−948. (30) Nagle, J. F.; Tristram-Nagle, S. Biochim. Biophys. Acta, Rev. Biomembr. 2000, 1469, 159−195.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/jacs.7b04724. Estimation of experimental sound speeds, dispersionrelation parameters, methods, and figures (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*[email protected] *[email protected] ORCID

Frauke Gräter: 0000-0003-2891-3381 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Julian Kappler, Matthias Schneider, Roland Netz, and Alfredo Alexander-Katz for helpful discussions. This work was supported with funds from the Deutsche Forschungsgemeinschaft, grant FOR1543 (to C.A.-S. and F.G.), the BIOMS postdoctoral program of the Heidelberg University (to C.A.-S.), and the Klaus Tschira Foundation (to F.G.).

■

REFERENCES

(1) Gauthier, N. C.; Masters, T. A.; Sheetz, M. P. Trends Cell Biol. 2012, 22, 527−535. (2) Diz-Muñoz, A.; Fletcher, D. A.; Weiner, O. D. Trends Cell Biol. 2013, 23, 47−53. 13591

DOI: 10.1021/jacs.7b04724 J. Am. Chem. Soc. 2017, 139, 13588−13591