Modelling the dynamics of vesicle reshaping and scission under osmotic shocks

Abstract

We study the effects of osmotic shocks on lipid vesicles via coarse-grained molecular dynamics simulations by explicitly considering the solute in the system. We find that depending on their nature (hypo- or hypertonic) such shocks can lead to bursting events or engulfing of external material into inner compartments, among other morphology transformations. We characterize the dynamics of these processes and observe a separation of time scales between the osmotic shock absorption and the shape relaxation. Our work consequently provides an insight into the dynamics of compartmentalization in vesicular systems as a result of osmotic shocks, which can be of interest in the context of early proto-cell development and proto-cell compartmentalisation.

---

1 Introduction

Lipid membranes are one of the most important building blocks in cellular biology, forming the envelope of the cell and defining the structure and function of intracellular units in eukaryotic cells. Given their physio-chemical properties, biological membranes present very particular reshaping capabilities that are intimately related to the proper functioning of the cellular units they constitute. Cell membrane morphology transformations in nature can occur in a number of different ways, which can also include creating of new units such as budding or division. Budding is the process by which a portion of a membrane detaches, usually in the form of a vesicle that is smaller than the mother membrane. Such processes are ubiquitous in cellular biology, playing important roles in vesicle secretion, protein transport, organelle synthesis and endocytosis.¹ In the division, a lipid vesicle or a cell splits into two closed units of comparable sizes. In more evolved cells this process is usually very specifically driven by a set of proteins (such as actomyosin,² ESCRT-III^3,4 or FtsZ⁵) composing the divisome machinery, responsible for exerting the required forces on the membrane. However, division of vesicles can also occur in spontaneous ways usually caused by environmental changes (shear stress,⁶ osmotic shocks^7–9), compositional changes to the membrane (protein addition^10,11) or even motility of nascent cells.¹²

In this context, vesicles are especially relevant as they represent a typical model system to study membrane reshaping events. Moreover, in the context of early life evolution and proto-cells development in particular, vesicles are thought to have played a crucial role, having evolved towards the existing forms of cellular life by compartmentalizing and developing more complex behaviors as a result.^13–15 There exist a number of experimental, theoretical and numerical works that examine the connection of distinct vesicle reshaping events to their physio-chemical conditions. Previous experimental efforts have focused on understanding the role of the physical properties (e.g. fluidity,¹⁶ multiphase domain coexistence^17,18) and composition (e.g. protein adsorption at high^10,19 and low¹¹ surface density) of membranes in reshaping Giant Unilamellar Vesicles (GUVs). Other experimental approaches have focused on the effect of osmotic shocks and other external forces (e.g. shear stress⁶) on vesicle morphology. It is thus known that hypertonic shocks (when the solute concentration is larger outside than inside the vesicle) induce strong deformations on GUVs with varying results depending on the specific composition and chemical properties of the membrane.^9,20,21 Division can hence occur when osmotic stress is coupled with phase separation in the lipid bilayer.^7,8 More important for cell compartmentalization and endocytosis, simple osmotic shocks can also drive inward budding of the membrane coupled with absorption of external medium into the newly created inner compartment.^{9,14,15,22,23} Finally, studies have also shown that hypotonic shocks result in pulsatile dynamics of stretch and release for the vesicle until equilibrium is reached between the two solutions.^24,25

Regarding theoretical approaches, most advances came at the end of the 20th century, with the introduction of the Helfrich model of bending free energy.²⁶ Subsequently, most experimentally observed vesicles conformations were described by applying free energy minimisation principles focusing primarily on the membrane curvature, the so-called curvature models.^27,28 Although very effective in describing equilibrium shapes, such theoretical approaches present some limitations. Probably the most important limitation comes from working with continuous surfaces, which makes it impossible to describe budding or membrane scission events, resulting in unphysical limit shapes. Furthermore, given that usually such models assume constant surface area, no stretching of the membrane is allowed, rendering any hypotonic shock scenario virtually inaccessible. Indeed, hypotonic shocks result in a swelling and bursting of the vesicle sometimes displaying cyclic relaxation dynamics^24,25 which cannot be replicated by a non-stretching surface description of the lipid bilayer that cannot break. Finally, by the very nature of the free energy minimisation approach, curvature models provide little insight into the relaxation dynamics of the vesicles.

An alternative approach are physiochemical kinetic models,^29–31 which do offer a solution for vesicle volume dynamics, but do not capture the shape transformations associated with it. Lastly, existing numerical efforts in trying to describe vesicle shape transformations are usually limited to Monte Carlo simulations of meshed surfaces governed by energy considerations associated with curvature.³² These are again successful in predicting most equilibrium shapes, but provide little help in overcoming the limitations presented by the theoretical models à la Helfrich. One exception to these examples, which manages to reproduce budding events, is the approach taken by Yuan et al.³³ The authors used molecular dynamics simulations to study the effects of volume changes resulting from osmotic shocks on various types of vesicles. One limitation of this study is that, by explicitly controlling the volume dynamics via a prescribed protocol, an artificial relaxation rate was introduced.^29,30

