Energetic and Entropic Motifs in Vesicle Morphogenesis in Amphiphilic Diblock Copolymer Solutions

1. Introduction

Self-assembly of block copolymers (BCPs) in solution has generated renewed interest following the pioneering discoveries of vesicular morphologies, referred to as polymersomes [1,2,3,4,5,6,7,8,9,10,11,12]. While traditional ionic surfactants with hydrocarbon chains can form vesicular structures, manipulation of their chemical environment is often required to facilitate vesiculation, e.g., by addition of a salt or co-surfactant, which limits their biocompatibility and biomedical applications. Further, non-ionic copolymer vesicles offer superior stability due to enhanced hydrophobicity as evidenced by their ultra-low critical micelle concentrations compared to those of traditional surfactants. Specifically, the change in chemical potential associated with the exchange of individual copolymer molecules between a molecular assembly and surrounding solvent is large enough to render kinetic time scales spanning hours to days [13]. Further, BCPs offer robust means of size control of the vesicles by varying polymer chain composition/length, addition of secondary spacer molecules, and microfluidics-based flow focusing techniques to produce droplets of controllable sizes: see review by Bleul et al. [1]. Polymersomes with size that ranges between several nanometers to tens of micrometers can be routinely manufactured. They form the building blocks of complex nanostructures that offer significant promise in detergency, environmental remediation, nanomedicine, and cancer theranostics, an approach that integrates diagnosis (e.g., by imaging) and therapeutics [1,14,15,16,17], catalysis and nanoreactor design [18,19,20,21] as well as development of biomaterials that mimic cellular structure and function [22,23,24].

Experimental studies have provided a wealth of information on and insights into the phase behavior of BCPs in solution. Further, they have offered several clues regarding the thermodynamic and kinetic pathways underlying the morphological diversity and structure transitions in such systems. However, a quantitative theoretical understanding of energetic and entropic conduits that facilitate spontaneous changes in aggregate morphology is still lacking due to many reasons. First, real time visualization of the self-assembly process leading to the formation of molecular aggregates in solution is well beyond the spatiotemporal resolution of advanced imaging technologies. Consequently, morphology characterization is achieved often by cryogenic Transmission Electron Microscopy (cryo-TEM) experiments although progress in NMR imaging has made it possible to track the trajectories of individually tagged molecules in native solution as they navigate within and across molecular aggregates [25,26]. Second, experimental observations of shape transitions are often rationalized based on elegant, yet simplified geometric arguments, such as the one based on the packing parameter, which was developed to explain spherical, cylindrical, and vesicular aggregates in amphiphilic surfactants [27]. However, experiments [7,8,9,10,11,12,28,29] and a recently reported systematic coarse-grained molecular dynamics (CGMD) study by the authors [30] have shown that self-assembly of BCPs in solution generates extraordinarily diverse and topologically complex morphologies. Quantitative characterization of structure evolution in such systems requires detailed understanding of the dynamics of molecule exchange between a self-assembled structure and solvent as well as internal reorganization of molecules that is often associated with morphology transitions. In this work, such dynamics are explored in detail in the context of vesiculation in amphiphilic diblock copolymer solutions which involves molecular aggregation in the early stages and internal reorganization of molecules at a constant aggregation number during the latter phase of vesiculation.

Quantitative understanding of entropic (e.g., segmental stretch and orientation of copolymer chains) and energetic (e.g., area of hydrophobic interfaces and potential energy of solvent-monomer interactions) changes underlying morphology transitions, especially rod to lamella as well as lamella to vesicle transitions seen in experiments [28], remains inadequate. At least two pathways have been proposed to explain vesicle formation in BCP solutions. One mechanism entails solvent diffusion into spherical micelles creating a polymer-free core region [1]. The other mechanism, which is widely accepted by scholars due to its perceived universality to phospholipids (biological membranes) and copolymers, is based on energetics arguments, i.e., a sufficiently large lamellar (bilayer) structure develops spontaneous curvature, bends, and folds into a vesicle (or polymersome). Specifically, it is envisioned that during such morphology change, the lamella discards the energy associated with the unfavorable hydrophobic interactions of the polymers along its edge with the surrounding solvent molecules at the expense of gaining curvature energy [1,13,27,29,31,32,33,34,35]. In experiments conducted by Fromherz et al. [23] on egg lecithin – taurochenodesoxycholate systems, electron microscopy revealed that metastable, discoid assemblies of phospholipids spontaneously closed to form vesicular morphologies. This process is thought to be initiated only when the perimeter of the lamella exceeds a threshold value or equivalently, when the driving force provided by the unfavorable hydrophobic edge energy, which is approximately equal to the product of the area of the lateral solvent-lamella interface and interfacial tension, is sufficiently large to overcome the bending energy of the bilayer [31]. Under such circumstances, geometric arguments, e.g., conservation of surface area of the copolymer aggregate, can provide a relationship between the linear dimensions of the lamella and vesicle produced: see for example Figure 5 in Antonietti and Förster [29].

