Kinetics of a Multilamellar Lipid Vesicle Ripening - ACS Publications

Feb 16, 2016 - vesicle (MLV) is effectively explored using a mesoscale coarse-grained ... and there is not yet a kinetics theory suitable to describe ...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCB

Kinetics of a Multilamellar Lipid Vesicle Ripening: Simulation and Theory Rui Xu†,§ and Xuehao He*,‡ †

Department of Polymer Science and Engineering, School of Chemical Engineering and Technology, Tianjin University, 300072 Tianjin, China ‡ Department of Chemistry, School of Science, Tianjin University, and Collaborative Innovation Center of Chemical Science and Engineering (Tianjin), 300072 Tianjin, China § State Key Laboratory of Separation Membranes and Membrane Processes, School of Material Science and Engineering, Tianjin Polytechnic University, Tianjin 300387, China S Supporting Information *

ABSTRACT: Lipid vesicle ripening via unimolecular diffusion and exchange greatly influences the evolution of complex vesicle structure. However, this behavior is difficult to capture using conventional experimental technology and molecular simulation. In the present work, the ripening of a multilamellar lipid vesicle (MLV) is effectively explored using a mesoscale coarse-grained molecular model. The simulation reveals that a small MLV evolves into a unilamellar vesicle over a very long time period. In this process, only the outermost bilayer inflates, and the inner bilayers shrink. With increasing MLV size, the ripening process becomes complex and depends on competition between a series of adjacent bilayers in the MLV. To understand the diffusion behavior of the unimolecule, the potentials of mean force (PMFs) of a single lipid molecule across unilamellar vesicles with different sizes are calculated. It is found that the PMF of lipid dissociation from the inner layer is different than that of the outer layer, and the dissociation energy barrier sensitively depends on the curvature of the bilayer. A kinetics theoretical model of MLV ripening that considers the lipid dissociation energy for curved bilayers is proposed. The model successfully interprets the MLV ripening process with various numbers of bilayers and shows potential to predict the ripening kinetics of complex lipid vesicles.

■

INTRODUCTION Vesicle ripening is a fundamental problem for understanding the stabilities of micellar structures. Vesicle ripening1,2 and fusion3,4 are two major vesicle destabilization mechanisms. Usually, the fusion of lipid or amphiphilic polymer vesicles rarely proceeds because the process must overcome a higher energy barrier.5 By contrast, vesicle ripening happens more easily through unimolecular diffusion that avoids the high energy barriers associated with vesicle fusion.6,7 The ripening of unilamellar vesicles (ULVs)8−18,12−14 has been confirmed through experiments applying fluorescent8,12,15,19−21 and deuterium13,22 labeling technologies. In the ripening process, surfactants or lipids in ULVs exchange through a common solution environment. For lipid vesicles, the lipid molecule motions include transbilayer (flip-flop) movement of the lipid between the inner and outer leaflets in a bilayer vesicle, desorption of the lipid from the outer layer of a donor vesicle, and binding of the lipid on the outer layer of an acceptor vesicle. Depending on the properties of the surfactants and lipids, different ripening products with unimodal or bimodal size distributions can be obtained at the final equilibrium state.1,2,23,24 Compared with the structure of ULVs, multilamellar vesicles (MLVs) possess an onion-like structure. An MLV can be regarded as a set of ULVs which are embedded one within the © XXXX American Chemical Society

next, from small to large. As complex vesicles, MLVs have special characteristics, such as a high efficiency in protection and entrapment of hydrophobic agents.25,26 MLVs have been used as microreactors,27−29 T2-MRI contrast agents,30 and carriers for controlled release of drugs,25,26,31,32 genes,33−35 or enzymes.36 The shape and stability of MLVs, depending on temperature,37,38 pressure,39,40 photo irradiation,41 shear rheology,42,43 and solvent environment, such as ionic strength,38 pH,44 and solvent polarity,45 have been studied. MLV ripening has different mechanisms than that of ULVs. The lipids in MLVs diffuse and exchange between the contiguous bilayers. Understanding the ripening of MLVs is helpful to effectively tune the structure and functionality of complex vesicles. Molecular simulation is a powerful tool to achieve dynamic information at the atomistic scale and understand the selfassembly process of lipids. Unfortunately, conventional atomistic molecular dynamics simulations are generally limited to hundreds of nanoseconds time scales and fail to reproduce the ripening phenomena, especially for large micellar system. It is still challenging to simulate vesicle ripening even when using Received: December 13, 2015 Revised: February 9, 2016

A

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ⎡ 1 ⎛ d ⎞6 12 ⎛ d ⎞0.5⎤ TT ⎜ m⎟ ⎥ Unonbonded = ε⎢ ⎜ m ⎟ − ⎢⎣ 11 ⎝ d ⎠ 11 ⎝ d ⎠ ⎥⎦

coarse-grained models with explicit solvent molecules, such as the Martini model.46 The large-scale coarse-grained models with implicit solvation47−49 are needed for the mesoscale simulation of vesicle ripening. According to our knowledge, however, simulations of MLV ripening are not yet reported, and there is not yet a kinetics theory suitable to describe MLV ripening. In this work, the ripening of MLVs is first explored by using an implicit solvent coarse-grained lipid model, which was recently developed by our group.47 To understand the dependence of lipid diffusion on the curvature of bilayers, the potentials of mean force (PMFs) of a single lipid across bilayer vesicles with different sizes are calculated via the umbrella sampling method. According to the PMF results, we propose a theoretical kinetics model to describe the MLV ripening process in which the activation energy of lipid dissociation is explicitly introduced. This kinetics model successfully reproduces the ripening of MLVs and predicts the evolution trajectory of each bilayer. It provides a theoretical framework for the exploration of ripening of complex vesicles and supports a deeper understanding of the microstructural formation and the evolution of complex lipid vesicles.