In the present work we focus on vesicle reshaping as a result of explicit osmotic shocks. In order to overcome limitations posed by the previous methods, where reshaping dynamics and vesicle breakage were inaccessible, we resort to using molecular dynamics simulations of coarse-grained model lipid bilayers under an explicit gradient of solute particles. This approach allows us to characterize the dynamic process leading to the different morphological transformations of vesicles. We find that hypertonic osmotic shocks can result in the engulfment of external material into inner vesicular compartments while hypotonic shocks on the contrary cause vesicle dilation and oscillatory bursting for sufficiently strong concentration differentials. We observe a separation of time scales in the dynamics of the shock absorption and shape transformations, where the vesicle first absorbs shocks by forming many local deformations and transient shapes, before relaxing to an equilibrium shape. Our study provides new insights into the dynamics leading to vesicle deformations and, in particular, vesicle endocytic events as a result of environmental osmotic stresses.

2 Simulation details

In order to study the effect of osmotic shocks on lipid vesicles we carry out coarse-grained (CG) molecular dynamics (MD) simulations where we explicitly consider the solute modeled as CG particles. The system consists of a spherical membrane surrounded by solute particles on the outside at a concentration ρ_out (see purple particles in Fig. 1) and encapsulating inner solute particles at a different concentration ρ_in (see cyan particles in Fig. 1). All particles are volume-excluded spheres of a diameter σ, σ being the MD unit of length, and mass m = 1, embedded in cubic box of constant volume V = L³, with L = 44σ and within periodic boundary conditions. See Fig. 1 for the snapshots of the system.

Fig. 1 Molecular dynamics simulation setup. (a) Simplified representation of the system: a spherical vesicle of CG lipids (yellow particles) in a cubic box with solute particles inside (cyan) at a concentration ρ_in and solute particles outside (purple) at a concentration ρ_out. (b) Cross-section snapshots of the simulated system at initial conditions (τ = 0 time steps) and equilibrium (τ = 3 × 10⁵ time steps) for representative parameters. The outer solute concentration is ρ_out = 0.18 and the concentration ratio is ν = 0.31. The complete transformation is shown in Supplementary Movie 2b (ESI†).

The membrane is modelled using the single particle thick model developed by Yuan et al.³⁴ which is capable of fission and fusion events and reproduces biologically relevant mechanical properties of membranes. In this model, bilayer sections of around 10 nm in width are represented by single beads interacting via an attractive potential. The bending energy is implemented via interactions between vectors on each bead, which represent the local normal of the membrane and are directly related to its curvature. The specific parameters of the model are chosen to encode for a fluid membrane of vanishing spontaneous curvature and bending rigidity of 15k_BT, k_B being the Boltzmann's constant and T being temperature. For more details on the interaction potential used in this model and the particular values used here refer to Section VI in the ESI.† All other interactions are modeled as simple volume exclusion repulsive interactions using a Weeks–Chandler–Anderson potential.

The system is always initialised in the same way: (i) N_mem = 4322 membrane particles define a relaxed spherical vesicle of mean radius R ∼ 17σ at the center of the box and (ii) solute particles take random positions on a hexagonal lattice at set concentrations ρ_in and ρ_out inside and outside of the vesicle, respectively, ensuring no overlap occurs with the membrane. In this way, the osmotic conditions of the system are fully defined by two parameters only: the outer solute concentration, ρ_out, and the ratio of concentrations ν = ρ_in/ρ_out. The latter quantity is often referred to as the theoretical reduced volume. Hence, setting ν < 1 defines hypertonic conditions while ν > 1 corresponds to hypotonic conditions. In this setup the vesicles studied have a diameter of a few hundred nanometers, which is within the range of physiologically relevant scales. However, to test the role of the vesicle size, we also explored larger model vesicles (see Section VII in the ESI†). We run simulations using the open source molecular dynamics package LAMMPS³⁵ with a Langevin thermostat (at reduced temperature T = 1 and damping coefficient of 1) within the NVE ensemble, E being the total energy of the system, for 3 × 10⁵ time steps, each of size 0.01 MD time units ( with ε = k_BT the MD energy unit).

3 Results

3.1 Characterising osmotic shocks