The abovementioned conceptualization of vesicle generation from bilayers is qualitatively supported by computer simulations that utilize solvent-free multi-particle collision dynamics and discrete representations of fluid membranes parameterized to manipulate particle diffusion and membrane rigidity [36,37]. Such discrete particle simulations can capture dynamic organization of fluid membranes and provide significantly more nuanced depictions of the pathways of morphology transitions compared to those offered by classical continuum mechanics theories based on energy minimization principles that govern bending/buckling of planar sheets to form curved and closed surfaces [38,39,40,41,42]. However, they lack critical details of the organization of constituent molecules within the aggregates, which is largely dictated by solvent-mediated interactions. In the context of aqueous solutions of copolymers, distributions of segmental stretch and water within the assemblies are highly heterogeneous and greatly dependent on the chain composition, confinement (e.g., water-hydrophilic segment interface in the inner and outer layers of a unilamellar vesicle), hydrogen bonding interactions, and/or solvent mobility [30,43]. Moreover, aggregate morphologies are dynamic with irregular shapes and diffuse interfaces [20,21,30]. Such molecular-scale aspects of shape transitions are addressed in the present study.

Experiments by Chen et al. [28] show the coexistence of rodlike micelles, lamellae, and vesicles in polystyrene-polyacrylic acid solutions in dioxane/water. These authors inferred that rod to vesicle transition involves an intermediate lamellar state by visualizing a series of structures by cryo-TEM. The transition from a solution that contains predominantly rods to one that is dominated by vesicles was facilitated by perturbing a solution of rodlike micelles at equilibrium, by increasing its water content by incremental amounts. The coexistence of rods, lamellae, and vesicles were quantified using turbidity data. Specifically, each morphology makes a distinctive contribution to the overall turbidity of the solution which was quantified using experimental measurements. Further, a kinetic model based on a two-state transition theory in which (a). rods are flattened to form a lamella, and (b). develops curvature and closes on itself to create a vesicle was proposed. Intermediate morphologies between spherical micelle and vesicles, such as cylindrical micelle, lamellar (discoid) and open vesicle (curved lamellar) have been observed in experiments for amphiphilic systems [44,45]. In cetyl-trimethyl-ammonium-bromide and 5-methyl salicylic acid solution, thermally reversible transition from vesicles to wormlike micelles has been observed [46]. In phospholipid solutions containing vesicles, spherical micelle formation can be triggered by the addition of surfactants [44,45]. The underlying pathway consists of opening of the vesicle to form a curved bilayer followed by bilayer to rodlike micelle transition and formation of spherical micelles from rods. In our earlier molecular dynamics study [30], co-existence of rods, bilayers and vesicles was predicted in diblock polymer solutions. In this work, we investigate the sequence of morphology transitions leading to vesicle formation in diblock copolymer solutions using CGMD simulations.

CGMD simulations conducted by Markvoort and coworkers [47,48] on phospholipid systems predict an increase in the overall potential energy associated with bilayer to vesicle transition, suggesting the existence entropic driving forces. Whether this is true in the case of copolymer solution is unknown. Configurational entropy associated with the distributions of chain stretch and orientation of hydrophobic/hydrophilic segments within BCP aggregates could influence their macroscopic mechanical properties. This in turn could have pronounced implications on shape transitions. Establishing quantitative relationships between molecular structure and mechanical properties of BCP aggregates is beyond the scope of the present study. However, we compute probability distribution functions of segmental stretch for rodlike, lamellar, and vesicular aggregates to glean insights into potential entropic driving forces underlying vesiculation.