HH/HT Unonbonded =

6 2 ⎛ dm ⎞ ε⎜ ⎟ 55 ⎝ d ⎠

(1)

(2)

Here, ε is the depth of the potential well, and dm is the distance at which the potential reaches its minimum and the force becomes zero. The two parameters (ε and dm are 2.1 kBT and 1.2 nm, respectively) control the bending rigidity and diffusivity of the lipid membrane. The model successfully reproduces typical phase behaviors and mechanical properties of a DPPC bilayer. In this model, the length unit is σ0 = 1 nm, the mass unit is m0 = 1 u (1/12 the mass of a 12C atom), and the energy unit is ε0 = kBT ≈ 2.62 kJ/mol at 315 K. The beads of headgroup and tails have the same mass, 245 u. The equilibrium area per lipid using our model is 0.648 nm2,47 close to 0.64 nm2 in experiment50 and MARTINI model.51 The bending rigidity of the model is 4.35 × 10−20 J,47 which is in agreement with the experimental value of 5.0 × 10−20 J.52 The stretching modulus measured in the model is 157.5 mN/m,47 smaller than the experimental value for DPPC 231 ± 20 mN/ m.50 The edge line tension is 7.2 pN47 in the range of experimental results of 6.5−30 pN for lipid bilayers.53−57 The lateral diffusion coefficient of 3.925 × 10−6 cm2 s−1 is faster by 2 orders of magnitude compared with the experimental value of 1.7−7.8 × 10−8 cm2 s−1 and the AA simulation of 2.6−3.9 × 10−8 cm2 s−1.47,58 The rate of flip-flop calculated in our model is 6.5 × 1011 s−1, which is 15 orders of magnitude faster than the value of 4.2 × 10−4 s−1 measured in experiments59,60 and the same quantitative order with the Cooke−Deserno model.48 Using the present coarse-grained model, the lipids flip-flop very easily due to the short chain lengths and weak head−tail repulsion. It provides a possibility to study the formation and ripening process of complex lipid vesicles via molecular simulation. Langevin dynamics are used in all simulations. The time step is set to 0.15 ps, and the friction constant in the Langevin thermostat is 0.5 ps. The simulation temperature is set to 315 K. To avoid artificiality at the cutoff position, a standard shift function (cubic spline function) is added to the nonbonded potential function with rshift = 2.25 nm, with a cutoff radius rcutoff of 2.772 nm.61 Past work has demonstrated that such a coarsegrained lipid model is highly effective in simulating the micellization of lipids and that the acceleration ratio is approximately 5−6 orders of magnitude compared with that of the conventional atomistic model.47 Because the model uses an implicit solvent force field, the diffusion of lipids in solution is very fast and suitable for studying micellar ripening, which is related to lipid diffusion and exchange. Detailed information on the CG model and parameter optimization can be found in our prior work.47 To study MLV ripening, the MLV is first prepared from a lipid droplet. The initial droplet is composed of randomly oriented lipids that form MLVs via in situ micellization in a short time. The number of bilayers of the MLV can be tuned by varying the initial size of the lipid droplet. After MLV formation, common molecular dynamics simulations are carried out in the NVT ensemble to explore the ripening process. All simulations are carried out using in-house software (Magic MD) with GPU acceleration, and visualizations are performed with VMD.62

■

SIMULATION AND RIPENING KINETICS THEORY Coarse-Grained Lipid Model and Simulation. In the present work, we use an implicit solvent coarse-grained lipid model, which we recently developed.47 In this model, each lipid molecule (dipalmitoyl-phosphatidylcholine, DPPC) in solution is mapped onto three coarse-grained beads, i.e., one hydrophilic head (H) and two hydrophobic tails (Ts) (see DPPC model in Figure 1). The force field of the coarse-grained lipid is

Figure 1. Schematic drawing of the lipid molecule mapping from the all-atom model (left) into the coarse-grained model (right). A lipid molecule in water is coarse-grained as three beads: a head bead (H) and two tail beads (T).

optimized according to experimental data which includes the phase behavior and mechanical properties of DPPC lipid bilayers. The bonded potentials are treated as harmonic potentials with Ubond = 0.5kbond(d − d0)2 and Uangle = 0.5kangle(θ − θ0)2, respectively. kbond and kangle are the force constants; d0 is the equilibrium bond length, and θ0 is the equilibrium angle. The parameters of bonded potentials are kbond = 1000 kJ mol−1 nm−2, d0 = 1 nm, kangle = 20 kJ mol−1 rad−2, and θ0 = 30°. The nonbonded interaction Unonbonded uses a variation of the Lennard-Jones potential form to incorporate long-range attractive interactions. B

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B dNi ,in

Ripening Kinetics Theory. A ripening kinetics theory is built for an MLV according to the lipid migration processes. The flip-flop and exchange motions of the lipids in the ith bilayer in MLV are schematically shown in Figure 2. Here, the

dt

= ki ,out → inNi ,out − ki ,in → outNi ,in + kion ,inNi ,in

Niaq V iaq

− kioff ,inNi ,in

(4)

Vaq 0

In the system containing a single MLV with n bilayers, is the volume of the aqueous phase outside the MLV, and its value is equal to the volume of the simulation box minus that of 3 the MLV (Vaq 0 = V − 4πr1/3 where V is the simulated space and r1 is the mean radius of the outermost MLV bilayer). Vaq n is the aqueous phase volume in the innermost cavity of the MLV (Vaq n = 4πr3n/3 where rn is the mean radius of the innermost bilayer of 3 3 MLV). The volume of solvent Vaq i equals 4π(ri − ri−1)/3, where ri and ri−1 are the radius of the ith bilayer and the (i − 1)th bilayer. Compared with the exchange motion of lipids, the flipflop of lipids between the hemileaflets of a vesicular bilayer is a fast process. For a large MLV, the thickness of the bilayer can be ignored, and the numbers of lipids in the inner and outer hemileaflets of a bilayer are equal, i.e., Ni,in = Ni,out = N∥i /2. N∥i is the number of lipids in the ith bilayer. The bilayer area is 4πr2 = N∥a/2, where r is the mean radius of a bilayer in MLV and a is the area per lipid. As a result, the rate equation (eq 5) of the ith bilayer for MLV is achieved.