First we make sure that explicitly simulating solute particles in a closed system with a coarse-grained model lipid vesicle properly reproduces osmotic shocks. For this purpose we measure the initial osmotic pressure measured across the vesicle (see Section I in the ESI† for details). Fig. 2a presents the dependence of the osmotic shock on the two system parameters: the outer particle absolute concentration, ρ_out, and the ratio of concentrations between the vesicle interior and exterior, ν = ρ_in/ρ_out. Fig. 2b shows the agreement between the measured osmotic pressure with the expected value of the osmotic pressure, as given by the equation Π = RTρ_out(ν − 1).²⁷ Indeed, beyond a transformation factor due to MD units, we find that we are truly applying an osmotic shock to the vesicle and its magnitude scales as expected theoretically.

Fig. 2 Osmotic shock at initial conditions. (a) Measured values of the initial osmotic pressure difference across the membrane (osmotic shock) Π (color coded) as a function of the two system parameters {ρ_out, ν} for hypertonic conditions ν < 1. Values are in MD units of pressure: [Π] = εσ⁻³. (b) Measured values of the osmotic shock (orange dots, in MD units) against the corresponding expected theoretical value Π/(RT) = ρ_out(ν − 1). A linear fit is shown by the dashed line. The left panel shows the average over 10 independent runs for each parameter set while the right panel shows all the individual measurements.

3.2 Equilibrium vesicle morphologies

We are next interested in investigating whether our system is able to reproduce the relaxation transformations observed experimentally on GUVs,^9,15 predicted theoretically,^26,28 or simulated by other means.^32,33 We thus look at the equilibrium shapes our simulated vesicles present in hypertonic conditions (larger solute concentration outside than inside the vesicle). While the spectrum of shapes and different morphologies lipid vesicles present when subject to hypertonic shocks is vast and not necessarily restricted to a few discrete options, we can nonetheless categorise our results into several distinct kinds of shapes that help analyze the results. Inspired by previous works using similar classifications^15,28,32 we distinguish six different equilibrium vesicle morphologies according to which we classify our results (see Fig. 3 and accompanying snapshots):

Fig. 3 Morphology diagram for hypertonic conditions. We distinguish six different characteristic equilibrium shapes of the vesicle, represented in the right panel. The resulting observed morphology for each point in the parameter space {ρ_out, ν} is color coded accordingly. Dotted lines indicate transitions. The cross and star symbols refer to the dynamics curves below (see Fig. 4 and 5). Snapshots are framed using the same color code as the diagram on the left. Particle coloring is used to distinguish multiple vesicles in the same system. All morphologies correspond to equilibrium configurations, reached after τ = 3 × 10⁵ time steps.

• Burst: at high outer solute concentration ρ_out and very low ratio of inner and outer concentrations ν, the osmotic shock causes the mother vesicle to burst, often yielding several smaller vesicles containing a mixture of both inner and outer solutions (Supplementary Movie 2a, ESI†).

• Inner bud: for a large subset of the parameter space where the ratio of inner to outer concentrations is still small (ν < 0.5), but the shock is not high enough to cause bursting, the membrane develops one or more inner sub-vesicles engulfing the solute composition of the environment (Supplementary Movie 2b, ESI†).

• Discocyte: as the solute concentration ρ_out is decreased, or the ratio ν is increased, the vesicle instead reaches a flat plate-like shape (Supplementary Movie 2c, ESI†).

• Elongated: as the concentration ratio ν is further increased (ν ∼ 0.6), the vesicle transitions to a prolate configuration, closer to an ellipsoid (Supplementary Movie 2d, ESI†).

• Globular: for moderate values of the concentration ratio (0.6 ≲ ν ≲ 0.75) the vesicle is distorted from a normal spherical shape but does not present a clear symmetry break (Supplementary Movie 2e, ESI†).

• Sphere: at high concentration ratio (ν ∼ 1) the vesicle conserves its sphere-like fluctuating shape independently of the solute concentration ρ_out (Supplementary Movie 2f, ESI†).

A few important conclusions can be taken from Fig. 3: (1) we recover the main equilibrium shapes – namely prolates (labeled here as globular or elongated), discocytes and stomatocytes – observed experimentally,^9,15 numerically,^32,33 and predicted theoretically^26–28 (2) we identify a large subset of the parameter space where the osmotic shock leads the vesicle to inner compartmentalization (labeled as inner bud in the figures), absorbing exclusively external solute particle into a smaller sub-vesicle inside the original one. Such a process is very similar to experimentally observed endocytic events that could be relevant to how proto-cells first developed compartmentalization as a form of reaction tanks or proto-organelles in the early stages of life evolution.^14,15 (3) The dependence of the final morphology proves to be highly non-trivial and does not only depend on the concentrations ratio ν, which would lead to a diagram consisting only of horizontal lines as equilibrium curvature models predict^26,27 (see Section V in the ESI†). We thus find that assuming an equilibration of the osmotic shock only via volume adaptation, as previously considered, is not sufficient to properly describe shape transformations in vesicles. In order to unveil the underlying process, we now turn to the dynamics of the system and study how the morphological changes of the vesicle (characterized by its shape, surface area and volume) dynamically relate to the relaxation of the osmotic pressure across the membrane.