The present study is built on our recently reported CGMD simulations of spontaneous self-assembly in model (PB-PEO) diblock copolymer solution leading to the formation of diverse morphologies [30]. The initial conditions of our simulations correspond to a homogeneous BCP solution at a prescribed concentration in which polymer molecules of given chain length and composition are randomly distributed. No pre-assembly is employed, and spatiotemporal evolution of structures is tracked in real time. In the past, mesoscopic simulations based on self-consistent field theory [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51], dissipative particle dynamics [52,53,54,55,56,57,58,59,60], and CGMD simulations with pre-assembled initial conditions [20,21,61] have been employed to study self-assembly of amphiphilic BCPs in solution, bilayer to vesicle transition, and kinetics of copolymer exchange between self-assembled aggregates in BCP solutions. Further, coarse-grained molecular dynamics simulations that account for solvent-mediated hydrophobic interactions implicitly via a suitably parameterized Lennard-Jones (LJ) potential function have been used to study differences in copolymer architecture (linear vs. bottlebrush) on self-assembly in solution [62] and shape deformation of PEO-PS polymersomes induced by osmotic pressure stress [63]. CGMD simulations have also been used to qualitatively study the self-assembly process of PS-PAA polymers [64], phospholipid [65], and DTAB and gemini surfactants [66]. In the present study, we employ the CGMD simulation framework reported in [30] which were adapted from our previous investigations of self-assembly and shape transitions in as well as shear/elongational flow dynamics and rheology of cationic surfactant solutions [67,68,69,70,71,72,73,74]. We use the MARTINI force field to describe all bonded and non-bonded interactions within the BCP solution [75]. Our simulations capture vesiculation from an initially homogeneous copolymer solution. This process, which is facilitated through multiple intermediate stages of micellization and structure reorganization, is analyzed to provide quantitative insights into the underlying energetics and entropic motifs. Simulation details and data analytics methods are discussed in Section 2. Results are discussed in Section 3. Conclusions and outlook are offered in Section 4.

2. Materials and Methods

2.1. Methods and Software

Lenard-Jones and harmonic potential functions are used to represent non-bonded and bonded pairwise interactions, respectively [54]. The parameter values used in this study are identical to those reported in Liu and Sureshkumar [30] in which morphological diversity in diblock BCP solutions mimicking PEO-PB architectures of varying chain length, chain composition and concentration was explored.

In the present simulations, 400,000 water beads, 40,000 antifreeze water beads and 9,000 polymer chains consisting of six PEO and six PB beads (14.4 wt%) are randomly placed in an initial cubical box of linear dimension 40 nm by using PACKMOL software. The simulations are performed by using the Gromacs 2020.2 MD software. The production runs are NPT simulations. They span 400 ns which is approximately 10% longer than the time required for a stable vesicular structure to evolve from an initially random configuration of copolymer chains. The molecular visualization is performed by VMD 1.9.3. Data analysis and graphing are performed using Matlab^TM 2023, Microsoft 365^TM, C++ and Originlab^TM 2017.

2.2. Data Analytics

Discrete spatiotemporal data generated from the CGMD simulations are analyzed to identify copolymer clusters. Subsequently, their geometric features such monomer number N, volume V, surface area A, and monomer density ρ ≡ N/V are computed. We also calculate solvent accessible surface area (SASA), which may be interpreted as the surface area of the water−aggregate interface. Specifically, SASA is measured as the area circumscribed by the motion of the center of a spherical probe that is slid along the periphery of the beads that constitute a cluster [76]. The algorithms utilized for the above calculations are described in detail in our earlier work [30].

End-to-end Distance and Body Gaussian Curvature

To quantitatively describe shape changes of the self-assembled structures, we calculate the end-to-end distance Q and body Gaussian curvature K_G of the entire assembly with respect to its center of mass. Note that K_G is a wholistic (body) curvature measure of the geometry of the aggregate rather than the integral of the pointwise Gaussian curvatures along its surface. This concept is illustrated in Figure 1 for two representative surfaces, one that is relatively flat along (Figure 1a–c) and the other curved (Figure 1d–f). The Cartesian (x, y, z) coordinate system used has the center of mass of the aggregate as its origin. The shortest linear dimension is along the y axis. Figure 1b shows the collections of discrete points (“point cloud”) that lie within -2.5 nm ≤ z ≤ 2.5 nm of the structure whose mesh representation is given in Figure 1a. This set of points are projected onto the x-y plane and binned into several intervals: see Figure 1b,e. We use 20 and 60 intervals for flat and curved surfaces, respectively as they provide sufficient discretization to calculate numerically converged results for the curvature. A structure is considered flat if its center of mass lies within its point cloud and curved if otherwise. The curvature in the x-y plane is calculated as the inverse of the radius of the circle (Figure 1c) that fits the centers of mass of the beads in each bin (red dots in Figure 1b), by using a robust fitting algorithm developed by Pratt [77]. We discard data for the first two (four) and last two (four) bins for flat (curved) structures, to minimize the uncertainty in the calculations arising from large fluctuations of the edges. The curvature is considered positive (negative) if the center of the best fit circle lies above (below) the center of mass of the structure. The end-to-end distance of the point cloud Q_p is calculated as the magnitude of the vector connecting the extreme points in the two outermost bins. The above procedure is repeated by rotating the cluster in increments of 10 degrees around the y-axis to a total of 180 degrees. The body Gaussian curvature K_G ≡ K₁K₂ where K₁ and K₂ denote the maximum and minimum curvatures. The end-to-end distance Q of the structure is defined as the maximum of the computed Q_p values.