Figure 2. Scheme of ripening process of the ith bilayer in a multilamellar lipid vesicle. It is related with the motions of lipids: flipflop of lipids between hemileaflets of the bilayers and exchange of lipids with the adjacent solution. The numbers of the outer and inner hemileaflets of the ith bilayer are Ni,out and Ni,in, and the numbers of free lipids in the aqueous phase adjacent the ith bilayer are Naq i−1 and Naq i . ki,out→in and ki,in→out are the flip-flop rate constants between on hemileaflets of the vesicle. kon i,out and ki,in are the binding rate constants from aqueous phase to the outer and inner hemileaflets of the ith off bilayer, and koff i,out and ki,in are the dissociation rate constants from the outer and inner hemileaflets of the ith bilayer to the adjacent aqueous phase, respectively.

aq ⎞ dNi Niaq 1⎛ −1 on Ni off off = ⎜kion ,out aq + k i ,in aq − k i ,out − k i ,in ⎟Ni dt 2⎝ Vi−1 Vi ⎠

Reasonably, the evolution of the number of lipids in the ith aqueous phase satisfies the following equation. aq ⎞ dNiaq 1⎛ on Ni = ⎜kioff − k ,in i ,in aq ⎟Ni dt 2⎝ Vi ⎠

shape of every bilayer is assumed to be spherical. In our theoretical framework, the coalescence mechanism between two bilayers is neglected for the ripening process. The (i, out) single layer is the outer layer and the (i, in) layer is the inner layer of the ith bilayer. The flip-flop of lipids takes place between two single layers of a bilayer. The exchange of lipids from a single layer to its adjacent aqueous phase can be considered as an interfacial reaction that depends on the area of the layer. For convenience, the size change of the layer is described by the number of lipids. The concentration of lipids solved in water is equal to the number density of free lipids in the divided space between two bilayers. So, the evolution of the number of lipids in the outer (i, out) leaflet of the ith bilayer with time can be written as dNi ,out dt

= ki ,in → outNi ,in − ki ,out → inNi ,out + kion ,outNi ,out − kioff ,outNi ,out

(5)

+

Niaq ⎞ 1 ⎛ off ⎜ki + 1,out − kion + 1,out aq ⎟Ni + 1 2⎝ Vi ⎠

(6)

The outermost (i = 0) and innermost (i = n) aqueous phases satisfy the following equations. aq ⎞ dN0aq 1 ⎛ off on N0 = ⎜k1,out − k1,out ⎟N1 dt 2⎝ V 0aq ⎠

(7)

dNnaq N aq ⎞ 1⎛ = ⎜knoff,in − knon,in naq ⎟Nn dt 2⎝ Vn ⎠

(8)

The detailed derivation of the rate equations above can be found in the Supporting Information. Umbrella Sampling. In the ripening kinetics model, the evolution of the number of lipids depends on the rate constants on off off kon i,out, ki,in, ki,out, and ki,in, which are related to the diffusion of lipid molecules between bilayers. To understand the dependence of the rate constants on the bilayer curvature, the potential of mean force (PMF) for a lipid across the vesicle membrane is extracted using umbrella sampling technology,63,64 which can overcome the energy barrier in the complex energy landscape. The umbrella potential acts on the center of mass of the free lipid molecule with a harmonic potential, and the force constant is 1000 kJ mol−1 nm−2. The free lipid is pulled across the vesicle bilayer from the inside out. The distances from the starting point and the destination to the bilayer are equal and are ∼6 nm. The sampling windows are created by pulling the single lipid molecule to different window locations, which use the umbrella potentials with a force constant of 1000 kJ mol−1 nm−2 in a 1.5 ns simulation. The window size is 0.2 nm. Then,

Niaq −1 V iaq −1 (3)

where Ni,out denotes the number of lipids in the outer leaflet of aq the ith bilayer. Vaq i−1 and Ni−1 are the space volume and the number of free lipids in the (i − 1)th aqueous phase, respectively. ki,in→out is the flip-flop rate constant of lipids from the inner layer to the outer layer in the ith bilayer, and the inverse is ki,out→in. koff i,out is the dissociation rate constant of lipids from the outer layer to the aqueous solution, and kon i,out is the binding rate constant of free lipids from the solution to the outer layer. Similarly, the evolution of the number of lipids in the inner leaflet of the ith bilayer is expressed as C

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B each window is equilibrated for 30 ns, followed by a 3 μs production simulation. Finally, the PMF profile is constructed from the biased distributions of the center of mass of the lipids, using the weighted histogram analysis method65 with a relative tolerance of 10−6. We calculate the PMFs of lipids across bilayers with different curvatures of 0.20, 0.18, and 0.16 nm−1. The corresponding lipid vesicles are composed of 1000, 1200, and 1600 lipids, respectively. It is noted that the small vesicles were used for calculating the PMFs to avoid the fluctuation errors inherent to large vesicle membranes. In all simulations the number density of lipids is fixed at 0.05 lipids per cubic nanometer. The umbrella sampling is carried out in the NVT ensemble using the GROMACS 4.5.5 software package.61

■