3.3 Surface stretching and volume adaptation are coupled to the shock magnitude

Consistent with previous work,^9,15,27,28 we find that, upon osmotic shock, the vesicle deforms by adapting its volume to reduce the pressure on its surface Π. The previous theoretical work assuming constant surface area and a perfect osmotic term cancellation (see Section V in the ESI†) predicts that the equilibrium shape and its reduced volume will depend solely on the ratio of the solute concentrations on the inner and outer side of the membrane, such that the reduced volume is v = ρ_in/ρ_out. To investigate this prediction we measure the reduced area at time τ, a_τ, defined as the ratio of the surface area at time τ to the initial surface area: a_τ = A_τ/A₀. Analogously, the measured reduced volume is defined as v_τ = V_τ/V₀ (see Sections I and II in the ESI† for more details on the computation of these quantities). As shown in Fig. 4a, we find that the relaxed equilibrium volume does not follow exactly the concentration ratio imposed by the osmotic shock, like theoretical considerations based purely on membrane bending energy and constant surface area would predict.^27,28 The reason behind this discrepancy is the fact that a part of the osmotic shock in our model is absorbed by membrane stretching. This is visible in Fig. 4b, where it is apparent that a change in volume is associated with the change in the surface area. Note that, because solute particles freely diffuse in the system, the dynamics is set by the diffusion timescale (τ_D = 100 time steps here), which corresponds to the average time needed for diffusive particles to travel one unit length, and sets the lower limit for the right regime of the system. This means that membrane and solute particles do not interact before this time and explains why no change is observed in Fig. 4b over the initial 100 time steps.

Fig. 4 Measured equilibrium reduced volume and relaxation dynamics of the volume and surface area. (a) Reduced volume measured at the end of the simulation, v_eq, plotted against the theoretically expected value v = ρ_in/ρ_out. Data points are colored according to the outer particle density ρ_out (see legend). The dashed line represents y = x and the dotted line y = 1. Each point corresponds to the equilibrium average over 5 independent runs. Error bars are negligible. (b) Measured reduced surface area a_τ and reduced volume v_τ for three representative cases (see legend). Inner bud: ρ_out = 0.18, ν = 0.31. Discocyte: ρ_out = 0.12, ν = 0.31. Tense sphere: ρ_out = 0.18, ν = 1.25. The simulation time τ on the x-axis is in logarithmic scale and expressed in time steps. The cross and star symbols refer to the corresponding points on the morphology diagram (see Fig. 3).

As Fig. 4a shows, the actual reduced volume does scale linearly with the theoretically expected one for low values (v < 1) but instead remains always above it and saturates around v = 1. For v < 1, the difference between the expected and observed values overall increases as the solute concentration decreases; for v > 1, instead, the measured reduced volume at equilibrium, v_eq, does not exactly collapse to one and actually exceeds it. This observation, together with the membrane area change displayed in Fig. 4b, strongly suggests that a part of the energy provided by the osmotic shock is transformed into area strain, coupling it to the volume evolution in a non-trivial way. This indicates that volume and area relaxation is an intricate process where the equilibrium volume does not only depend on the ratio of inner and outer concentrations (v_eq ≠ v = ρ_in/ρ_out) but rather on the actual magnitude of the osmotic shock.

Membrane stretching is rarely accounted for in theoretical models, which does not allow such models to explore hypotonic conditions. Molecular dynamics simulations can access these conditions as shown by the magenta curves in Fig. 4b and 5, corresponding to a tense sphere: the vesicle maintains its spherical shape but increases its volume via stretching of the membrane. Nevertheless, when assessing the simulation results it is important to consider how the mechanical properties of the simulated membrane compare to those of biologically relevant membranes. In the case simulated here, the compressibility modulus of the membrane is estimated to be around 1 mN m⁻¹ (assuming the membrane bead size σ ∼ 10 nm)³⁴ which is around an order of magnitude lower than, for instance, the three-beads-per-lipid CG lipid membrane model³⁶ and even further from experimentally measured values for single bilayer vesicles of varying composition (ranging from ∼60 to ∼1700 mN m⁻¹).^37,38 However, the rupture tension (tension beyond which the membrane breaks) presents similar disparities, resulting in critical area strains a_c (normalized surface area change beyond which the membrane breaks) that are consistent with the experimental ones. Indeed, for the model used in this work a_c ∼ 0.09³⁴ (very similar to that of the Cooke et al. membrane model³⁶ of a_c ∼ 0.08) while experiments have measured typical lipid bilayer vesicles to have a critical area strain of a_c ∼ 0.02–0.05.^37,38 This suggests that our simulations are somewhat overestimating membrane stretching (9% instead of 5%) but does not invalidate the role of membrane stretching in the response to osmotic shock. In fact, membrane stretching in response to osmotic shock is known to be crucial for the equilibration of osmotic pressures at the cell interface via mechanosensitive channels in unicellular organisms.^39,40