Information Entropy of Aggregates

To quantify the changes in the configurational entropy associated with the packing of copolymers within various intermediate assemblies leading to vesiculation, we calculate the probability distribution function (pdf) of the extension of the hydrophobic (PB) and hydrophilic (PEO) segments of the polymer chains. The associated structure parameter s ≡ q_max − q, where q_max is the maximum end-to-end distance between the first and last beads of the PB or PEO segment within a structure. The pdf of s is normalized such that

$$\frac{1}{q_{max}-q_{min}}{\int }_0^{s_m}p\left(s\right)ds=1,$$

where s_m = q_max − q_min. Given p(s), we calculate the information entropy H of the structure using its standard definition as H(s) ≡ ∑ p(s) ln (p(s)) where the discrete summation is carried out for the ensemble of copolymers within a cluster.

3. Results and Discussion

We study the evolution of nearly spherical polymersomes made up of symmetric BCP chains which consist of 6 PB and 6 PEO beads corresponding to a PEO mole fraction of 50%. The BCP concentration is 14.4 wt%. The chain lengths used in the simulations translate to a polymer molecular weight of 590 Da. Hence, these are relatively short polymers with molecular weights appreciably greater than those of surfactants with long hydrocarbon chains, e.g., cetyl-trimethyl ammonium chloride, a C-16 cationic surfactant studied extensively in experiments and our previous CGMD simulations, has a molecular weight of 320 Da.

To quantify the geometry of the structure, we tracked the number of monomers (N), body Gaussian curvature (K_G) and end-to-end distance of the self-assembled structure (Q) as functions of time. Representative results are shown in Figure 2 and Figure 3. Vesiculation is facilitated through the following pathway: formation of well-defined spherical aggregates from a homogeneous copolymer solution, merger of spherical aggregates to form rodlike micelles (t < 100 ns, Figure 2A–F), fusion of rods to form longer wormlike micelles (100 ns < t < 120 ns, Figure 2F–H), flattening of flexible cylindrical structures to form rectangular lamellae (bilayers) which further reorganize into disk-shaped lamellae (120 ns < t < 250 ns, Figure 2H–N), bending and curving of disklike lamellae leading to cavity formation, and closure of the cavity to form vesicles (250 ns < t < 375 ns, Figure 2N–Q). A video (Video S1) of this process is provided in SI.

The early stage of vesiculation is characterized by copolymer aggregation into spherical micelles. Each one of the three step changes seen in N in Figure 3 corresponds to a merger event. First, small spherical aggregates consisting of ~1000 monomers coalesce to form a rodlike micelle. Second, two rodlike micelles approach each other at an obtuse angle (see Figure 2F) and combine into a longer wormlike micelle (WM) by an endcap opening mechanism. Specifically, it is energetically favorable for cylindrical micelles to merge as this process results in the loss of the curvature energy associated with their endcaps [72]. Third, the wormlike micelle accommodates a small spherical aggregate and grows larger. The first two mergers occur at terminal positions and are accompanied by an increase in Q while the third merger occurs on the side of the rodlike micelle resulting in little change in Q. For the rodlike structures shown in Figure 2E–H, positive maximum body curvature appears along the major axis, and negative minimum body curvature appears along the minor axis. Hence, K_G is negative.

