Vesicle Morphogenesis in Amphiphilic Triblock Copolymer Solutions

: Coarse-grained molecular dynamics simulations are employed to investigate the spatiotemporal evolution of vesicles (polymersomes) through the self-assembly of randomly distributed amphiphilic BAB triblock copolymers with hydrophilic A and hydrophobic B blocks in an aqueous solution. The vesiculation pathway consists of several intermediate structures, such as an interconnected network of copolymer aggregates, a cage of cylindrical micelles, and a lamellar cage. The cage-to-vesicle transition occurs at a constant aggregation number and practically eliminates the hydrophobic interfacial area between the B block and solvent. Molecular reorganization underlying the sequence of morphology transitions from a cage-like aggregate to a vesicle is nearly isentropic. The end-to-end distances of isolated copolymer chains in solution and those within a vesicular assembly follow lognormal probability distributions. This can be attributed to the preponderance of folded chain configurations in which the two hydrophobic end groups of a given chain stay close to each other. However, the probability distribution of end-to-end distances is broader for chains within the vesicle as compared with that of a single chain. This is due to the swelling of the folded configurations within the hydrophobic bilayer. Increasing the hydrophobicity of the B block reduces the vesiculation time without qualitatively altering the self-assembly pathway.


Introduction
Block copolymers (BCPs) have been a focal point of research for several decades.They have been engineered to produce various morphologies by altering the combinations of chemically distinct polymer segments that interact selectively with their environment, such as an aqueous/organic solvent or a polymer matrix.The ability to tailor morphological features by manipulating polymer chain length, composition, and monomer chemistry has led to their applications in areas such as advanced material manufacturing, catalysis, emulsification, environmental remediation, targeted drug delivery, gene therapy, and medical diagnostics [1][2][3][4][5][6][7][8][9][10][11][12].
The self-assembly of BAB triblock copolymer solutions, in which B and A represent the solvophobic and solvophilic blocks, respectively, was investigated by Kotaka et al. [43], who reported that poly(methylmethacrylate)-polystyrene-poly(methylmethacrylate) BCPs in a toluene-p-cymene (pCY) solvent mixture with high pCY content form highly branched structures.Balsara et al. [44] suggested that micellization in BAB triblock solutions with solvophilic A blocks could lead to morphologies with B-rich cores surrounded by loops of A segments.A series of quasi-elastic light-scattering experiments conducted on poly(vinylpyridine)-polystyrene-poly(vinylpyridine) solutions in toluene, which preferentially solvates the polystyrene block, confirmed only the presence of individual micelles [44].However, certain BCP-solvent combinations promote the formation of networks in which adjacent solvophobic spherical aggregates are connected by extended A strands.Such branched structures are arguably thermodynamically favorable since they avoid the entropic penalty associated with loop formation.In subsequent experiments with different triblock polymer-solvent systems, network structures facilitated by extended solvophilic segments have been directly observed or inferred [45][46][47][48][49][50][51][52].Overall, the evolution and stability of such structures depend on the entropic and energetic contributions to the system's free energy, which are, in turn, dependent on the molecular architecture and polymer-solvent interactions.
Depending on the chain length, the mole fraction of A/B and monomer chemistry, BAB, and ABA copolymers can form different nanostructures, such as lamellae, star-like or flower-like micelles, and vesicles [53][54][55].In general, BAB architectures with solvents with a relative preference for the A block readily form bilayers and vesicles (polymersomes) [56][57][58][59][60][61].Hence, they are excellent candidates for the synthesis of vesicles that are extensively used to encapsulate therapeutic agents for targeted drug delivery [62,63] or stimuli-responsive nanoparticles for applications in medical imaging [64].In comparison, ABA polymers tend to readily form network structures and are often used to produce thermoreversible gels [65][66][67].
Self-consistent field theory [30][31][32], dissipative particle dynamics (DPD) [33][34][35][36][37][38][39][40][41], and coarse-grained molecular dynamics (CGMD) simulations with prescribed or biased initial conditions [28,29,42] have been employed to study bilayer-to-vesicle transitions and the kinetics of copolymer exchange between self-assembled aggregates in BCP solutions.Further, CGMD simulations that account for solvent-mediated hydrophobic interactions have been implicitly used to study differences in copolymer architecture (linear vs. bottlebrush) with respect to self-assembly in solution [68] and the shape deformation of polymersomes induced by osmotic pressure stress [69].Mattice et al. used Monte Carlo simulations to study self-assembly in ABA copolymers [70][71][72] and found that they form a rich variety of structures.In their simulations, the formation of branched structures in ABA triblock copolymer solutions was observed under conditions in which the corresponding AB diblock copolymers would provide steric stabilization for a polymer colloid.Chen et al. [73] explored the self-assembly pathways in coil-rod-coil triblock copolymers in a rod-selective solvent using DPD simulations and constructed a phase diagram of the predicted morphologies.They observed the formation of different types of micelles while changing the relative block fractions and rod length.These authors also identified a kinetic process leading to vesiculation that involved intermediate structures such as rodlike micelles and lamellae, as observed in CGMD simulations of AB BCPs [27].Han et al. [74] reported two mechanisms of vesicle formation through ABA polymers with Monte Carlo simulations, which are facilitated by either solvent diffusion into a spherical micelle or the deformation of a lamellar structure.Kangarlou et al. [75] performed CGMD simulations of polyethylene oxide molecules end-capped with cholesterol in water by using the MARTINI force field and predicted the formation of flower-like micelles.Song et al. [76] conducted a Monte Carlo simulation study on the effect of the hydrophilicity of the B block on the self-assembly of cyclic AB and BAB polymers.They found that increasing the compatibility between the B block and solvent results in a sequence of morphological changes in the following order: network structures, cylindrical micelles, disklike (lamellar) and toroidal-like micelles, vesicles, and multicompartment vesicles.
Despite the insights gained from the abovementioned studies, simulations with nearmolecular-scale resolutions that capture the real-time evolution of the self-assembly process in the presence of explicit solvent-mediated interactions in triblock copolymer solutions are lacking.As isolated copolymers dispersed within a solvent are brought together by thermodynamic forces to form a molecular assembly, their individual configurations undergo significant changes due to inter-chain interactions.Our previous work on vesicle formation in AB BCPs [27] revealed that, once a molecular aggregate of a critical size is formed, the isentropic reorganization of polymer chains within a constrained geometrical environment directs the pathway of morphology evolution, guided by the motif of reducing hydrophobic contacts.CGMD simulations of vesiculation in BAB BCPs reported in this work track configurational changes at the single-chain level and, hence, provide a fuller understanding of the energetic and entropic mechanisms underlying the self-assembly process.
The equilibrium configurations of BAB and AB copolymers in solvents with greater affinity for the A block are qualitatively different.A comparison of the probability distribution functions of the magnitude of the end-to-end vector of AB and BAB polymers of the same chain length and the A/B molar ratio is provided in Figure S2 in the SI, along with depictions of representative chain topologies.Specifically, the most probable configuration for AB copolymers corresponds to a semi-stretched state.In comparison, BAB copolymers exhibit more complex configurational dynamics resulting in dominant equilibrium topologies that resemble hairpins, rings, and "ρ" shapes, in which a hydrophobic end group stays close to another hydrophobic portion of the chain.Such differences in chain architecture/topology could have a pronounced influence on vesiculation pathways.This is explored in this work by performing CGMD simulations of vesicle formation in BAB copolymer solutions.
The CGMD simulations reported in this work are adapted from previous studies on self-assembly, shape transitions, and the rheology of surfactant solutions [26,27,[77][78][79][80][81][82][83][84].They utilize the MARTINI coarse-graining approach and force fields [85].These simulations provide a detailed picture of vesicle formation from an initially homogeneous copolymer solution, capturing various intermediate (transient) morphologies.We present the simulation methodologies and data analysis techniques in Section 2, the results in Section 3, and the conclusions in Section 4.

Methods and Software
The force field constraints of the molecular models used in this work were based on the requirements of the MARTINI force field [85].
Stretching and bending interactions between bonded beads were modeled using weak harmonic potentials V b and V θ , respectively, provided by Equations ( 1) and (2) below: where b 0 and θ 0 represent the equilibrium bond distance and equilibrium bond angle at the minimum energy configuration, respectively, and K b and K θ represent constants associated with bond-stretching and -bending energies, respectively.The non-bonded interactions between beads i and j are described by a 6-12 Lennard−Jones (LJ) potential represented by Equation ( 3) can also be written as where In the above equations, r represents the distance between beads i and j, ε ij is the depth (minimum) of the potential well (energy function), and σ ij is the distance at which V LJ is zero.ε ij and σ ij are also referred to as the interaction strength and cutoff length, respectively.
The MARTINI force field requires the LJ potential to be shifted smoothly to 0 between 0.9 and 1.2 nm rather than using a hard cutoff to avoid singularities in the computation of the inter-particle forces.This is achieved by using a potential−switch function, S V , represented by ij S (12) − C (6) ij S (6) , (6) where α = 6 or 12.
In Equation (7), parameters A, B, and C are defined by and where the values r shi f t = 0.9 nm and r cut = 1.2 nm are used.The force field parameter values used in this study are listed in the Supplementary Information as Tables S1-S3.
In the present simulations, water (four-molecule cluster, represented by W), PB monomer (represented by B), and PEO monomer (represented by A) are modeled by bead types P4, C4, and SNda, respectively [85,86].In total, 400,000 water beads; 40,000 antifreeze water beads (modeled by bead type WF) [87]; and 3883 polymer chains consisting of five PB, ten PEO, and five PB beads (10.7 wt%) were randomly placed in an initial cubical box with a linear dimension of 40 nm using the PACKMOL 20.14.0 software [88].The spring constant, K b , for the harmonic bond-stretching potential of the PEO (A) segment is an order of magnitude greater than that for the PB (B) block.Further, K θ for the A block is twice that of the B block.Resultantly, the middle (solvophilic A) block is stiffer than the side (solvophobic B) segments.Hence, the BAB system considered in our work has certain similarities to the mechanical properties of a coil-rod-coil copolymer, which has been shown by DPD simulations to form vesicular structures [73].The simulations were performed using the Gromacs 2020.2MD software.The reference pressure and temperature were 1 bar and 300 K, respectively.First, simulations were conducted to minimize the energy of the initial system by using the steepest descent algorithm.This was followed by a short (2 ns in this work) NVT simulation to equilibrate the system at a desired temperature to ensure algorithmic stability.Subsequently, a sufficiently long NPT simulation (production run) was carried out (600 ns in this work).In the subsequent discussions, the state of the system at the start of the NPT simulations was referred to as the initial condition (t = 0).More details about the simulation are provided in the Supplementary Information as Section S1.The molecular visualization was performed using VMD 1.9.3.Data analysis and graphing were performed using Matlab TM 2023, Microsoft 365 TM , Visual C++ 2022, and Originlab TM 2022.

Data Analytics
Discrete spatiotemporal data generated from the CGMD simulations were analyzed to identify copolymer clusters.Any two beads with center-to-center separations less than the cutoff distance were regarded as part of a cluster (aggregate).Further, within any bead in each cluster, there exists at least one bead at a distance less than the cutoff value.In this study, we used 0.5 nm as the cutoff distance.This selection is based on the cutoff distance of 0.47 nm for non-bonded LJ interactions between PB and PEO monomers.We also calculated solvent-accessible surface area (SASA), which may be interpreted as the surface area of the solvent-polymer interface.Specifically, SASA was measured as the area circumscribed by the motion of the center of a spherical probe that slides along the periphery of the polymer beads that constitute a cluster [89].

Information Entropy (H)
To quantify the changes in the configurational entropy associated with the packing of copolymers within various intermediate assemblies leading to vesiculation, we calculated the probability distribution function (pdf ) of the extension of the polymer chains.First, the end-to-end distance of each chain (Q) in the cluster was calculated.Given p(Q), we calculated the information entropy, H, of the structure using its standard definition as [90], where the discrete summation was carried out for the ensemble of copolymers within a cluster.The end-to-end vector is the most widely used structure parameter in polymer kinetic theory [91].Specifically, p(Q) quantifies the distribution of polymer stretch, which depends on the polymer configuration.A broader Q distribution points to greater configurational diversity and, hence, larger degrees of configurational freedom (entropy).This is discussed further in Section 3.1.1.

Pair Correlation Function (g(r))
Pair correlation function (g(r)) describes how the density of B, A, or W beads varies as a function of distance from the center of mass of a cluster along the radial direction.A larger value of g signifies more tightly packed beads.In general, g shows (local density)/(overall density).g(r) is represented by where < ρ(r) > is the average density of beads at a distance, r, from the center of mass; ρ = 3N m / 4πr 3 max is the average cluster density; N m is the total number of monomer beads contained in the cluster; r max is the distance of the farthest bead from the center of mass of the cluster; V(r) = (4/3)π (r + ∆r) 3 − r 3 ; and Equation ( 12) means beads at the inner boundary are included.while those at the outer boundary are excluded.In the calculation, we used ∆r = 0.1 nm.

Morphology Evolution
Figure 1 illustrates the evolution of morphology during vesiculation in BAB copolymers in the NPT simulation.The structure evolution during the NVT simulation is provided in Figure S1 as the Supplementary Information.The polymers form core (B)-corona (A) micelles that are interconnected by bridges of elongated A segments: see Figure 1a.This structure is unstable due to its relatively large solvent-polymer interfacial area.Hence, it coarsens by the merger between neighboring micelles and, to a lesser extent, the incorporation of small polymer aggregates from the solution (Figure 1b).The coarsening process continues until a cage-like morphology with a circumferential mesh of cylindrical micelles is established (Figure 1c,d).The cylindrical micelles rearrange and coalesce with adjacent structures to form a large lamellar (bilayer) shell with holes (Figure 1e).This structure subsequently evolves into a fully closed bilayer (vesicle) (Figure 1f,g).As seen in Figure 2, the closure of the bilayer is accompanied by a reduction in the volume (or, equivalently, the radius of gyration, R g ) of the structure.In the early stages of structure evolution, the aggregation number, N, grows rapidly due to copolymer aggregation.The step changes in N seen in the later stages of morphology evolution are caused by the addition of small polymer aggregates into the primary structure.Overall, once the cage-like structure shown in Figure 1d is formed, the vesiculation process is facilitated largely through its reorganization at a nearly constant aggregation number.Figure 1h shows the radial distribution functions (rdf s) of water, as well as monomers A and B, for the vesicle depicted in Figure 1g.It can be inferred that the vesicle has a water core radius of approximately 6 nm that does not contain any polymers.The vesicle wall has a diffuse double-layer structure with the rdf s exhibiting two well-defined A peaks on either side of the B-rich inner layer.
Equation ( 12) means beads at the inner boundary are included.while those at the outer boundary are excluded.In the calculation, we used ∆r = 0.1 nm.

Morphology Evolution
Figure 1 illustrates the evolution of morphology during vesiculation in BAB copolymers in the NPT simulation.The structure evolution during the NVT simulation is provided in Figure S1 as the Supplementary Information.The polymers form core (B)-corona (A) micelles that are interconnected by bridges of elongated A segments: see Figure 1a.This structure is unstable due to its relatively large solvent-polymer interfacial area.Hence, it coarsens by the merger between neighboring micelles and, to a lesser extent, the incorporation of small polymer aggregates from the solution (Figure 1b).The coarsening process continues until a cage-like morphology with a circumferential mesh of cylindrical micelles is established (Figure 1c,d).The cylindrical micelles rearrange and coalesce with adjacent structures to form a large lamellar (bilayer) shell with holes (Figure 1e).This structure subsequently evolves into a fully closed bilayer (vesicle) (Figure 1f,g).As seen in Figure 2, the closure of the bilayer is accompanied by a reduction in the volume (or, equivalently, the radius of gyration, Rg) of the structure.In the early stages of structure evolution, the aggregation number, N, grows rapidly due to copolymer aggregation.The step changes in N seen in the later stages of morphology evolution are caused by the addition of small polymer aggregates into the primary structure.Overall, once the cage-like structure shown in Figure 1d is formed, the vesiculation process is facilitated largely through its reorganization at a nearly constant aggregation number.Figure 1h shows the radial distribution functions (rdfs) of water, as well as monomers A and B, for the vesicle depicted in Figure 1g.It can be inferred that the vesicle has a water core radius of approximately 6 nm that does not contain any polymers.The vesicle wall has a diffuse doublelayer structure with the rdfs exhibiting two well-defined A peaks on either side of the Brich inner layer.We studied the changes in the distribution of copolymer chain configurations during the vesiculation process.Figure 3a-g show the pdfs of the end-to-end distance, Q, of the polymer chains within the structures depicted in Figure 1a-g.The probability distribution function, p(Q), for the interconnected micelle morphology shown in Figure 1a exhibits two prominent peaks, namely, the one at the larger Q value (≈4.9 nm), corresponding to the stretched polymers that form the bridges between the micelles, and the second at the lower Q value (≈1.7 nm), representing the folded configurations within the micelles.As the self-assembly process progresses, the peak at the larger Q value becomes smaller and eventually disappears, implying that folded configurations become the predominant ones within the assemblies.There is a very small, albeit persistent, third peak at Q ≈ 0.5 nm, which is close to the cutoff distance for non-bonded interactions in the simulations.This corresponds to a relatively small number of ring-like polymer chains.In Figure 3g, p(Q), obtained from CGMD simulations of a single chain in solution, is also shown for comparison.Polymer chains that constitute the vesicle exhibit a broader Q distribution compared with that of an isolated (single) chain in solution.The mutual affinity of B-type monomer pairs is greater than those of A-B and B-solvent pairs.This allows for the opening of the hairpin-like configurations of the individual polymer chains within the B-rich inner layer of the vesicle, thereby resulting in greater end-to-end distances.We studied the changes in the distribution of copolymer chain configurations during the vesiculation process.Figure 3a-g show the pdfs of the end-to-end distance, Q, of the polymer chains within the structures depicted in Figure 1a-g.The probability distribution function, p(Q), for the interconnected micelle morphology shown in Figure 1a exhibits two prominent peaks, namely, the one at the larger Q value (≈4.9 nm), corresponding to the stretched polymers that form the bridges between the micelles, and the second at the lower Q value (≈1.7 nm), representing the folded configurations within the micelles.As the selfassembly process progresses, the peak at the larger Q value becomes smaller and eventually disappears, implying that folded configurations become the predominant ones within the assemblies.There is a very small, albeit persistent, third peak at Q ≈ 0.5 nm, which is close to the cutoff distance for non-bonded interactions in the simulations.This corresponds to a relatively small number of ring-like polymer chains.In Figure 3g, p(Q), obtained from CGMD simulations of a single chain in solution, is also shown for comparison.Polymer chains that constitute the vesicle exhibit a broader Q distribution compared with that of an isolated (single) chain in solution.The mutual affinity of B-type monomer pairs is greater than those of A-B and B-solvent pairs.This allows for the opening of the hairpin-like configurations of the individual polymer chains within the B-rich inner layer of the vesicle, thereby resulting in greater end-to-end distances.

Interfacial Area and Information Entropy
The variations in the SASA and H, normalized with respect to their initial values, accompanying vesiculation are shown in Figure 4.The polymer-solvent interfacial area decreases rapidly at the early stages of self-assembly.However, once the lamellar cage is formed at t ≈ 60 ns, the SASA does not change appreciably.The total SASA at equilibrium is approximately 16% of its initial value.As expected, the decrease in the relative SASA is greater for the hydrophobic segments (≈99%) than that for the hydrophilic ones (≈74%).In other words, the self-assembly and reorganization of the copolymers are guided by the motif of almost entirely shielding the hydrophobic monomers from the solvent, i.e., to reduce the unfavorable interfacial energy of the hydrophobic interactions at the expense of gaining curvature energy.A qualitatively similar energetic motif of vesicle formation has been inferred from CGMD simulations of AB copolymer solutions [27].However, the vesiculation pathway in the AB solutions consists of topologically simpler morphologies (spherical/rod-like micelles and rectangular/discoid lamella) compared with those depicted in Figure 1.The complex intermediate structures observed in BAB solutions may be a consequence of the remarkable diversity inherent to the single-chain configurations produced by the selective A-B segmental interactions with the solvent, as well as the relative differences in A-B and A-A pair affinities in triblock systems, as illustrated by Figure S2 in the Supplementary Information.

Interfacial Area and Information Entropy
The variations in the SASA and H, normalized with respect to their initial values, accompanying vesiculation are shown in Figure 4.The polymer-solvent interfacial area decreases rapidly at the early stages of self-assembly.However, once the lamellar cage is formed at t ≈ 60 ns, the SASA does not change appreciably.The total SASA at equilibrium is approximately 16% of its initial value.As expected, the decrease in the relative SASA is greater for the hydrophobic segments (≈99%) than that for the hydrophilic ones (≈74%).In other words, the self-assembly and reorganization of the copolymers are guided by the motif of almost entirely shielding the hydrophobic monomers from the solvent, i.e., to reduce the unfavorable interfacial energy of the hydrophobic interactions at the expense of gaining curvature energy.A qualitatively similar energetic motif of vesicle formation has been inferred from CGMD simulations of AB copolymer solutions [27].However, the vesiculation pathway in the AB solutions consists of topologically simpler morphologies (spherical/rod-like micelles and rectangular/discoid lamella) compared with those depicted in Figure 1.The complex intermediate structures observed in BAB solutions may be a consequence of the remarkable diversity inherent to the single-chain configurations produced by the selective A-B segmental interactions with the solvent, as well as the relative differences in A-B and A-A pair affinities in triblock systems, as illustrated by Figure S2 in the Supplementary Information.The variations in H along the vesiculation pathway (≈9%) are modest compared with those in the SASA.When the cage-like structure (Figure 1d) is formed from the interconnected micelles (Figure 1a), the number of stretched chains is significantly reduced.This decreases H, as seen in Figure 4.While a cage-like structure (≈40 ns) reorganizes into a vesicle (≈150 ns), H increases slightly.This can be attributed to the opening of the folded chain configurations due to increased attractive B-B pair interactions within the bilayer.The small decrease in H for t > 150 ns is due to the rearrangement of the copolymers within the closed bilayer forming a nearly spherical and compact vesicle.During this process, R g also decreases by a very small amount (from 11.7 nm to 11.3 nm) (Figure 2).

Effect of Solvent Affinity of the Hydrophobic Segment on Self-Assembly
We further investigated the effect of the variation in hydrophobicity of the B block on the self-assembly mechanism.In the above simulations, bead type C4 from the MARTINI database was used to model the B monomers.The set of force field parameters for this bead type, provided in Table S1 in the Supplementary Information, mimics a polybutadiene monomer.As mentioned in Section 2, the bead type SNda used for segment A models a polyethylene oxide monomer.To simulate variations in solvent affinity, we keep the bead type of the A segment unchanged and change that of the B segment.Specifically, we perform a set of simulations of self-assembly in 5B10A5B copolymer solutions in which the A monomer has the characteristics of type SNda and the B monomer is modeled by the bead type Cn, n = 1, 2, 3, 4, 5.The hydrophilicity of the Cn bead type is in the following order: C1 < C2 < C3 = C4 < C5.The difference between bead types C3 and C4 is that the latter is more attractive to A beads.The relevant force field parameters of the Cn beads are provided in Table 1.Simulations showed that the conclusions reached regarding the pathway of morphology evolution, as well as variations in SASA and H presented in Section 3.1, are qualitatively unchanged for bead types C1-C4; see Figure S3 in the SI for an illustration of the self-assembly process for the C1-type copolymer, as well as Figure 5 for the variations in N, SASA, and H for the various bead-types.However, the kinetics of self-assembly and The variations in H along the vesiculation pathway (≈9%) are modest compared with those in the SASA.When the cage-like structure (Figure 1d) is formed from the interconnected micelles (Figure 1a), the number of stretched chains is significantly reduced.This decreases H, as seen in Figure 4.While a cage-like structure (≈40 ns) reorganizes into a vesicle (≈150 ns), H increases slightly.This can be attributed to the opening of the folded chain configurations due to increased attractive B-B pair interactions within the bilayer.The small decrease in H for t > 150 ns is due to the rearrangement of the copolymers within the closed bilayer forming a nearly spherical and compact vesicle.During this process, R g also decreases by a very small amount (from 11.7 nm to 11.3 nm) (Figure 2).

Effect of Solvent Affinity of the Hydrophobic Segment on Self-Assembly
We further investigated the effect of the variation in hydrophobicity of the B block on the self-assembly mechanism.In the above simulations, bead type C4 from the MARTINI database was used to model the B monomers.The set of force field parameters for this bead type, provided in Table S1 in the Supplementary Information, mimics a polybutadiene monomer.As mentioned in Section 2, the bead type SNda used for segment A models a polyethylene oxide monomer.To simulate variations in solvent affinity, we keep the bead type of the A segment unchanged and change that of the B segment.Specifically, we perform a set of simulations of self-assembly in 5B10A5B copolymer solutions in which the A monomer has the characteristics of type SNda and the B monomer is modeled by the bead type Cn, n = 1, 2, 3, 4, 5.The hydrophilicity of the Cn bead type is in the following order: C1 < C2 < C3 = C4 < C5.The difference between bead types C3 and C4 is that the latter is more attractive to A beads.The relevant force field parameters of the Cn beads are provided in Table 1.Simulations showed that the conclusions reached regarding the pathway of morphology evolution, as well as variations in SASA and H presented in Section 3.1, are qualitatively unchanged for bead types C1-C4; see Figure S3 in the SI for an illustration of the self-assembly process for the C1-type copolymer, as well as Figure 5 for the variations in N, SASA, and H for the various bead-types.However, the kinetics of self-assembly and vesiculation are influenced by the hydrophilicity of the B segment: the less hydrophilic the beads, the faster the vesicle formation process.For instance, as seen in Table 1, the vesicle formation time predicted for the copolymers with the C1-type B segment is nearly one-half of that for the C3/C4-type cases.Further, the kinetics of vesiculation for the C3and C4-types are statistically indistinguishable, as evidenced by data provided in Table 1 and Figure 5.
Colloids Interfaces 2024, 8, x FOR PEER REVIEW 10 of 16 vesiculation are influenced by the hydrophilicity of the B segment: the less hydrophilic the beads, the faster the vesicle formation process.For instance, as seen in Table 1, the vesicle formation time predicted for the copolymers with the C1-type B segment is nearly one-half of that for the C3/C4-type cases.Further, the kinetics of vesiculation for the C3and C4-types are statistically indistinguishable, as evidenced by data provided in Table 1 and Figure 5. Figure 6 shows the p(Q) of a single chain in solution at equilibrium and chains within the equilibrium morphologies for copolymers with B monomers of different C1-C5 characteristics.As the hydrophilicity of the B segment increases, the distribution of chain lengths becomes broader, as expected.Moreover, in all cases, the pdfs of the chains within the molecular assemblies are broader and exhibit a peak shift to larger Q values as compared with those of single chains.This is due to the swelling of chains within the B-rich regions of the bilayers.Further, the pdfs show a very good degree of fit to lognormal distributions, as previously reported for AB copolymers [27].Figure 6 shows the p(Q) of a single chain in solution at equilibrium and chains within the equilibrium morphologies for copolymers with B monomers of different C1-C5 characteristics.As the hydrophilicity of the B segment increases, the distribution of chain lengths becomes broader, as expected.Moreover, in all cases, the pdfs of the chains within the molecular assemblies are broader and exhibit a peak shift to larger Q values as compared with those of single chains.This is due to the swelling of chains within the B-rich regions of the bilayers.Further, the pdfs show a very good degree of fit to lognormal distributions, as previously reported for AB copolymers [27].
The BCP architecture with the most hydrophilic B block (C5 beads) does not form vesicles but instead forms multilayered spherical micelles.The primary driving force for self-assembly is hydrophobic interactions that lead to the shielding of the B blocks from the solvent (water).However, morphology selection depends on the relative affinities of the A and B blocks with the solvent and with each other.The enhanced hydrophilicity of the C5 beads makes them significantly more compatible with the A segment.This inhibits the tendency of microphase separation and promotes the formation of multilayered micelles with mixed domains, as evidenced by the rdfs shown in Figure 7h. Figure 7 shows the self-assembly pathway for the C5-type copolymers.Interconnected copolymer aggregates are formed at the early stages of self-assembly (Figure 7a).Subsequently, adjacent micelles aggregate to form short rods (Figure 7b).Up to this stage, the morphology evolution parallels those seen in vesicle-forming BAB solutions (Figure 1 and Figure S3 in the Supplementary Information).However, owing to the enhanced hydrophilicity, the rodlike micelles of the C5-type BCPs merge with adjacent structures to form longer wormlike micelles, with spatial distributions shown in Figure 7c,d.Fragments that break away from the largest (primary) structure cause a decrease in R g .The primary structure further shrinks to become a spherical multilayer micelle (Figure 7e-g).The rdf s plotted in Figure 7h clearly show a B-rich region in the micelle core surrounded by a shell that contains both A and B monomers inside and hydrophilic A segments at the solvent-micelle interface.Figure 6 shows the p(Q) of a single chain in solution at equilibrium and chains within the equilibrium morphologies for copolymers with B monomers of different C1-C5 characteristics.As the hydrophilicity of the B segment increases, the distribution of chain lengths becomes broader, as expected.Moreover, in all cases, the pdfs of the chains within the molecular assemblies are broader and exhibit a peak shift to larger Q values as compared with those of single chains.This is due to the swelling of chains within the B-rich regions of the bilayers.Further, the pdfs show a very good degree of fit to lognormal distributions, as previously reported for AB copolymers [27].The BCP architecture with the most hydrophilic B block (C5 beads) does not form vesicles but instead forms multilayered spherical micelles.The primary driving force for self-assembly is hydrophobic interactions that lead to the shielding of the B blocks from the solvent (water).However, morphology selection depends on the relative affinities of the A and B blocks with the solvent and with each other.The enhanced hydrophilicity of the C5 beads makes them significantly more compatible with the A segment.This inhibits the tendency of microphase separation and promotes the formation of multilayered micelles with mixed domains, as evidenced by the rdfs shown in Figure 7h. Figure 7 shows the self-assembly pathway for the C5-type copolymers.Interconnected copolymer aggregates are formed at the early stages of self-assembly (Figure 7a).Subsequently, adjacent micelles aggregate to form short rods (Figure 7b).Up to this stage, the morphology evolution parallels those seen in vesicle-forming BAB solutions (Figure 1 and Figure S3 in the Supplementary Information).However, owing to the enhanced hydrophilicity, the rod-like micelles of the C5type BCPs merge with adjacent structures to form longer wormlike micelles, with spatial distributions shown in Figure 7c,d.Fragments that break away from the largest (primary) structure cause a decrease in Rg.The primary structure further shrinks to become a spherical multilayer micelle (Figure 7e-g).The rdfs plotted in Figure 7h clearly show a B-rich region in the micelle core surrounded by a shell that contains both A and B monomers inside and hydrophilic A segments at the solvent-micelle interface.

Conclusions
CGMD simulations of the real-time evolution of vesicles-starting with an initial condition of randomly distributed stretched BAB copolymer chains in an explicit solventhave been performed to identify intermediate self-assembled states along the vesiculation pathway.The vesiculation process is facilitated by the rapid formation of interconnected spherical aggregates; the merger of spherical aggregates to form a cage-like assembly made up of cylindrical micelles that reorganize into lamellar (bilayer) structures; and the closure of the lamellar cage to form a vesicle.The SASA decreases monotonically during vesiculation, with a near (99%) elimination of the hydrophobic (B-solvent) interface.The transition from a cage-like morphology to vesicle morphology occurs at a practically constant aggregation number with negligible changes in configuration entropy.Such a vesiculation mechanism facilitated by the isentropic reorganization of an aggregate of a critical size (aggregation number) has also been observed in CGMD simulations of AB copolymer simulations [27].However, the intermediate structures predicted in BAB systems are

Conclusions
CGMD simulations of the real-time evolution of vesicles-starting with an initial condition of randomly distributed stretched BAB 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 the rapid formation of interconnected spherical aggregates; the merger of spherical aggregates to form a cage-like assembly made up of cylindrical micelles that reorganize into lamellar (bilayer) structures; and the closure of the lamellar cage to form a vesicle.The SASA decreases monotonically during vesiculation, with a near (99%) elimination of the hydrophobic (B-solvent) interface.The transition from a cage-like morphology to vesicle morphology occurs at a practically constant aggregation number with negligible changes in configuration entropy.Such a vesiculation mechanism facilitated by the isentropic reorganization of an aggregate of a critical size (aggregation number) has also been observed in CGMD simulations of AB copolymer simulations [27].However, the intermediate structures predicted in BAB systems are topologically complex in comparison with the regular transient shapes (spherical/rod-like, rectangular/discoid lamella [25,27]) seen during vesicle morphogenesis in AB diblock polymers.The endto-end distance of a single chain at equilibrium and that of individual chains within the equilibrated vesicular structure follow lognormal probability distributions.However, the latter distribution is broader, suggesting that the increased hydrophobicity (the increase in the number of proximal B monomer pairs) within the bilayer results in the swelling of the hairpin-shaped BAB copolymer chains.The simulations were used to probe the influence of monomer chemistry on morphology evolution.Increasing the hydrophobicity of the B segments does not influence the self-assembly pathway.However, it makes the kinetics of vesiculation faster.Decreasing the hydrophilicity of the B segments prevents the formation of vesicles and promotes a molecular assembly process that results in multi-layered spherical micelles.

Figure 1 .
Figure 1.Structure evolution: (a) 0 ns (initial condition of the NPT simulation): a network of interconnected spherical aggregates; (b) 10 ns: a network of ellipsoidal micelles; (c) 20 ns; (d) 30 ns: cage of cylindrical micelles; (e) 50 ns: a lamellar cage; (f) 150 ns: precursor to a vesicle; (g) 600 ns: a vesicle at equilibrium.Panel (h) shows the pair correlation functions of the vesicle shown in panel (g).The

Figure 1 .
Figure 1.Structure evolution: (a) 0 ns (initial condition of the NPT simulation): a network of interconnected spherical aggregates; (b) 10 ns: a network of ellipsoidal micelles; (c) 20 ns; (d) 30 ns: cage of cylindrical micelles; (e) 50 ns: a lamellar cage; (f) 150 ns: precursor to a vesicle; (g) 600 ns: a vesicle at equilibrium.Panel (h) shows the pair correlation functions of the vesicle shown in panel (g).The black, blue, and red lines represent water, A, and B, respectively.The smaller illustration next to panel (a) is an exploded view of the micelle network.The smaller illustrations next to panels (b-g) show the cross-section of the corresponding structure.
black, blue, and red lines represent water, A, and B, respectively.The smaller illustration next to panel (a) is an exploded view of the micelle network.The smaller illustrations next to panels (b-g) show the cross-section of the corresponding structure.

Figure 2 .
Figure 2. Structure evolution: aggregation number (N) and radius of gyration (Rg) versus time.The step changes seen in N correspond to the merger of a small copolymer aggregate with the primary structure.For snapshots 4-6, cross-sections are also shown.The dotted line shows the time at which a fully enclosed vesicular structure is formed.

Figure 2 .
Figure 2. Structure evolution: aggregation number (N) and radius of gyration (Rg) versus time.The step changes seen in N correspond to the merger of a small copolymer aggregate with the primary structure.For snapshots 4-6, cross-sections are also shown.The dotted line shows the time at which a fully enclosed vesicular structure is formed.

Figure 3 .
Figure 3.The probability distribution function of the end-to-end distance, Q, of the polymer chains at different times corresponding to the structures (a-g) shown in Figure 1: (a) 0 ns, (b) 10 ns, (c) 20 ns, (d) 30 ns, (e) 50 ns, (f) 150 ns, and (g) 600 ns.Typical chain configurations are shown in the inset of panels (a,g).The black points and the line in panel (g) show p(Q) for a single chain at equilibrium.

Figure 3 .
Figure 3.The probability distribution function of the end-to-end distance, Q, of the polymer chains at different times corresponding to the structures (a-g) shown in Figure 1: (a) 0 ns, (b) 10 ns, (c) 20 ns, (d) 30 ns, (e) 50 ns, (f) 150 ns, and (g) 600 ns.Typical chain configurations are shown in the inset of panels (a,g).The black points and the line in panel (g) show p(Q) for a single chain at equilibrium.

Figure 4 .
Figure 4. Relative solvent-accessible surface area and relative information entropy versus time.The dotted line shows the time at which a fully enclosed vesicular structure is formed.

Figure 4 .
Figure 4. Relative solvent-accessible surface area and relative information entropy versus time.The dotted line shows the time at which a fully enclosed vesicular structure is formed.

Figure 6 .
Figure 6.p(Q) vs. Q for BAB copolymer solutions with varying degrees of hydrophilicity in the B block (C1 < C2 < C3 = C4 < C5) for (a) a single chain in solution and (b) the self-assembled structures

Figure 6 .
Figure 6.p(Q) vs. Q for BAB copolymer solutions with varying degrees of hydrophilicity in the B block (C1 < C2 < C3 = C4 < C5) for (a) a single chain in solution and (b) the self-assembled structures Figure 6.p(Q) vs. Q for BAB copolymer solutions with varying degrees of hydrophilicity in the B block (C1 < C2 < C3 = C4 < C5) for (a) a single chain in solution and (b) the self-assembled structures at equilibrium.Cases C1-C4 form vesicles, whereas C5 forms spherical micelles.Lines represent lognormal fits to the simulation data shown in the symbols.

Figure 7 .
Figure 7. Structure evolution for C5-type B segment: (a) 0 ns, (b) 10 ns, (c) 20 ns, (d) 30 ns, (e) 50 ns, (f) 150 ns, and (g) 600 ns.Panel (h) shows the radial distribution functions of the micelle shown in panel (g).The black, blue, and red lines represent water, A, and B, respectively.The smaller illustrations next to panels (b-g) show the cross-sections of the corresponding large structures.

Figure 7 .
Figure 7. Structure evolution for C5-type B segment: (a) 0 ns, (b) 10 ns, (c) 20 ns, (d) 30 ns, (e) 50 ns, (f) 150 ns, and (g) 600 ns.Panel (h) shows the radial distribution functions of the micelle shown in panel (g).The black, blue, and red lines represent water, A, and B, respectively.The smaller illustrations next to panels (b-g) show the cross-sections of the corresponding large structures.

Table 1 .
Time required for the formation of a fully enclosed vesicular structure for chains with hydrophobic segments with different degrees of hydrophobicity.

Table 1 .
Time required for the formation of a fully enclosed vesicular structure for chains with hydrophobic segments with different degrees of hydrophobicity.