Taken together, the above observations indicate that the absorption of osmotic shocks by volume adaptation is more involved than just a volume equilibration at a constant surface area. Indeed, we find that membrane stretching plays a non-negligible role and absorbs part of the energetic shock. Moreover, the magnitude of the volume change depends on the magnitude of the shock, which in turn is dependent on the solute concentration: for a same concentration ratio, the larger the absolute solute concentration the larger the deformation.

3.4 Osmotic shock adaptation is a two-step process: absorption and dissipation

Next we turn to the dynamics of the relaxation process. We observe that the system displays a striking separation of scales in the relaxation process towards equilibrium, as visible in Fig. 5a. This figure shows the measured osmotic pressure and vesicle volume (see Section I in the ESI†) for three representative cases: (i) a hypertonic osmotic shock resulting in an inner sub-vesicle (green line, indicated with a cross on the diagram in Fig. 3), (ii) a hypertonic shock resulting in a discocyte (purple line, indicated with a star on the diagram in Fig. 3) and (iii) a hypotonic shock resulting in a tense sphere (pink line).

Fig. 5 Shock relaxation dynamics. (a) Measured osmotic pressure difference across the membrane Π (in MD units: [Π] = εσ⁻³) and reduced volume v_τ relaxation curves for three representative cases (see legend). Inner bud: ρ_out = 0.18, ν = 0.31. Discocyte: ρ_out = 0.12, ν = 0.31. Tense sphere: ρ_out = 0.18, ν = 1.25. The green and purple dashed vertical lines indicate the time of maximum bending energy for the “Inner bud” and “Discocyte” scenarios respectively (see panel (b)). (b) Bending energy evolution for the same three representative cases (see legend). The dashed line is the mean bending free energy of a relaxed spherical vesicle of the same size. Inset: Snapshots of the vesicles at different points in time for each scenario. Particles are colored accordingly. Time is in simulation time steps. The cross and star symbols refer to the corresponding points on the morphology diagram (see Fig. 3).

Taken together, these results show the volume and pressure dynamics are coupled and occur on a much shorter time scale than the shape relaxation, since equilibrium for both pressure and volume is reached at times shorter than τ ≪ 10⁵ time steps in all cases. Interestingly, such volume dynamics display striking similarities with the osmotic shock relaxation curves presented in the work by Gabba et al., where the authors have observed an analogous decaying profile in the normalized volume evolution of impermeable vesicles subject to hypertonic osmotic shocks, both in vitro²⁹ and in a physiochemical kinetic model.³⁰ The hypotonic scenario resulting in a tense sphere is of particular interest, as it presents an incomplete osmotic pressure relaxation: Π_eq ∼ 0.25εσ⁻³. This is caused by the vesicle reaching its maximal allowed volume (for maximal stretching of the lipid membrane, a_c ∼ 0.09) before the pressure difference is completely cancelled. However, this remaining osmotic stress is not sufficient to overcome the energy barrier for bursting the vesicle, which results in the observed tense sphere morphology. For hypotonic shocks of larger magnitude we observe vesicle bursting instead of the tense sphere morphology. This is a distinct feature of hypotonic scenarios, indicating that vesicles are inherently better prepared to endure hypertonic shocks. A possible way of adaptation to such conditions is the successive opening and closing of pores to release the tension imposed by the solution, which has been observed experimentally.²⁴ However, our simulations do not recover such cyclic dynamics, for which we believe a more detailed membrane model would be necessary.

Tracing the behaviour of the mechanical energy of the vesicle during the relaxation process, measured by the reduced bending free energy F̄ in time (refer to Section III in the ESI† for a detailed definition and computation details), allows us to gain further insight into the relaxation process (Fig. 5b). For the discocyte and inner bud scenarios (purple and green curves, respectively), we observe a sharp maximum in the bending energy which is concurrent with the volume equilibration (τ ∼ 1.2 × 10⁴ time steps and τ ∼ 4.8 × 10⁴ time steps respectively, as shown by the two vertical lines on Fig. 5a). This indicates that the osmotic shock absorption induces very strong deformations on the vesicle which are then smoothed out over time as the vesicle relaxes to its minimum energy. Notice that for the green curve in Fig. 5b this implies an inward budding event which results in a sudden drop of the energy. For the tense sphere the vesicle does not change its shape and remains spherical throughout the simulation, hence the bending energy fluctuates around F̄ = 1. Refer to Supplementary Movies 2b, 2c and 2g (ESI†) for a better visualization of the relaxation process.