After the third merging process, the maximum and standard deviation in N is 0.28% and 0.06%, respectively. Consequently, vesicle formation from the rodlike structure occurs through a series of structure re-organization at a practically constant aggregation number n = N/12 as seen from Figure 3. For 140 ns ≤ t ≤ 180 ns, the reordering of the copolymers within the flexible WM results in the formation of a rectangular lamella (RL). WM to RL transition is accompanied by a continuous decrease in Q. Early in the process, due to the initial angle of approach between the two rodlike micelles during the merging process shown in Figure 2G, a topological twist is formed within the RL, thus creating a saddle point: see row 3 in Table 1. For 180 ns < t ≤ 215 ns, the cluster rearranges to a flat circular disk lamella (DL) to increase symmetry. During this transition, K_G approaches 0 and Q decreases further. For 215 ns < t ≤ 375 ns spontaneous curvature fluctuations of the DL initiate a curling process. Subsequently. the structure begins to bend and form a spherical cavity (C). This is accompanied by a reduction in Q and an increase in K_G. For t > 375 ns, a stable spherical bilayer structure (vesicle, V) is established by the closure of the cavity. Simulations predict that the bending process is characterized by a larger density of monomers (ethylene oxide beads) in the internal leaflet: see Figure 4.

Table 1 summarizes the five prominent structures observed along the vesiculation pathway, namely wormlike micelle (WM), rectangular lamella (RL), disklike lamella (DL), cavity (C) and vesicle (V). Images consisting of blue and red domains depict the entire copolymer, while those with red surfaces reveal only the PB segments. The table also shows the particle pair correlation functions, g(r), and the probability distribution functions, p(s), for the structure parameters corresponding to the segmental stretches of PB (q_PB) and PEO (q_PEO). Variations in g(r) for PB (hydrophobic), PEO (hydrophilic), and water are shown in red, blue, and black lines, respectively. These functions quantify molecular density profiles within the structures and provide information on the relative sizes of the hydrophilic and hydrophobic domains. A WM is devoid of water inside and has a PB core surrounded by PEO shell exposed to water. RL, DL, and C all exhibit a PEO-PB-PEO bilayer structure. Vesicles possess a water core and a PEO-PB-PEO bilayer shell. It is intriguing that p(s) variations for all structures are skewed toward lower s (larger segmental stretch q) values and better fitted to a log-normal rather than Gaussian (normal) distributions, discounting the second minor peak seen for PEO. This second peak in the tail of the pdfs corresponds to ring-like configurations primarily located at the overlapping PEO-PB interfaces: see paragraph below. Diffuse (overlapping) interfaces predicted by our simulations are consistent with those observed in previous CGMD simulations of flexible diblock BCPs [4,63]. The average segmental stretch is expectedly greater for the hydrophilic component (PEO). The p(s) of PB features relatively long and low-probability tails, which is a consequence of large chain length fluctuations in the less extended (coiled) PB configurations.

To understand the log normal variations in s, we conducted simulations to track the equilibrium configuration dynamics of a single chain in solution. Analysis of the data collected for the instantaneous values of q_PEO and q_PB showed that the chain stretch parameter s follows log-normal distribution for PEO and PB segments as well as the entire diblock chain, as shown in Figure 5. Consequently, the pdfs of s for assemblies of individual polymer chains are also log-normal. The qualitative agreement of p(s) to a log normal distribution might have implications in terms of entropy maximization. Specifically, a log normal distribution maximizes the information entropy of a random variable with prescribed values of both mean and variance [78]. The average segmental stretch and configurational fluctuations are bounded by the ring-like and fully stretched states of the segments. The skewness and kurtosis in the distributions appear to depend on the aggregate shape. This suggests that the details of molecular packing within the aggregates are topology dependent. Further investigations are needed to quantify such structure-topology relationships.

To investigate whether vesicle formation is entropy-driven, we calculated the information entropy, H along the vesiculation pathway within the constant N regime. Figure 6a shows that H values for the PB and PEO segments remain practically unchanged during the vesiculation process, with H_PEO and H_PB = 4.62 ± 0.10 nat and 5.02 ± 0.02 nat, respectively. Figure 6a also shows that H_PB exhibits relatively larger fluctuations. Upon examining q_max and q_min (Figure 6b), we find that for PEO, q_max and q_min remain practically constant along the vesiculation pathway. However, for PB, q_max shows negligible changes while q_min undergoes relatively large fluctuations. This observation can be attributed to the fact that q_max is related to the fully stretched segment length, which does not change substantially. On the other hand, q_min is associated with the shortest coiled configurations. Based on their chemical structures, the parameter K_θ in the bending potential function for PB is smaller than that for PEO (K_θ = 25 kJ/mol and 50 kJ/mol for PB and PEO, respectively) [75,79]. A lower K_θ value results in greater fluctuations of the coiled configurations. In the absence of systematic statistically significant variations in the information entropy, we conclude that entropic driving forces are not responsible for vesiculation in the system studied.