RESULTS AND DISCUSSION The initial lipid droplet is constructed by randomly packing lipids into a spherical space using per lipid volume of 1.5 nm3. The spherical space determines the size of lipid droplet. The number of bilayers of MLVs is tuned by the size of the initial lipid droplet. The large lipid droplet as an initial state can selfassemble spontaneously into MLVs in a short time. When the size of the initial droplet is 8000, the number of bilayers of MLV equals n = 2 (see Figure 3). When the droplet size is 2 × 104, an MLV with three bilayers is obtained (see Figure 4). As the droplet size increased to 5 × 104, n = 4 is found (see Figure 5). Figure 3a shows the evolution process of the transformation of a MLV (n = 2) into an ULV (n = 1). There are four stages found: MLV formation, MLV ripening, merging of the innermost bilayer, and ULV formation. The corresponding

Figure 4. (a) Structural evolution of a lipid droplet with 2 × 104 lipids. Each snapshot represents the cross-sectional structures at different times. The head and tail of the lipid molecule are blue and yellow, respectively. (b) Variations in the number of lipids in the 1st, 2nd, and 3rd bilayers of the MLV with time. The sudden changes in lipid numbers with time indicate coalescence of the innermost bilayer into the second-innermost bilayer.

changes of the number of lipids in each bilayer are shown in Figure 3b. In the first stage, the lipid molecules in the droplet rapidly rearrange to form an MLV structure with a small spherical micelle in the center (see Figure 3a from 0 to 0.03 μs) due to the repulsive interaction between hydrophobic tails and hydrophilic headgroups or solvent. Then, the inner spherical micelle merges to the second bilayer of MLV, and a two-bilayer vesicle is formed at 0.25 μs. The corresponding number of lipids in bilayers shows a rapid change at stage I (Figure 3b). After MLV formation, MLV ripening proceeds at stage II from 0.25 to 9 μs. At this stage, the driving force for MLV ripening comes from the difference of curvature energy between the two bilayers. The first bilayer grows and the second bilayer gradually shrinks over time as shown in Figure 3b from 0.25 to 3 μs, which makes the distance between these two bilayers increase. It is about 8 nm at 3 μs. As the inner vesicle becomes smaller and smaller, it can move more easily inside the cavity of the outer vesicle. As a result, it collides with the inner layer of the outer vesicle occasionally (see Figure 3a at 9 μs) which accelerates the transfer of lipids from the inner vesicle to the outer vesicle. When the number of lipids in the small vesicle decreases to about 600, the vesicle structure is difficult to maintain. The small vesicle has a very high interfacial energy due to its high curvature and can easily merge with the outer vesicle. It is shown that the numbers of lipids of the two

Figure 3. (a) Structural evolution of a lipid droplet with 8000 lipids. Each snapshot represents the cross-sectional structures at different times. The head and tail of the lipid molecule are colored blue and yellow, respectively. (b) Variations in the number of lipids in the 1st and 2nd bilayers of the MLV with time. D

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

kinetics model is proposed in the Ripening Kinetics Theory section. In this ripening kinetics model, the rate constant is critical to determine the basic rules of MLV ripening. First, we calculate the potential of mean forces (PMFs) of a lipid across different curved bilayers to understand the diffusion behavior of lipids in curved bilayers. Figure 6 shows how the PMFs depend on the

Figure 5. (a) Structural evolution of a lipid droplet with 5 × 104 lipids. Each snapshot represents the cross-sectional structures at different times. The head and tail of the lipid molecule are blue and yellow, respectively. (b) Variation in the number of lipids in each bilayer with time. The sudden changes of lipid numbers with time indicate coalescence of the innermost bilayer into the second-innermost bilayer. Figure 6. Potential of mean forces for a lipid transferring across (a) a flat bilayer with a curvature of 0 nm−1 (512 lipids) and (b) vesicles with curvatures of 0.20 nm−1 (1000 lipids), 0.18 nm−1 (1200 lipids), and 0.16 nm−1 (1600 lipids) are calculated from adaptive biasing force after umbrella sampling. The value of x axis is the distance of the migrating lipid (a) from the middle plane of the bilayer and (b) from the mass center of the vesicle. The starting location of the PMF profile at the left is chosen at the first inflection point in panel b.

bilayers change suddenly at about 12 μs and an ULV is formed finally at 13.5 μs. Figure 4 shows the evolution of a three-bilayer MLV composed of 2 × 104 lipids. The in situ micellization of the large droplet is similar to that of the small droplet with 8000 lipids. A three-bilayer vesicle is formed after 0.14 μs (see Figure 4a). When the innermost vesicle is small enough, it merges into the second bilayer, and then a larger two-bilayer vesicle is formed (see Figure 4a at 3.6 μs). Due to computational limitations of simulation time, the vesicle still maintains a twobilayer structure at the conclusion of 34 μs. The four-bilayer MLV with 5 × 104 lipids also shows the similar ripening phenomena (Figure 5a). The fourth bilayer disappears after about 6.3 μs. During its ripening process, the outermost vesicle becomes larger and larger, and the inner bilayers shrink over time (Figure 5b). It is noteworthy that only the outermost bilayer inflates and all other inner bilayers shrink together for the small MLV. The ripening phenomenon is a very slow process which proceeds through the flip-flop and exchange motions of lipids from the inner bilayers to the outermost bilayer. Usually, it is difficult to observe and analyze the ripening process even from experiment. So far, only the slow diffusion of lipids in the ripening processes of ULVs has been reported in experiments.1,2,66 To understand the ripening mechanism of MLV, a theoretical