In the case of the vesicle that relaxes to the discocyte conformation, the equilibrium value of F̄ ∼ 1.8 agrees very well with the value previously predicted by the analytical models that consider only vesicle curvature.^27,28 However for the conditions when the inner bud is formed we observe a three stage evolution: first a large amount of energy is absorbed by developing a very strong deformation (peak of F̄ ∼ 3), which then relaxes as F̄ equilibrates around F̄ ∼ 2 by developing a large invagination that is still connected to the parent membrane (free energy value predicted by curvature models for a stomatocyte conformation^27,28). Contrary to what curvature models can predict, the vesicle does not remain in this local minimum and instead evolves further towards F̄ < 1 by developing an inner compartment as the invagination buds off. This inner bud scenario agrees very well with experimentally observed vesicle budding under osmotic stress.^9,15 Such transformations are however inaccessible to equilibrium curvature models where vesicles are not allowed to break.

The separation of scales described in this section can thus be understood as an initial volume adaptation which absorbs the osmotic shock and a subsequent shape transformation which dissipates the accumulated tension. The first step therefore corresponds to solvent flux across the membrane while the second is achieved by lipid rearrangement. In our simulations the solvent is modelled implicitly (solute particles follow Brownian dynamics), hence the first relaxation is diffusion limited while the second one is controlled by the membrane physical properties (i.e. bending rigidity, fluidity, spontaneous curvature, compressibility modulus).

4 Conclusions

In this work we have adopted a novel approach in the study of vesicle osmotic shock adaptation by combining molecular dynamics with a coarse-grained lipid membrane model³⁴ and explicit solute particles. This allows to overcome some of the limitations of previous works while still recovering the main known features of this physical process.

In particular, here we demonstrate that with this approach we have good control over the magnitude of the osmotic shock vesicles are subject to, which scales appropriately with solute concentrations inside and outside the lipid body (see Fig. 2). Moreover, we show that this method allows to investigate known equilibrium vesicle morphologies and how they relate to the osmotic perturbation. Indeed, not only do we recover the main vesicle morphologies predicted^27,32 and observed experimentally⁹ in previous works but we can also precisely map them to the relevant parameters of the system (see Fig. 3).

Furthermore, using MD simulations allows for membrane scission, bursting and budding events, typically inaccessible to continuum modeling approaches. This enables us to investigate biologically relevant topologies such as the inner bud (or vesicle-inside-vesicle), where the perturbed vesicle develops an inner compartment engulfing external solute to accommodate the osmotic shock. In this work we thus define a subset of the parameter space for which endocytosis occurs as a result of osmotic shocks (see Fig. 3), possibly an important element of proto-cell development in the early stages of life.^9,14,15 Likewise, working with a coarse grained MD model allows us also to explore hypotonic shocks, which result in tense spheres with increased surface areas.

Another important feature of this approach is that it allows tracking of the dynamics of the relaxation process. While molecular dynamics have already been considered in previous studies of vesicular osmotic responses,³³ the present work offers important advantages in that the relaxation of the osmotic shock is purely the result of the molecular interactions while in previous works it was imposed externally. We thus avoid any assumptions regarding the relaxation dynamics and work with membrane transformations defined only by the imposed osmotic shock. This allows us to dynamically characterize the transformation process by following the evolution of relevant physical magnitudes such as the bending free energy, reduced volume or pressure differential across the lipid bilayer (see Fig. 5). As a result, we find that the vesicle evolution upon osmotic upshift is a two-step process where the membrane first absorbs the shock by deforming away form its spherical shape (volume and surface area changes), which results in a sharp increase of the bending energy, then dissipated by reaching the observed equilibrium morphology. Overall, the dynamic insight thus obtained provides valuable information complementing the shortcomings of classical approaches to the problem.

In the present work we have restricted our results to model membranes displaying physiologically relevant physical properties in simple scenarios. However, the combination of a coarse-grained membrane model with explicit solute particles, which allows to precisely control different aspects of the membrane physics, could easily be applied to systems of growing complexity (e.g. membranes with non-zero spontaneous curvature or in different fluid phases). Possible further applications of the approach taken here could thus be to study multiple endocytic events in fluctuating and heterogeneous environments, which could give rise to complex internal vesicle structures of different compositions. The model is also well suited for studying vesicle bursting cycles observed in hypotonic environments²⁴ or the different morphologies observed for multiphasic⁸ and protein charged vesicles,^11,19 where the dynamical nature of MD simulations should prove useful.