We now examine the energetics-based mechanism proposed for lamella to vesicle transition [27,34]. Towards this end, we compare the relative change in the surface area △S ≡ 1 - S_v/S_d associated with the transition where S_v and S_d denote the surface area exposed to the solvent for the vesicle and disklike lamella, respectively. For an ideal lamella of radius R_d and thickness L, S_d = 2πR_d²(1 + L/R_d). If such a lamella bends, curves, and folds to form an ideal vesicle of radius R_v, the resulting surface area S_v = 4π [R_v² + (R_v – L)²] where R_v ≈ R_d/2. However, this assumes perfect geometric shapes and that packing density is conserved during the structure transition process. To assess the impact of such non-idealities on △S, we also compute S_d and S_v directly from the simulations. The procedure is illustrated in Figure 7 for the ideal (Figure 7a) and actual (Figure 7b) cases for R_d = 9.9 nm and L = 4.7 nm. Figure 7b,c show that the bilayer thickness L remains unchanged during the lamella to vesicle transition, as assumed in classical theoretical models. However, the radius of the resulting vesicle is significantly larger than that predicted by ideal geometric conditions, i.e., simulations yield R_v/R_d = 0.73. This implies that the size of the water core, and consequently the surface area, of the vesicles formed are larger than those predicted by models that use ideal geometric representations of the self-assembled structures. This translates to a smaller value of △S ≈ 0.18 predicted by the simulations vs. △S ≈ 0.66 calculated for idealized structures.

Figure 8 shows the change in solvent accessible surface area (SASA) along the vesiculation pathway. Note that SASA is different from the macroscopic surface area calculated from the linear geometric dimensions shown in Figure 7: specifically, as explained in [30], SASA computations incorporate interfacial area variations on the scale of the monomer bead size. Figure 8a shows that the PEO-water interface accounts for nearly 95% of the cumulative SASA. Moreover, SASA decreases continuously during rod to RL and DL to vesicle transitions. However, as the RL reorganizes into a DL, SASA is unchanged. This observation is consistent with the changes in the packing density during vesiculation shown in Figure 9. Specifically, during the transition from RL to DL the packing density and SASA are unaffected. Furthermore, although the SASA associated with PB segments is much smaller compared to that of PEO, the relative change in SASA during rod to vesicle transition of the PB segments (≈ 0.3) is significantly greater than that of PEO segments (≈ 0.1). This observation supports the edge energy hypothesis [1,13,27,29,31,32,33,34,35] that the unfavorable edge energy of the hydrophobic boundaries is the driving force for vesicle formation. The quantitative differences in SASA between theoretical and simulation predictions can be attributed to the shape imperfections of the aggregates and changes in the packing density. Figure 9 clearly shows that the polymers within the vesicle are more compactly packed as compared to those in the lamella. The reduction in SASA manifests as lower energies of polymer-solvent interaction as shown in Figure 10. Specifically, the (non-bonded) interaction energy of either PB or PEO segments with water is seen to become less negative during vesiculation.

Finally, the robustness of the above conclusions was established by performing simulations with different initial configurations of the system: see Figures S1 and S2. Further, the isentropic pathway of vesiculation was seen in simulations in which bond stretching and bending force parameters were varied with 10% of their values.

4. Conclusions

CGMD simulations of real-time evolution of vesicles, starting with an initial condition of randomly distributed diblock copolymer chains in an explicit solvent, have been performed to identify intermediate self-assembled states along the vesiculation pathway. The vesiculation process is facilitated by rapid formation of spherical aggregates, the merger of spherical aggregates to form rodlike assemblies, fusion of rods to form wormlike micelles, flattening of a WM to create a rectangular lamella, reorganization of the RL into a disk-shaped lamella, bending and curving of the DL leading to cavity formation, and closure of the cavity to form a vesicle. Analysis of information entropy shows that vesiculation from a wormlike micelle is an isentropic process. SASA decreases during vesiculation, with a more significant relative surface area reduction occurring in the PB-water interface. These predictions support the edge energy hypothesis that vesiculation is accompanied by a loss in the energy of hydrophobic interactions between solvent and copolymers along the lateral interface of the lamella, at the expense of gaining curvature energy. However, simulations show that the rod-lamella-vesicle transition is accompanied by an increase in copolymer packing density. Hence, the change in the surface area accompanying vesiculation predicted by the simulations is significantly lower than theoretical estimates based on idealized geometric representations of lamellar and vesicular morphologies.