bilayer curvature. It is found that lipid dissociation from the bilayer must overcome an energy barrier. This consistent with the results of McLean and Phillips.14 The rate-limiting step for the overall exchange process between ULVs is desorption in the absence of protein. In Figure 6, it is shown especially that the energy barrier for lipid dissociation from the inner layer of a bilayer is larger than that from the outer layer. As the bilayer curvature increases, the energy barriers for lipid dissociation decrease. This indicates that lipids dissociate from bilayers with high curvature more easily. According to the calculated PMF off results, the rate constants koff i,out and ki,in can be expressed as A exp[−ΔE/RT], where A is the pre-exponential factor, ΔE is the activation energy, R is the universal gas constant, and T is off temperature. The rate constants koff i,out and ki,in depend on the curvature of the bilayer, because the dissociated energy barrier of the lipids increases with decreasing bilayer curvature. By E

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B on contrast, kon i,out and ki,in are constants for the binding process with no energy barrier. Furthermore, the energy barrier for the lipid flip-flop motion is smaller than the dissociation energy calculated via the coarse-grained lipid model without solvent. This indicates that neglecting the faster lipid flip-flop motion process in eqs 5−8 is reasonable. However, using our model, the lipid dissociation energy from a planar bilayer is 22.5 kJ mol−1, and the lipid flip-flop energy barrier is 2.5 kJ mol−1, which are both less than 75−80 kJ mol−1 calculated from allatom simulations.64 The highly coarse-grained model without solvent may sacrifice some authenticity, but the lipids in the model have higher diffusion and flip-flop rates which provide the possibility to effectively reveal the vesicle ripening processes. The expression for the dissociation energy depending on the bilayer curvature can be written as a Gaussian equation

⎡ (c − x )2 ⎤ c ⎥ ΔE = ΔE0 − αin/out exp⎢ − 2w 2 ⎦ ⎣

(9)

where αin/out is the coefficient of the inner/outer layer of a bilayer, and xc and w are constants. As the curvature c tends to zero, the dissociation energy decreases to a minimum ΔE0 which is the dissociation energy of a lipid from a planar bilayer. Using the equation of dissociation energy to fit the PMF results, α for the outer and inner layers are 4.5 and 2.2 nm2 kJ mol−1, respectively; xc is 0.2 nm−1, and w is 0.04 nm−1 (see Figure 7). The curvature c is equal to 1/r, and can be written as

8πN /m , where m is 0.65 nm2. Figure 8. (a) Kinetic model fit to the two-bilayer MLV (8000 lipids) ripening simulation using the present theoretical model. (b) The ripening process of a system with a larger ULV and a small ULV. The two ULVs grow by size disproportionation; i.e., the larger vesicle gets bigger, and the small vesicle gets smaller until the vesicular structure is unable to be maintained, which is consistent with experimental results.1

Kinetics Theory section and the kinetics parameters fitted above, we built a kinetics model for ULV systems with different sizes (the detailed equations are shown in Supporting Information section II). The result demonstrates the size disproportionation mechanism in the aging process of ULVs. The larger ULV gets bigger, and the small ULV gets smaller until it is unable to maintain the vesicular structure (see Figure 8b). Our kinetics model is in agreement with the experimental result.1 The kinetics model for MLV aging can be used to predict the ripening process of larger two-bilayer vesicles (see Figure 9a) or MLVs with more bilayers such as MLVs with four bilayers and eight bilayers (see Figure 9b,c). Usually, such systems are also difficult to study using CG model simulations. Figure 9a shows the larger two-bilayer vesicle of 2 × 104 lipids (which is the same as the MLV at 3.6 μs in Figure 4a) will become a larger ULV as lipids migrate from the inner vesicle to the outmost one. The process is similar to the simulation result seen for the small two-bilayer vesicle in Figure 3a. For the four-bilayer vesicle, the initial size is set according to the simulation data at 0.25 μs in Figure 5b; i.e., the curvatures of the four bilayers are 0.039, 0.052, 0.077, and 0.147 nm−1, respectively. The corresponding numbers of lipids are 2.6 × 104, 1.45 × 104, 6500, and 1800, respectively. The lipid number of the

Figure 7. Dependence of the lipid dissociation energy from the inner and outer layers on bilayer curvature. The fitting yields parameter values in eq 9 as ΔE0 = 22.5 kJ mol−1, xc = 0.2 nm−1, and w = 0.04 nm−1, αout = 4.5 nm2 kJ mol−1, and αin = 2.2 nm2 kJ mol−1.

Figure 7 shows that the dissociation energies of lipids from the inner layer of the first bilayer are always larger than the outer layer of the second bilayer; that is, the lipids prefer to transfer from the second bilayer to the first one. The kinetics equations of MLV ripening can be numerically solved.67 The corresponding parameters kon and Aoff for the coarse-grained lipid MLV can be obtained by fitting the simulation result for the two-bilayer vesicle of 8000 lipids (at stages II and III from 0.27 to 11.9 μs in Figure 3b). As a result, w is 0.13 nm−1, and the fitting values of kon and Aoff are 1.7 × 104 nm3 s−1 and 107.7 s−1, respectively. The fitting result is shown in Figure 8a. On the basis of the kinetics model described in the Ripening F

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

shrink monotonically by lipid diffusion. When the number of lipids of a bilayer is unable to maintain the structure of vesicle, it will merge into its adjacent outside bilayer. After 140 μs, a five-bilayer vesicle is formed. However, as distance increases between the outermost bilayer and the second outermost bilayer, the rate of lipid number decrease of the secondary outer bilayer exceeds that of the third bilayer. Also, then, a MLV with a double-bilayer structure is obtained. Double-bilayer vesicles consisting of single-component 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) were observed in experiment and simulations.68 The ripening process after the double-bilayer structure formation can be further predicted when the doublebilayer is regarded as a pseudobilayer clamping an inversed bilayer in the ripening kinetics (see Figure 9c after 171.8 μs). The inner bilayers of MLV disappear layer by layer, and the MLV finally becomes an ULV.

■