Conflicts of interest

There are no conflicts to declare.

Note added after first publication

This article replaces the version published on 25th February 2021, which contained typesetting errors in Fig. 4 and 5 in the PDF version. No changes have been made to the HTML version.

Acknowledgements

We acknowledge support from the Royal Society (C. V. C. and A. Š.), the Medical Research Council (C. V. C. and A. Š.), and the European Research Council (Starting grant “NEPA” 802960 to A. Š.). We thank Johannes Krausser and Ivan Palaia for fruitful discussions.

Notes and references

1. J. H. Hurley, E. Boura, L. A. Carlson and B. Róycki, Membrane budding, Cell, 2010, 143(6), 875–887.
2. T. H. Cheffings, N. J. Burroughs and M. K. Balasubramanian, Actomyosin Ring Formation and Tension Generation in Eukaryotic Cytokinesis, Curr. Biol., 2016, 26(15), R719–R737,  DOI:10.1016/j.cub.2016.06.071.
3. L. Harker-Kirschneck, B. Baum and A. ela Šarić, Changes in ESCRT III filament geometry drive membrane remodelling and fission in silico, BMC Biol., 2019, 17(1), 1–8.
4. G. T. Risa, F. Hurtig, S. Bray, A. E. Hafner, L. Harker-Kirschneck and P. Faull, et al., Proteasome-mediated protein degradation resets the cell division cycle and triggers ESCRT-III mediated cytokinesis in an archaeon, bioRxiv, 2019, 774273 . Available from: https://www.biorxiv.org/content/10.1101/774273v1.
5. Y. Liao, S. Ithurbide, J. Löwe and I. G. Duggin, Two FtsZ proteins orchestrate archaeal cell division through distinct functions in ring assembly and constriction, bioRxiv, 2020, 133736 . Available from: https://www.biorxiv.org/content/10.1101/2020.06.04.133736v1.abstract?{%}3Fcollection=.
6. T. F. Zhu and J. W. Szostak, Coupled growth and division of model protocell membranes, J. Am. Chem. Soc., 2009, 131(15), 5705–5713.
7. M. Andes-Koback and C. D. Keating, Complete budding and asymmetric division of primitive model cells to produce daughter vesicles with different interior and membrane compositions, J. Am. Chem. Soc., 2011, 133(24), 9545–9555.
8. Y. Dreher, J. P. Spatz and K. Göpfrich, Protein-free division of giant unilamellar vesicles controlled by enzymatic activity, bioRxiv, 2019, 881557 . Available from: http://biorxiv.org/content/early/2019/12/30/2019.12.30.881557.
9. W. Zong, Q. Li, X. Zhang and X. Han, Deformation of giant unilamellar vesicles under osmotic stress, Colloids Surf., B, 2018, 172, 459–463 . Available from: http://www.sciencedirect.com/science/article/pii/S092777651830585X.
10. W. T. Snead, C. C. Hayden, A. K. Gadok, C. Zhao, E. M. Lafer, P. Rangamani and J. C. Stachowiak, Membrane fission by protein crowding, Proc. Natl. Acad. Sci. U. S. A., 2017, 114(16), E3258–E3267.
11. J. Steinkühler, R. L. Knorr, Z. Zhao, T. Bhatia, S. M. Bartelt and S. Wegner, et al., Controlled division of cell-sized vesicles by low densities of membrane-bound proteins, Nat. Commun., 2020, 11(1), 1–20.
12. M. Lluch-Senar, E. Querol and J. Piñol, Cell division in a minimal bacterium in the absence of ftsZ, Mol. Microbiol., 2010, 78(2), 278–289.
13. J. W. Szostak, The Narrow Road to the Deep Past: In Search of the Chemistry of the Origin of Life, Angew. Chem., Int. Ed., 2017, 56(37), 11037–11043.
14. S. Svetina, Vesicle budding and the origin of cellular life, ChemPhysChem, 2009, 10(16), 2769–2776.
15. W. Zong, S. Ma, X. Zhang, X. Wang, Q. Li and X. Han, A Fissionable Artificial Eukaryote-like Cell Model, J. Am. Chem. Soc., 2017, 139(29), 9955–9960,  DOI:10.1021/jacs.7b04009 , PMID: 28677973.
16. R. Mercier, P. Domínguez-Cuevas and J. Errington, Crucial Role for Membrane Fluidity in Proliferation of Primitive Cells, Cell Rep., 2012, 1(5), 417–423.
17. T. Baumgart, S. T. Hess and W. W. Webb, Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension, Nature, 2003, 425(6960), 821–824.
18. T. Baumgart, S. Das, W. W. Webb and J. T. Jenkins, Membrane elasticity in giant vesicles with fluid phase coexistence, Biophys. J., 2005, 89(2), 1067–1080.
19. J. C. Stachowiak, C. C. Hayden and D. Y. Sasaki, Steric confinement of proteins on lipid membranes can drive curvature and tubulation, Proc. Natl. Acad. Sci. U. S. A., 2010, 107(17), 7781–7786.
20. M. Ohno, T. Hamada, K. Takiguchi and M. Homma, Dynamic behavior of giant liposomes at desired osmotic pressures, Langmuir, 2009, 25(19), 11680–11685.
21. M. Saleem, S. Morlot, A. Hohendahl, J. Manzi, M. Lenz and A. Roux, A balance between membrane elasticity and polymerization energy sets the shape of spherical clathrin coats, Nat. Commun., 2015, 6, 6249.
22. A.-L. Bernard, M.-A. Guedeau-Boudeville, L. Jullien and J.-M. di Meglio, Raspberry vesicles, Biochim. Biophys. Acta, Biomembr., 2002, 1567, 1–5 . Available from: http://www.sciencedirect.com/science/article/pii/S000527360200617X.
23. E. Boroske, M. Elwenspoek and W. Helfrich, Osmotic shrinkage of giant egg-lecithin vesicles, Biophys. J., 1981, 34(1), 95–109.
24. M. Chabanon, J. C. Ho, B. Liedberg, A. N. Parikh and P. Rangamani, Pulsatile Lipid Vesicles under Osmotic Stress, Biophys. J., 2017, 112(8), 1682–1691,  DOI:10.1016/j.bpj.2017.03.018.
25. K. Oglecka, P. Rangamani, B. Liedberg, R. S. Kraut and A. N. Parikh, Oscillatory phase separation in giant lipid vesicles induced by transmembrane osmotic differentials, eLife, 2014, 3, 1–18.
26. W. Helfrich, Elastic Properties of Lipid Bilayers: Theory and Possible Experiments, Zeitschrift für Naturforschung C, 1973, 28, 11–12.
27. U. Seifert, Configurations of fluid membranes and vesicles, 1997, vol. 46, pp. 13–137.
28. U. Seifert, K. Berndl and R. Lipowsky, Shape transformations of vesicles: Phase diagram for spontaneous-curvature and bilayer-coupling models, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 1182–1202,  DOI:10.1103/PhysRevA.44.1182.
29. M. Gabba, J. Frallicciardi, J. van't Klooster, R. Henderson, Ł. Syga and R. Mans, Weak Acid Permeation in Synthetic Lipid Vesicles and Across the Yeast Plasma Membrane, Biophys. J., 2020, 118(2), 422–434.
30. M. Gabba and B. Poolman, Physiochemical Modeling of Vesicle Dynamics upon Osmotic Upshift, Biophys. J., 2020, 118(2), 435–447,  DOI:10.1016/j.bpj.2019.11.3383.
31. T. Ruiz-Herrero, T. G. Fai and L. Mahadevan, Dynamics of Growth and Form in Prebiotic Vesicles, Phys. Rev. Lett., 2019, 123(3), 38102,  DOI:10.1103/PhysRevLett.123.038102.
32. G. Gueguen, N. Destainville and M. Manghi, Fluctuation tension and shape transition of vesicles: Renormalisation calculations and Monte Carlo simulations, Soft Matter, 2017, 13, 6100–6117.
33. H. Yuan, C. Huang and S. Zhang, Dynamic shape transformations of fluid vesicles, Soft Matter, 2010, 6, 4571–4579.
34. H. Yuan, C. Huang, J. Li, G. Lykotrafitis and S. Zhang, Oneparticle-thick, solvent-free, coarse-grained model for biological and biomimetic fluid membranes, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011905,  DOI:10.1103/PhysRevE.82.011905.
35. http://lammps.sandia.gov .
36. I. R. Cooke, K. Kremer and M. Deserno, Tunable generic model for fluid bilayer membranes, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 2–5.
37. R. Kwok and E. Evans, Thermoelasticity of large lecithin bilayer vesicles, Biophys. J., 1981, 35, 637–652.
38. D. Needham and R. S. Nunn, Elastic deformation and failure of lipid bilayer membranes containing cholesterol, Biophys. J., 1990, 58, 997–1009.
39. E. S. Haswell, R. Phillips and D. C. Rees, Mechanosensitive channels: what can they do and how do they do it?, Structure, 2011, 19, 1356–1369.
40. A. Paraschiv, S. Hegde, R. Ganti, T. Pilizota and A. Šarić, Dynamic clustering regulates activity of mechanosensitive membrane channels, Phys. Rev. Lett., 2020, 124, 048102.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm02012e

---