CONCLUSIONS Vesicle ripening is an important problem, although it is difficult to study through experimental and atomistic simulation methods. A mesoscale coarse-grained lipid model is used to simulate the ripening of MLVs. The ripening process for a small MLV is discovered, in which only the outermost bilayer inflates and all the inner bilayers shrink. More complex ripening processes are also found for larger MLVs. The PMFs of a single lipid across curved bilayers are calculated using umbrella sampling to understand the ripening mechanism. The result shows that lipid dissociation from the bilayer has to overcome an energy barrier which depends on the bilayer curvature. While the lipid dissociation energy from the inner layer of a bilayer is higher than that from the outer layer of the bilayer, the aggregation of lipid from solution into the bilayer has no energy barrier. According to the characteristics of the MLV structure and PMFs, a ripening kinetics theoretical model for MLV is built in terms of the kinetic processes. The bilayer curvature dependence of the lipid dissociation energy is considered for the rate constants. The kinetics model is consistent with the simulation and can be applied to predict the ripening process of larger MLVs and other complex vesicle structures. The predicted results show the structure of MLVs are unstable and tend to change into ULVs. The size of the outermost bilayer of a MLV increases monotonically with time, and the others decrease monotonically. The innermost bilayer will merge into the second-innermost bilayer when its size is unable to maintain the structure of the vesicle. When the distance of the divided space between the outermost bilayer and the second bilayer is very large, the lipid number decay rate of the second bilayer is larger than that of the third bilayer, and a double-bilayer structure is observed. Assuming the doublebilayer is a pseudobilayer clamping an inversed bilayer in the kinetics model, the MLV becomes an ULV finally.

Figure 9. Predictions of the ripening processes of (a) a larger twobilayer MLV, (b) a four-bilayer MLV, and (c) an eight-bilayer MLV. Sudden changes of lipid number with time indicate coalescence of the innermost bilayer into the second-innermost bilayer.

outermost bilayer increases, and the lipid numbers of the other bilayers decrease until the innermost bilayer merges into the third bilayer (about 22 μs, see Figure 9b). The coalescence of the innermost bilayer into the third bilayer leads to a sudden change in the number of lipids. Then, the four-bilayer vesicle transforms into a three-bilayer vesicle. Analogously, the threebilayer vesicle becomes a two-bilayer vesicle after 75 μs. Finally, the two-bilayer vesicle transforms to a larger ULV at 140 μs. The more complex ripening of a larger MLV with eight bilayers is also calculated using present theory (see Figure 9c). The initial distance of the adjacent two bilayers is set as 6 nm, and the size of the innermost bilayer of the initial MLV is 6.2 nm in radius. It is shown that the number of lipids in the outermost bilayer still increases monotonically and all the inner vesicles

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b12193. Detailed kinetics equations for the ripening of lipid MLVs and two ULVs with different size (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +86-22-27404303. G

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Notes

(22) Wimley, W. C.; Thompson, T. E. Exchange and Flip-Flop of Dimyristoyl Phosphatidylcholine in Liquid-Crystalline, Gel and TwoComponent, Two-Phase Large Unilamellar Vesicles. Biochemistry 1990, 29, 1296−1303. (23) Kremer, J. M. H.; Kops-Werkhoven, M. M.; Pathmamanoharan, C.; Gijzeman, O. L. J.; Wiersema, P. H. Phase Diagrams and the Kinetics of Phospholipid Exchange for Vesicles of Different Composition and Radius. Biochim. Biophys. Acta, Biomembr. 1977, 471, 177−188. (24) Shioi, A.; Hatton, T. A. Model for Formation and Growth of Vesicles in Mixed Anionic/Cationic (SOS/CTAB) Surfactant Systems. Langmuir 2002, 18, 7341−7348. (25) Guinedi, A. S.; Mortada, N. D.; Mansour, S.; Hathout, R. M. Preparation and Evaluation of Reverse-Phase Evaporation and Multilamellar Niosomes as Ophthalmic Carriers of Acetazolamide. Int. J. Pharm. 2005, 306, 71−82. (26) Maestrelli, F.; Gonzalez-Rodriguez, M. L.; Rabasco, A. M.; Mura, P. Effect of Preparation Technique on the Properties of Liposomes Encapsulating Ketoprofen-Cyclodextrin Complexes Aimed for Transdermal Delivery. Int. J. Pharm. 2006, 312, 53−60. (27) Gauffre, F.; Roux, D. Studying a New Type of Surfactant Aggregate (″Spherulites″) as Chemical Microreactors. A First Example: Copper Ion Entrapping and Particle Synthesis. Langmuir 1999, 15, 3738−3747. (28) Zhang, K.; Gao, L.; Chen, Y.; Yang, Z. Onionlike Spherical Polymer Composites with Controlled Dispersion of Gold Nanoclusters. Chem. Mater. 2008, 20, 23−25. (29) Meyre, M. E.; Lambert, O.; Desbat, B.; Faure, C. Synthesis of Stable, Gold-Particle-Containing Onion-Type Multilamellar Vesicles. Influence of Particle Size on the Onions’ Internal Structure. Nanotechnology 2006, 17, 1193−1201. (30) Meyre, M. E.; Raffard, G.; Franconi, J. M.; Duguet, E.; Lambert, O.; Faure, C. Production of Magnetic Multilamellar Liposomes as Highly T2-Efficient MRI Contrast Agents. Nanomedicine 2011, 7, 18− 21. (31) Cohen, B. E. Concentration- and Time-Dependence of Amphotericin-B Induced Permeability Changes across ErgosterolContaining Liposomes. Biochim. Biophys. Acta, Biomembr. 1986, 857, 117−122. (32) Agnihotri, S. A.; Soppimath, K. S.; Betageri, G. V. Controlled Release Application of Multilamellar Vesicles: A Novel Drug Delivery Approach. Drug Delivery 2010, 17, 92−101. (33) Pott, T.; Roux, D. DNA Intercalation in Neutral Multilamellar Membranes. FEBS Lett. 2002, 511, 150−154. (34) Mignet, N.; Brun, A.; Degert, C.; Delord, B.; Roux, D.; Helene, C.; Laversanne, R.; Francois, J. C. The Spherulites: A Promising Carrier for Oligonucleotide Delivery. Nucleic Acids Res. 2000, 28, 3134−3142. (35) Freund, O.; Mahy, P.; Amedee, J.; Roux, D.; Laversanne, R. Encapsulation of DNA in New Multilamellar Vesicles Prepared by Shearing a Lyotropic Lamellar Phase. J. Microencapsulation 2000, 17, 157−168. (36) Bernheim-Grosswasser, A.; Ugazio, S.; Gauffre, F.; Viratelle, O.; Mahy, P.; Roux, D. Spherulites: A New Vesicular System with Promising Applications. An Example: Enzyme Microencapsulation. J. Chem. Phys. 2000, 112, 3424−3430. (37) Valerio, J.; Lameiro, M. H.; Funari, S. S.; Moreno, M. J.; Melo, E. Temperature Effect on the Bilayer Stacking in Multilamellar Lipid Vesicles. J. Phys. Chem. B 2012, 116, 168−178. (38) Brooks, M. S.; Moggridge, G. D. The Effect of Additives and the Stability of Multilamellar Vesicles in a Commercial Surfactant System. Chem. Eng. Res. Des. 2006, 84, 139−146. (39) Mountcastle, D. B.; Biltonen, R. L.; Halsey, M. J. Effect of Anesthetics and Pressure on the Thermotropic Behavior of Multilamellar Dipalmitoylphosphatidylcholine Liposomes. Proc. Natl. Acad. Sci. U. S. A. 1978, 75, 4906−4910. (40) Reis, O.; Winter, R.; Zerda, T. W. The Effect of High External Pressure on DPPC-Cholesterol Multilamellar Vesicles: A Pressure-

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the National Natural Science Foundation of China (Nos. 20974078, 91127046, and 21274107). The computation is finished in the High Performance Computing Center of Tianjin University and National Supercomputer Center in Tianjin, and the calculations were performed on TianHe-1 (A).

■

REFERENCES

(1) Madani, H.; Kaler, E. W. Aging and Stability of Vesicular Dispersions. Langmuir 1990, 6, 125−132. (2) Olsson, U.; Wennerstrom, H. On the Ripening of Vesicle Dispersions. J. Phys. Chem. B 2002, 106, 5135−5138. (3) Smeijers, A. F.; Markvoort, A. J.; Pieterse, K.; Hilbers, P. A. J. A Detailed Look at Vesicle Fusion. J. Phys. Chem. B 2006, 110, 13212− 13219. (4) Chernomordik, L. V.; Kozlov, M. M. Mechanics of Membrane Fusion. Nat. Struct. Mol. Biol. 2008, 15, 675−683. (5) Kozlovsky, Y.; Kozlov, M. M. Stalk Model of Membrane Fusion: Solution of Energy Crisis. Biophys. J. 2002, 82, 882−895. (6) Fourcade, B.; Mutz, M.; Bensimon, D. Experimental and Theoretical Study of Toroidal Vesicles. Phys. Rev. Lett. 1992, 68, 2551−2554. (7) Ou-Yang, Z. Anchor Ring-Vesicle Membranes. Phys. Rev. A: At., Mol., Opt. Phys. 1990, 41, 4517−4520. (8) Storch, J.; Kleinfeld, A. M. Transfer of Long-Chain Fluorescent Free Fatty Acids between Unilamellar Vesicles. Biochemistry 1986, 25, 1717−1726. (9) Kol, M. A.; de Kroon, A. I. P. M.; Killian, J. A.; de Kruijff, B. Transbilayer Movement of Phospholipids in Biogenic Membranes. Biochemistry 2004, 43, 2673−2681. (10) Wimley, W. C.; Thompson, T. E. Transbilayer and Interbilayer Phospholipid Exchange in Dimyristoylphosphatidylcholine/Dimyristoylphosphatidylethanolamine Large Unilamellar Vesicles. Biochemistry 1991, 30, 1702−1709. (11) Martin, F. J.; MacDonald, R. C. Phospholipid Exchange between Bilayer Membrane Vesicles. Biochemistry 1976, 15, 321−327. (12) Roseman, M. A.; Thompson, T. E. Mechanism of the Spontaneous Transfer of Phospholipids between Bilayers. Biochemistry 1980, 19, 439−444. (13) Jones, J. D.; Thompson, T. E. Mechanism of Spontaneous, Concentration-Dependent Phospholipid Transfer between Bilayers. Biochemistry 1990, 29, 1593−1600. (14) McLean, L. R.; Phillips, M. C. Mechanism of Cholesterol and Phosphatidylcholine Exchange or Transfer between Unilamellar Vesicles. Biochemistry 1981, 20, 2893−2900. (15) Bai, J.; Pagano, R. E. Measurement of Spontaneous Transfer and Transbilayer Movement of Bodipy-Labeled Lipids in Lipid Vesicles. Biochemistry 1997, 36, 8840−8848. (16) Koynova, R.; MacDonald, R. C. Lipid Transfer between Cationic Vesicles and Lipid-DNA Lipoplexes: Effect of Serum. Biochim. Biophys. Acta, Biomembr. 2005, 1714, 63−70. (17) Jähnig, F. Lipid Exchange between Membranes. Biophys. J. 1984, 46, 687−694. (18) Gerelli, Y.; Porcar, L.; Lombardi, L.; Fragneto, G. Lipid Exchange and Flip-Flop in Solid Supported Bilayers. Langmuir 2013, 29, 12762−12769. (19) Nichols, J. W.; Pagano, R. E. Kinetics of Soluble Lipid Monomer Diffusion between Vesicles. Biochemistry 1981, 20, 2783−2789. (20) Needham, D.; Zhelev, D. Lysolipid Exchange with Lipid Vesicle Membranes. Ann. Biomed. Eng. 1995, 23, 287−298. (21) Arvinte, T.; Hildenbrand, K. N-NBD-L-α-Dilauroylphosphatidylethanolamine a New Fluorescent Probe to Study Spontaneous Lipid Transfer. Biochim. Biophys. Acta, Biomembr. 1984, 775, 86−94. H

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Tuning Fourier Transform Infrared Spectroscopy Study. Biochim. Biophys. Acta, Biomembr. 1996, 1279, 5−16. (41) Wang, D.; Dong, R.; Long, P.; Hao, J. Photo-Induced Phase Transition from Multilamellar Vesicles to Wormlike Micelles. Soft Matter 2011, 7, 10713−10719. (42) Pommella, A.; Caserta, S.; Guido, S. Dynamic Flow Behaviour of Surfactant Vesicles under Shear Flow: Role of a Multilamellar Microstructure. Soft Matter 2013, 9, 7545−7552. (43) Gentile, L.; Mortensen, K.; Rossi, C. O.; Olsson, U.; Ranieri, G. A. Multi-Lamellar Vesicle Formation in a Long-chain Nonionic Surfactant: C16E4/D2O System. J. Colloid Interface Sci. 2011, 362, 1−4. (44) Freund, O.; Amedee, J.; Roux, D.; Laversanne, R. In Vitro and in Vivo Stability of New Multilamellar Vesicles. Life Sci. 2000, 67, 411− 419. (45) Koley, P.; Pramanik, A. Multilayer Vesicles, Tubes, Various Porous Structures and Organo Gels through the Solvent-Assisted SelfAssembly of Two Modified Tripeptides and Their Different Applications. Soft Matter 2012, 8, 5364−5374. (46) 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. (47) Xu, R.; Wang, Z.; He, X. Mesoscale Simulation of Vesiculation of Lipid Droplets. Chin. J. Chem. Phys. 2014, 27, 663−671. (48) Cooke, I. R.; Deserno, M. Solvent-Free Model for SelfAssembling Fluid Bilayer Membranes: Stabilization of the Fluid Phase Based on Broad Attractive Tail Potentials. J. Chem. Phys. 2005, 123, 224710. (49) Noguchi, H. Solvent-Free Coarse-Grained Lipid Model for Large-Scale Simulations. J. Chem. Phys. 2011, 134, 055101. (50) Nagle, J. F.; Tristram-Nagle, S. Structure of Lipid Bilayers. Biochim. Biophys. Acta, Rev. Biomembr. 2000, 1469, 159−195. (51) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750−760. (52) Petrache, H. I.; Gouliaev, N.; Tristram-Nagle, S.; Zhang, R.; Suter, R. M.; Nagle, J. F. Interbilayer Interactions from HighResolution X-Ray Scattering. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1998, 57, 7014−7024. (53) Wang, Z. J.; Deserno, M. A. Systematically Coarse-Grained Solvent-Free Model for Quantitative Phospholipid Bilayer Simulations. J. Phys. Chem. B 2010, 114, 11207−11220. (54) Genco, I.; Gliozzi, A.; Relini, A.; Robello, M.; Scalas, E. Electroporation in Symmetric and Asymmetric Membranes. Biochim. Biophys. Acta, Biomembr. 1993, 1149, 10−18. (55) Karatekin, E.; Sandre, O.; Guitouni, H.; Borghi, N.; Puech, P. H.; Brochard-Wyart, F. Cascades of Transient Pores in Giant Vesicles: Line Tension and Transport. Biophys. J. 2003, 84, 1734−1749. (56) Zhelev, D. V.; Needham, D. Tension-Stabilized Pores in Giant Vesicles: Determination of Pore Size and Pore Line Tension. Biochim. Biophys. Acta, Biomembr. 1993, 1147, 89−104. (57) Kolb, D. A.; Weber, G. Quantitative Demonstration of the Reciprocity of Ligand Effects in the Ternary Complex of Chicken Heart Lactate Dehydrogenase with Nicotinamide Adenine Dinucleotide Oxalate. Biochemistry 1975, 14, 4471−4476. (58) Böckmann, R. A.; Hac, A.; Heimburg, T.; Grubmüller, H. Effect of Sodium Chloride on a Lipid Bilayer. Biophys. J. 2003, 85, 1647− 1655. (59) Kol, M. A.; de Kroon, A. I. P. M.; Rijkers, D. T. S.; Killian, J. A.; de Kruijff, B. Membrane-Spanning Peptides Induce Phospholipid Flop: A Model for Phospholipid Translocation across the Inner Membrane of E. coli. Biochemistry 2001, 40, 10500−10506. (60) Liu, J.; Conboy, J. C. 1,2-Diacyl-Phosphatidylcholine Flip-Flop Measured Directly by Sum-Frequency Vibrational Spectroscopy. Biophys. J. 2005, 89, 2522−2532. (61) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (62) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38.

(63) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187−199. (64) Tieleman, D. P.; Marrink, S. J. Lipids out of Equilibrium: Energetics of Desorption and Pore Mediated Flip-Flop. J. Am. Chem. Soc. 2006, 128, 12462−12467. (65) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (66) Higuchi, W. I.; Misra, J. Physical Degradation of Emulsions Via the Molecular Diffusion Route and the Possible Prevention Thereof. J. Pharm. Sci. 1962, 51, 459−466. (67) Mathematica, version 8.0, software for technical computation; Wolfram Research: Champaign, IL, 2010. (68) Sakashita, A.; Imai, M.; Noguchi, H. Morphological Variation of a Lipid Vesicle Confined in a Spherical Vesicle. Phys. Rev. E 2014, 89, 040701.

I

DOI: 10.1021/acs.jpcb.5b12193 J. Phys. Chem. B XXXX, XXX, XXX−XXX