Langevin Dynamics Simulations of the Exchange of Complex Coacervate Core Micelles: The Role of Nonelectrostatic Attraction and Polyelectrolyte Length.

Complex coacervate core micelles (C3Ms) are promising encapsulators for a wide variety of (bio)molecules. To protect and stabilize their cargo, it is essential to control their exchange dynamics. Yet, to date, little is known about the kinetic stability of C3Ms and the dynamic equilibrium of molecular building blocks with micellar species. Here we study the C3M exchange during the initial micellization by using Langevin dynamics simulations. In this way, we show that charge neutral heterocomplexes consisting of multiple building blocks are the primary mediator for exchange. In addition, we show that the kinetic stability of the C3Ms can be tuned not only by the electrostatic interaction but also by the nonelectrostatic attraction between the polyelectrolytes, the polyelectrolyte length ratio, and the overall polyelectrolyte length. These insights offer new rational design guides to aid the development of new C3M encapsulation strategies.

Introduction

Complex coacervate core micelles (C3Ms) have a core that consists of complexes of oppositely charged polyelectrolytes and are therefore well-suited to encapsulate hydrophilic (bio)molecules. The core formation relies on the associative phase separation of the oppositely charged polyelectrolytes from the water phase. The phase rich in polyelectrolytes is called the complex coacervate. Macroscopic phase separation is prevented by a neutral, water-soluble block that is attached to at least one the two polyelectrolyte types. This neutral block forms the corona around the complex coacervate core. Molecules that prefer to go to the complex coacervate phase can be incorporated in the C3M core and can be protected from the outside by the surrounding corona. This makes the C3Ms promising encapsulators for different types of (bio)molecules. In fact, C3Ms have already been studied as encapsulators for many applications,[1] especially for drug and gene delivery.[2−4]

To design good C3M encapsulators, it is essential to understand their exchange dynamics. First, the exchange dynamics between C3Ms determines the rate with which cargo in the core is exposed to the surroundings and thus the level of protection the encapsulation vehicle offers. In addition, the structure of the C3Ms sometimes depends on their preparation pathway,[5−8] which means that kinetic effects can govern the C3M structure and thus their encapsulation properties.

For amphiphilic diblock copolymer micelles the importance of exchange dynamics is widely recognized, and their exchange dynamics has been thoroughly studied.[9−18] These micelles consist of a single macromolecular species featuring an soluble and an insoluble block. Often, two different mechanisms are distinguished to describe the exchange of these micelles, based on the theoretical framework developed by Dormidontova.[11] The first mechanism is unimer exchange. Here, one polymer (or a few polymers) splits off and is inserted into another micelle. The second mechanism is fission followed by fusion. In that case, the micelle splits into two parts of both substantial sizes, which can subsequently fuse with another micelle. For the fission, both parts that are formed still have a corona structure. For the expulsion, however, the expelled part contains only one or two soluble blocks, which is not enough to form a micelle corona. Because of this difference in corona both mechanism have different rate-limiting steps and therefore different dependencies on the system parameters. A change of one of the system parameters can thus change which of the two mechanisms dominates or whether they both occur. As a result, micelle exchange rates can have a complex dependence on system parameters like core block length,[12−16] corona block length,[12,16,17] polymer concentration,[12,14] chain flexibility,[16] and interfacial tension between core and solvent.[12,14,18]

Although the theory developed for amphiphilic diblock copolymer micelle exchange provides a good starting point to describe C3M exchange, it cannot describe the C3M exchange completely. The two micelle types differ in the interactions that drive the core formation. For amphiphilic diblock copolymer micelles, the core formation is usually driven by hydrophobic attraction, while for C3Ms the core formation is the result of electrostatic attraction between the oppositely charged polyelectrolytes enabling the release of counterions. The difference in interaction can be partly accounted for because polyelectrolyte complexes have already been thoroughly studied. For example, it has been shown that the dynamics within the complex coacervate phase depends on both the polyelectrolyte length and the salt concentration.[19] This can help to describe the relaxation within the C3M core. In addition, the interfacial tension of certain complex coacervates[20] and the strength of an ionic bond have been measured,[21] which can help to describe the release of polyelectrolytes from the C3M core. However, the fact that the core formation is based on the attraction between two different block types instead of one also introduces additional tuning parameters that are absent for amphiphilic diblock copolymer micelles. Instead of only varying the length of one core block, the block length of the negative and positive polymer block can be varied independently from each other. Furthermore, the choice to attach the corona block to both core blocks or to only one of the two can alter the micelle properties.[22] It is thus insufficient to focus only on amphiphilic diblock copolymer micelles to understand the exchange of C3Ms.

Unfortunately, the exchange of C3Ms is much less studied, and the few studies done on this subject[23,24] had an indirect way of interpreting their results. The authors mixed C3Ms labeled with a donor fluorophore with C3Ms labeled with an acceptor fluorophore and took the rate of increase in Förster resonance energy transfer (FRET) as a measure for the micelle exchange rates. In this way, they found that the exchange rate depends on the polyelectrolyte length and charge stoichiometry[23] and that the exchange rate is fast for C3Ms containing proteins.[24] Subsequent interpretations on the exchange mechanisms were mainly based on the observed exchange rates and not on direct observations. Although these studies give very useful insights into the time scales at which micelle exchange can occur, they are limited in the mechanistic descriptions they can provide.

In this paper we aim to provide insights into the molecular mechanisms of exchange in C3Ms. To this end, we exploit coarse-grained dynamics simulations. This type of simulations has already shown its value in the studies on amphiphilic diblock copolymer micelle exchange.[12,14−16] In addition, coarse-grained dynamics simulations have also been used to study complex coacervation,[25,26] the formation of a single C3M,[27] and the static properties of multiple C3Ms.[28−30] We use the coarse-grained simulations to follow the initial micellization kinetics of multiple C3Ms. In this way, we obtain new mechanistic insights into the exchange of C3Ms, and we identify ways to improve the kinetic stability of C3Ms. This information about the kinetic stability can complement earlier on the static C3M stability,[31−34] which is usually expressed as the critical salt concentration where the micelles fall apart. In particular, we show that also the kinetic stability of C3Ms can be improved by tuning the nonelectrostatic interactions and the polyelectrolyte length ratio.

Methods

We used the Kremer–Grest bead–spring model to describe flexible chains in a good solvent and included electrostatic interactions and explicit ions to account for the electrostatic nature of the process we aim to describe. In the model, the polymers are represented by multiple beads connected with springs, while the counterions are represented by single beads. The solvent is modeled implicitly. All beads have the same diameter σ and the same mass m. The springs represent polymer bonds and are modeled with a finitely extensible nonlinear elastic (FENE) potential with a bond stiffness k of 30kBT/σ2 and a maximum bond extension distance r of 1.5σ. The electrostatic interactions between the beads are modeled with a Coulomb potential. Unless otherwise stated, we modeled the nonelectrostatic interactions between equally and oppositely charged monomers with a Lennard-Jones potential with a cutoff distance of 2.5σ. We varied the strength of the nonelectrostatic attraction by changing the minimum of the Lennard-Jones potential εLJ in the range 0.05kBT–0.25kBT. In contrast to the nonelectrostatic interaction between equally and oppositely charged monomers, the nonelectrostatic interaction between all other monomer–monomer, monomer–ion, and ion–ion combinations was purely repulsive. For this repulsive interaction we used the Weeks–Chandler–Andersen (WCA) potential with an interaction strength ε = 1kBT. A graphical overview of all used potentials can be found in the Supporting Information (SI1).

This coarse-grained bead–spring representation has often been used to model polyelectrolytes and can be mapped to experimentally realistic systems via the Bjerrum length lB. The Bjerrum length sets the length scale of the electrostatic interaction. In pure water at room temperature the Bjerrum length is 0.71 nm. In our model we used a Bjerrum length of lB = 2.5σ, which means that the bead diameter σ roughly corresponds to ∼0.3 nm. The average polyelectrolyte bond distance is 0.97σ for the simulation parameters that we used. This means that the charge separation distance in our simulation is ∼2.8 Å, which is close to the distance between two adjacent side groups of a polymer carbon backbone. We thus simulate polyelectrolytes where every side group is charged, such as the strongly charged polyelectrolyte polystyrenesulfonate.

The C3Ms in our simulations are formed from a combination of coarse-grained homopolymers, diblocks, and counterions. The homopolymers are negatively charged and have a length Nneg, where N represents the number of monomers. The diblock consists of a positively charged block with length Npos and a neutral block with length Nneu. We varied the lengths of the negative polymer and the diblock polymer between the different simulations, but we chose the parameters such that the number of negative monomers was always equal to the number of positive monomers. In addition, we fixed the total number of charged monomers at 24000, and we kept the ratio between the positive and neutral block length at 2:5. Similar ratios have also been used in experimental studies of C3Ms.[34,35] In our simulations this ratio ensured that the neutral block was long enough to prevent macrophase separation and on the other hand short enough to allow the formation of micelles instead of only free soluble complexes. We added only counterions to the simulation box; no additional salt ions were added. We note that in solvents of experimental systems always some additional salt ions are present. However, since we use a periodic box size of L = 235σ every time, the counterion concentration is relatively large compared to, for example, the ion concentration of distilled water, and these solvent ions thus can be neglected.

We started our simulation by placing the homopolymers, diblocks, and counterions randomly in the simulation box and then used Langevin dynamics simulations to simulate the formation of the C3Ms in time. We used γ = 1m/τ as drag coefficient and Δt = 0.005τ as simulation time step where is the time unit in the system. We saved the configuration of the simulation every 2500 steps.

To perform the simulations, we used the GPU-optimized molecular dynamics software package HOOMD-Blue.[36−39] The Coulomb interactions were calculated by using the particle–particle particle-mesh (PPPM) Ewald summation method[38] with a real space cutoff distance of 2.5σ. The neighbor lists were generated by using the linear bounding volumes hierarchies (LBVHs) method.[39] We used visual molecular dynamics (VMD)[40] to visualize the simulations.

To analyze the simulation data, we first identified the micelles with the data clustering algorithm DBSCAN. In particular, we used the algorithm as implemented in the Python package scikit-learn[41] with a maximum allowable neighborhood radius of 2σ and a minimum neighborhood points number of three for a point to be a core point. To avoid that two near micelles were identified as a single micelle, we based the clustering algorithm on the polyelectrolyte coordinates and did not take the neutral block into account. We provided a precomputed sparse array as neighbor array for the DBSCAN algorithm. To obtain this array, we used the KDTree neighbor algorithm from scikit-learn. Because this algorithm does not take the periodic boundaries into account, we first added the surrounding periodic boundary images, used the KDTree algorithm to calculate neighbor list for the original simulation box and its periodic images together, and then converted this to a periodic neighbor list for the original simulation box. We performed the micelle cluster identification for every tenth saved configuration (125τ). The intermediate saved configurations were analyzed if the micelle composition changed within these ten steps.

The clustering algorithm yielded the micelle compositions for every time step, and this was used to analyze the micelle exchange. A decrease in micelle size was counted as a split event and an increase in micelle size was counted as a merge event. The discrimination between expulsion and fission and between insertion and fusion was based on the size of the smallest cluster involved in the exchange event: if this cluster contained fewer than five polyelectrolytes, the cluster contained too little diblocks to form a corona structure, and therefore the split event was called expulsion and the merge event was called insertion. Otherwise, the split event was called fission and the merge event was called fusion.

Results and Discussion

C3M Formation and C3M Exchange Mechanisms

Complex coacervate core micelles rapidly form when we mix coarse-grained homopolymers, diblocks, and counterions together in the simulation box (Figure ). First, the polyelectrolytes have a relatively stretched configuration and are surrounded by their counterions. The oppositely charged polyelectrolytes rapidly form complexes upon which they decrease in size and release their counterions (Supporting Information SI 2). At the end of the simulation, we observe clear C3M structures, where the positive and negative polymers together form the core of the micelle, while the neutral blocks form the surrounding corona (Figure b and Supporting Information SI 3). The initial assembly is particularly fast. Afterward the micelle growth levels off, and at the end of the simulation the average micelle size fluctuates around the same value (Figure c).

Figure 1

Initial micellization kinetics of C3Ms for Nneg = 20, Npos = 20, Nneu = 50, and εLJ = 0.15kBT. Snapshots of the begin (a) and the end (b) of the simulation. Homopolymers are depicted blue, the positive block in red, and the neutral block in gray. Counterions are not shown. (c) Increase of the average micelle size, expressed as the average aggregation number of positive blocks per micelle Nagg,pos, in time. (d) Histograms of the C3M size distribution, expressed as total number of polymers per micelle Nagg, at different time points during the simulations.

Although the micelle growth approaches a plateau at the end of the simulation, the micelles probably do not reach complete equilibrium yet. At the end of the simulation, the micelles still have a broad size distribution that is not centered around one optimum value (Figure d). A broad size distribution has been observed for C3Ms at large salt concentrations, but at low salt concentrations these C3Ms are more monodisperse.[34,42] At larger salt concentrations the C3Ms can form wormlike micelles where the length of the micelles can be easily varied without large changes in the free energy. At low salt concentrations the C3Ms form spherical micelles. For the spherical micelles, a change in aggregation number changes the free energy of the micelle. In equilibrium, the micelles will adapt their most favorable configuration, and the micelles thus will have sizes centered around the size with the lowest free energy. Because in our simulations we did not add additional salt ions and the micelles are spherical, we expect that the equilibrium C3M size distribution will be centered around one optimal size. At the end of our simulations, however, we still observe a broad size distribution with multiple maxima. Simulations of amphiphilic diblock copolymer micelles have shown that the broad size distribution can indeed evolve to a clear bimodal size distribution after longer simulation times with one size corresponding to the unimers and one to the micelles.[12] However, C3Ms sometimes need days to fully equilibrate,[8] which would take far too long to simulate with Langevin dynamics. Therefore, we focus here on the initial micellization kinetics of C3Ms and not on their equilibrium dynamics.

In addition to following the initial average micelle growth, the simulations also allow us to directly follow the individual micelle exchange (Figure ). We can observe both expulsion and insertion events and fission and fusion events. The small time and length scales make it impossible to directly observe these events in real experiments. Our simulations thus complement the experiments and can give a deeper insight into the mechanisms that underlie the C3M exchange in the early micellization stage. We note that the C3M exchange mechanisms can be different in a later stage, as is the case for amphiphilic diblock copolymer micelles.[12,43] However, since the interactions between the different monomers remain the same, our observations on the initial exchange can still help to better understand the equilibrium exchange.

Figure 2

Two mechanisms by which micelle exchange occurs in the simulations. (a) Only one or a few polyelectrolytes are expelled from the micelle and inserted into another micelle. (b) The micelle splits in two parts of both substantial sizes (fission), which can combine with other micelles to form a new micelle (fusion). Images are snapshots from simulations with the same color coding as in Figure . Counterions are not shown.

The C3M expulsion mechanism (Figure a) is slightly different from the one of amphiphilic diblock copolymer micelles, and we ascribe this difference to the difference in core interactions. In amphiphilic diblock copolymer micelles, the unfavorable interactions of the core block with the solvent drives the core formation. For C3Ms, however, the core formation is mainly driven by the electrostatic attraction between the oppositely charged core blocks. In particular, strongly and oppositely charged polyelectrolytes form complexes mainly because the entropy increases due to the release of counterions.[26,44] This entropy increase is less when a part of the charge of the polyelectrolytes is not compensated by the oppositely charged polyelectrolytes, since in that case fewer counterions are released. Neutral complexes are thus preferred. As a result, we observe that small neutral complex are expelled from C3Ms instead of the unimers that are usually expelled from amphiphilic micelles. This expulsion of neutral complexes was already predicted to describe the exchange kinetics of C3Ms loaded with fluorescent proteins[24] and is now confirmed by our simulations.

Effect of Nonelectrostatic Attraction Strength

Although the electrostatics plays the most important role in the formation of polyelectrolyte complexes, the nonelectrostatic interactions can also have an effect. The critical salt concentration and binding strength of polyelectrolyte complexes strongly depend on the polyelectrolytes that are used.[45,46] Also for C3Ms both the critical salt concentration and their structure strongly depend on the type of polyelectrolyte.[28,34,47] This shows that apart from the number of charges on the polyelectrolyte, the polymer chemistry also plays a role in polyelectrolyte complexation. Examples of factors that increase the nonelectrostatic attraction between the polyelectrolytes are the hydrophobicity of the polymers and the presence of amine groups or aromatic groups.[46]

To test whether the nonelectrostatic attraction also affects the exchange kinetics of the C3Ms, we varied the nonelectrostatic attraction between the polyelectrolytes by varying εLJ. An increase in εLJ results in a decrease in the repulsion between likely charged monomers at distances where the repulsion energy is close to the thermal energy (Supporting Information SI 1). For the oppositely charged monomers the change in εLJ has a more complex effect, and it also changes the location of the minimum of the combined Lennard-Jones and Coulomb potential (SI 1). The largest nonelectrostatic attraction strength that we used was εLJ = 0.25kBT. This is still a weak attraction and still corresponds to a polymer in good solvent: for the Lennard-Jones bead–spring model the transition to poor solvent conditions occurs at εLJ = 0.33kBT.[48] Even for the largest nonelectrostatic attraction strength in our simulations, the electrostatics thus remains the main driving force to form micelles and not the nonelectrostatic attraction: when the electrostatic interactions are turned off, the micelles fall apart (SI 4). In this way we ensured that we specifically studied C3Ms instead of repeating the studies on amphiphilic micelles where nonelectrostatic attraction drives micelle formation.

Even a small increase in the nonelectrostatic attraction largely decreases the number of exchange events, especially in the later stage of the micelle formation (Figure ). The first 104τ, the insertion, and fusion rates are still comparable for the different nonelectrostatic attraction strengths (SI 5). The nonelectrostatic attraction thus does not play a large role in the early assembly. In this stage, the merge events strongly outnumber the split events. This indicates rapid micelle growth, which we indeed observed in Figure c. The fusion outnumbers the fission longer than the insertion events outnumber the expulsion events. The fast early assembly stage is thus followed by a stage where the micelle growth occurs at the expense of smaller micelles, while the dimer population remains approximately constant (as also shown in SI 2.2–2.3). In this stage, the majority of the chains in the micelles originates from fusion events for εLJ = 0.15kBT and εLJ = 0.25kBT, while for εLJ = 0.05kBT insertion remains the dominant mechanism of micelle growth (SI 6). The situation of the stronger nonelectrostatic attractions resembles to what has been earlier observed in simulations for amphiphilic diblock copolymer micelles: these amphiphilic micelles also had a stage where their growth was mainly governed by fusion of small aggregates.[12] At the end of our simulations, both the insertion rate is similar to the expulsion rate and the fusion rate is similar to the fission rate, indicating a slow micelle growth. In this stage, the occurrence of all exchange processes depend strongly on the nonelectrostatic attraction.

Figure 3

Cumulative number of insertion (a), expulsion (b), fusion (c), and fission (d) events for a nonelectrostatic attraction strength between the polyelectrolytes of εLJ = 0.05kBT, εLJ = 0.15kBT, and εLJ = 0.25kBT. In all cases, Nneg = 20, Npos = 20, and Nneu = 50.

For the parameter set of Figure , the expulsion of dimers occurs often, which indicates that dimers can be easily formed. At the end of the simulation, we indeed observe a large population of dimers, which decreases with increasing nonelectrostatic attraction (Figure ). A similar decrease in dimer population with increasing nonelectrostatic attraction was observed by Šindelka et al.,[29] although they more strongly increased the nonelectrostatic attraction and directly went from good to poor solvent conditions. Dimers can be formed because the enthalpic penalty to expel a neutrally charged dimer from a complex coacervate is small and can be counterbalanced by the entropy that the dimer gains when it is expelled from the complex coacervate. For low polyelectrolyte concentrations complex coacervation even does not occur, and instead only neutral globules are formed.[49,50] As already pointed out by Šindelka et al.,[28] a stronger nonelectrostatic attraction introduces another enthalpic attraction that counteracts the entropy increase when the dimer is released. In this way, the dimer formation is diminished.

Figure 4

Snapshots (a–c) and histograms (d–f) of the C3M size distribution at the end of the simulation for a nonelectrostatic attraction strength between the polyelectrolytes of εLJ = 0.05kBT (a, d), εLJ = 0.15kBT (b, e), and εLJ = 0.25kBT (c, f). In all cases, Nneg = 20, Npos = 20, and Nneu = 50. The color coding of the simulation snapshots is the same as in Figure . Counterions are not shown.

Up to now, we have concluded only that a stronger nonelectrostatic attraction decreases both the insertion/expulsion and fission/fusion events, but we have not yet further quantified this decrease. To make a more quantitative description, we have plotted the number of split events in the time range 2 × 106τ–5 × 106τ versus the nonelectrostatic attraction strength (Figure ). We selected the split events in the last part of the simulations because these events are mainly caused by the micelle exchange kinetics. The merge events and the early split events are affected by both the micelle exchange kinetics and the micelle growth. The simultaneous occurrence of two different processes complicates the analysis, and therefore we decided to focus on the later split events.

Figure 5

Effect of the nonelectrostatic attraction strength on the expulsion events (filled symbols) and fission events (open symbols) in the time range 2 × 106τ–5 × 106τ, for Npos = 20, Nneg = 20, and Nneu = 50 and for Npos = 30, Nneg = 30, and Nneu = 75. Data points are the average of two simulations.

The number of expulsion events seems to decrease exponentially with increasing nonelectrostatic attraction for a polyelectrolyte length of 20. This suggests a thermally activated process where the breaking of nonelectrostatic attraction contributes to the energy barrier. The rate k of a thermally activated process can be described by the Arrhenius equation. To break the interactions of n monomers, the Arrhenius equation is given by k = A exp(−nEa/kBT), where Ea is the activation energy to break the interactions of a single monomer and A is a constant. The total number of core block monomers is 40 for dimers with Npos = Nneg = 20. In the simplest description, all these monomers contribute to the activation energy and the nonelectrostatic activation energy per monomer is just given by εLJ. This is not the case: the slope of the natural logarithm of the number of expulsion events versus the nonelectrostatic attraction strength εLJ is −25.3 ± 0.51 and not −40.

Two factors can contribute to the difference between the observed slope and the slope of the simplified description. First, the activation energy per monomer is not given by εLJ but first has to be multiplied with a numerical prefactor. Not every monomer is positioned from one other monomer at exactly the distance of the minimum of the Lennard-Jones potential. The distance between the monomers can deviate, and a monomer might also have interactions with more than one monomer. To correct for this, the numerical prefactor is needed. This numerical prefactor was also used to describe the equilibrium exchange of amphiphilic diblock copolymer micelles.[14] Second, some of the monomers might not contribute to the activation energy barrier. For example, if the dimer is expelled as a compact globule of N monomers, only the outer N2/3 monomers will contribute to the activation energy barrier (in the Arrhenius equation n = N2/3). In fact, simulation snapshots of free dimers show a configuration in between a linear chain and a compact globule (SI 7). This would mean that N monomers contribute to the energy barrier with the exponent a in between 2/3 and 1.

Increasing the polyelectrolyte lengths to 30 largely decreases the number of expulsion events. For the larger nonelectrostatic attractions, the number of expulsion events in the simulations even becomes too low for reliable statistics. We note that for εLJ = 0.05kBT the decrease in expulsion rate is much larger than we expect for a thermally activated process based on only nonelectrostatic attraction. For a thermally activated process, the maximum decrease based on only nonelectrostatic attraction occurs when all monomers N contribute to the activation energy barrier. In that case, based on the data for Npos = Nneg = 20, the expulsion rate would depend on N as k = A exp(−0.63NεLJ). For an increase of both polyelectrolyte lengths from 20 to 30, the expected decrease factor is thus k/k = exp(−0.63·40·εLJ/kBT)/exp(−0.63·60·εLJ/kBT) = 2. In the simulations, however, the expulsion rate is more than 60 times decreased. This much larger decrease could mean that the dimer expulsion is not a simple thermally activated process, although the expulsion rate seems to decrease exponentially with increasing εLJ for the polyelectrolyte lengths of 20. Alternatively, this larger decrease could also mean that apart from the nonelectrostatic attraction also additional factors play a role in the dimer expulsion.

Other factors that can affect the dimer expulsion are the corona block and the electrostatic interactions. In simulations of amphiphlic diblock copolymer micelles the expulsion rate slightly decreased when the corona block:core block length ratio was increased.[12,16] The change in micelle exchange rates was ascribed to the change in micelle aggregation number that occurred by increasing the corona:core block length ratio. Here we tried to minimize the effect of the corona block by keeping the corona:core block length ratio fixed. If the corona block length has any effect for this fixed ratio, we would expect that a longer corona block increases the expulsion rate: when the dimer is expelled, the corona block gains entropy because it does not longer have to be in a stretched configuration. This entropy increase will be larger for longer corona blocks. The expected expulsion rate increase is opposite to the large decrease that we observe in simulations. This suggests that the corona block has a minor effect on the dimer expulsion. If the expulsion is indeed a thermally activated process, the large decrease thus has to be ascribed to the electrostatic interactions. The electrostatic interactions might affect the expulsion rate because the electrostatic bonds might first need to rearrange before a neutral dimer can be expelled. A larger number of monomers requires that more electrostatic bonds are rearranged and thus decreases the expulsion rate.

So far we have mainly described the expulsion. The fission requires a slightly adapted description. For example, the fission rate also decreases with increasing nonelectrostatic attraction, but this decrease is smaller than for the expulsion rate. We ascribe this difference to a larger change in surface energy for the expulsion compared to fission. For both split events, the total surface increases, but for expulsion this decrease will be larger because the expelled dimer has a large surface-to-volume ratio. Effectively, more nonelectrostatic bonds thus have to be broken for an expulsion events, resulting in a stronger dependence on nonelectrostatic attraction.

The different dependencies on nonelectrostatic attraction seems to result in a transition from expulsion dominated exchange to fission dominated exchange in our simulations. This transition is better visible if we correct the number of split events for the number of polymers that is involved per split event (SI8). An insertion/expulsion dominant exchange was expected by Nolles et al.,[24] who studied the exchange of C3Ms loaded with fluorescent proteins. They expected that the dense corona of the C3Ms would prevent their merging. Here we see that at least in the initial micellization stage micelle fusion can occur. For longer polyelectrolytes with a relatively large nonelectrostatic attraction the fusion/fission is even the dominant exchange mechanism. However, the total number of exchange events is low in our simulations of the long polyelectrolytes with larger nonelectrostatic attraction. More events are needed to confirm that fusion/fission is really the dominant mechanism in these cases. In addition, in this initial micellization stage many small aggregates are present, which fuse more easily than the large micelles that are mainly present in equilibrium.[12,14] Domination of the fission/fusion exchange in the initial micellization period thus does not necessarily mean that this exchange mechanism is also dominant once the micelles have reached equilibrium.

Effect of Asymmetric Nonelectrostatic Attraction

Up to now, we have assumed that the nonelectrostatic attraction is the same for all charged monomers, but this is usually not the case. The two polyelectrolyte types in the C3M always have a different polymer chemistry; otherwise, they cannot be oppositely charged. This different chemistry will result in asymmetric nonelectrostatic interactions: the nonelectrostatic attraction between the negative polyelectrolytes can be different from the one between the positive polyelectrolytes. In addition, the nonelectrostatic attraction between like-charged polyelectrolytes can differ from the nonelectrostatic attraction between oppositely charged polyelectrolytes. Recently, it has been suggested that differences in the latter asymmetry might explain how the stability of C3Ms depends on the positive polyelectrolyte block: Marras et al.[34] found that micelles formed from DNA and a poly(ethylene glycole)–poly(lysine) were more stable than the micelles where the poly(lysine) was replaced by the more hydrophobic poly((vinylbenzyl)trimethylammonium). One of their explanations was that poly(lysine) might form hydrogen bonds with the DNA, while poly((vinylbenzyl)trimethylammonium) might mainly have nonelectrostatic interaction with itself and not with the DNA. This would mean that nonelectrostatic attraction between only the like-charged polyelectrolytes stabilizes the C3Ms less than when also the oppositely charged polyelectrolytes nonelectrostatically attract each other.

We adapted our simulations to test whether a lower nonelectrostatic attraction between the oppositely charged polyelectrolytes indeed results in a lower C3M stability. In addition, we aimed to see how this lower attraction affects the C3M exchange. In the adapted simulation only the positive monomers nonelectrostatically attracted only the other positive monomers, while all other nonelectrostatic interactions between the charged monomers were purely repulsive.

The C3Ms are largely destabilized when they lack nonelectrostatic attraction between the oppositely charged polyelectrolytes and the negative polyelectrolytes (Figure ). In this asymmetric case, both the dimer population and the exchange rate are larger than for C3Ms with a completely symmetric nonelectrostatic attraction of εLJ = 0.05kBT, even though the nonelectrostatic attraction between the positive polyelectrolytes was relatively large in the asymmetric case (εLJ = 0.25kBT). A nonelectrostatic attraction between only one of the polyelectrolytes is thus insufficient to stabilize the C3Ms, and this might indeed explain why the micelles with the hydrophobic poly((vinylbenzyl)trimethylammonium) were less stable than the micelles with poly(lysine). In addition, this shows that it is insufficient to consider only the properties of the individual polyelectrolytes to design stable C3Ms. Instead, also the polyelectrolyte ability to nonelectrostatically interact with the other polyelectrolyte has to be taken into account.

Figure 6

Effect of nonelectrostatic attraction between only the positive polyelectrolytes, while the other polyelectrolyte nonelectrostatic attraction is purely repulsive. In all cases, Nneg = 20, Npos = 20, and Nneu = 50. (a) Snapshot and (b) histogram of the micelle size distribution at the end of the simulation for a nonelectrostatic positive polyelectrolyte attraction of εLJ = 0.25kBT. (c) Cumulative number of expulsion and (d) fission events for the same nonelectrostatic attraction between all polyelectrolytes of εLJ = 0.05kBT or a nonelectrostatic attraction only between the positive polyelectrolytes of εLJ = 0.25kBT.

Effect of Polyelectrolyte Length Asymmetry

Above we have described how the extraction of charge neutral complexes of few molecules is the main mechanism of exchange. This is facilitated when both homopolymer and charged block of the diblock are of the same length. When the block lengths are incommensurate, extracting a strictly charge neutral complex is challenging as it requires a much larger number of molecules to exit the micelle simultaneously. We thus explore how the block-length asymmetry can be used as an additional handle to tune the kinetic micelle stability. This information can complement earlier experimental studies that have shown that the length of the homopolymer (or equivalent) affects the static stability of C3Ms.[32,33]

A small change in the polyelectrolyte length ratio largely affects the expulsion rate (Figure a–c). The expulsion rate is the largest when both polyelectrolytes have an equal length and largely decreases when the length of only the negative homopolymer is changed to give an polyelectrolyte length ratio that differs from 1. For Nneg/Npos = 0.75 and Nneg/Npos = 1.25, charged complexes split off (SI 3). This introduces an additional free energy penalty compared to the neutral complexes that can be formed for Nneg/Npos = 1. As a result, the expulsion rate has a maximum at equal polyelectrolyte lengths. The maximum is the clearest for a small nonelectrostatic attraction. This shows that the net charge of the expelled dimer is most important when no additional nonelectrostatic attraction helps to prevent the dimer formation.

Figure 7

Effect of the polyelectrolyte length ratio Nneg/Nneg on the expulsion (a–c) and fission (d–f) events for a nonelectrostatic attraction strength of εLJ = 0.05kBT (a, d), εLJ = 0.15kBT (b, e), and εLJ = 0.25kBT (c, f). The different polyelectrolyte length ratios were obtained by fixing the positive diblock length at Npos = 20 and varying the negative homopolymer length.

The fission is also affected by a change in the negative homopolymer length (Figure d–f), but in a different way than the expulsion. The fission rate seem to decrease with increasing negative homopolymer length, instead of having a maximum at Nneg/Npos = 1. Only for εLJ = 0.25kBT is this trend no longer visible, which is probably due to the low number number of fission events that occurred at this nonelectrostatic attraction strength. The decrease in fission rate with increasing negative homopolymer length indicates that for fission the total length of the polyelectrolytes is more important than the length ratio. We can explain this by the fact that we used relatively small variations in this ratio. Therefore, neutral complexes can still be formed when multiple polyelectrolytes are combined together, which happens during fission. For shorter polyelectrolytes, the rearrangement of all the noncovalent bonds of one polyelectrolyte is easier. As a result, we observe the fastest fission rate for the shortest negative homopolymers. For low nonelectrostatic attraction, the increase in fission rate for decreasing the negative polymer length to 15 is not enough to compensate for the simultaneous decrease in expulsion rate. This shows that sometimes a counterintuitive situation can occur where a decrease of the polyelectrolyte length results in a decrease of the micelle exchange.

A small length imbalance already has a large effect on the exchange, and this effect becomes even larger when we change the length ratio further to Npos/Nneg = 20/150 (Figure ). In particular, the expulsion rate is largely diminished. For the long homopolymer length of Nneg = 150, the expulsion can occur only when a complex with a large net charge is formed. These uncompensated charges are unfavorable, and therefore the expulsion events nearly disappear. The penalty for the formation of largely charged complexes is also reflected in the micelle size distribution (Figure a,b) where no small complexes can be observed. On the other hand, still many of the micelle sizes that we observe have a net charge (SI 3.2), since only the aggregation numbers that are a multiple of 17 correspond to neutrally charged micelles. Also, the complexes that split off all had a net charge (Figure d). The formation of charged complexes is thus governed not only by the absolute number of charges but also by the number of polyelectrolytes over which these charges can be distributed. If the charges can be distributed over more polyelectrolytes, the formation seems to be easier. This results in the formation of larger charged complexes and in a larger contribution of fission than of expulsion.

Figure 8

Micelle size and exchange characteristics for Npos = 20, Nneg = 150, and Nneu = 50 and εLJ = 0.15kBT. (a) Snapshot and (b) histogram of the micelle size distribution at the end of the simulation. (c) Cumulative number of expulsion and fission events. (d) Histogram of the net charge of the smallest complex formed in a split event.

Finally, we note that the slow expulsion for the long homopolymer length shows that only the exchange times are insufficient to determine the exchange mechanisms. This way of interpretation has, for example, been used by Holappa et al.[23] They measured the Förster resonance energy transfer (FRET) after mixing C3Ms labeled with an acceptor or a donor fluorophore for polyelectrolyte length ratios of 4.3 and 1.8. They ascribed the fast increase in FRET signal to expulsion/insertion and the slow increase to fusion/fission. However, here we observe that for a long homopolymer the expulsion/insertion is initially slower than the fission/fusion instead of the other way round. Although we expect the fusion to slow down a bit once equilibrium has been reached due to the decrease of the number of small aggregates, we do not know yet whether this decrease is large enough to become slower than the expulsion. To check whether the slowest exchange rate observed in the experiments by Holappa et al. indeed corresponds to fission/fusion, additional experiments has to be performed, for example, testing how the exchange rate depends on concentration.[14]

Conclusions

In conclusion, we have directly observed the C3M exchange by using Langevin dynamics simulations, and we have shown that the formation of neutral complexes plays an essential role. For polyelectrolytes of equal length, neutral dimers are expelled. Although the electrostatic attraction plays the major role, also the nonelectrostatic attraction between the polyelectrolytes can be used to tune the exchange. The dimer exchange can be largely diminished by slightly increasing the nonelectrostatic attraction between both polyelectrolytes. Also, the fission of the polyelectrolytes can be diminished in this way, but this decrease will be smaller. An increase in nonelectrostatic attraction between both polyelectrolytes is thus most effective for the insertion/expulsion dominated exchange. It is essential to realize that increasing the nonelectrostatic attraction between only one of the two polyelectrolytes is not effective to make the C3Ms more kinetically stable. Rather than solely focusing on the single polyelectrolyte properties, new studies should thus also pay attention to the interaction between the two polyelectrolytes. Another way to tune the C3M exchange is by changing the length of the polyelectrolytes. The expulsion rate can be easily decreased by using oppositely polyelectrolytes of unequal length, especially by making one of the polyelectrolytes very long. The fission rate can be decreased by increasing the total polyelectrolyte length. These insights into the C3M exchange can be used to develop new C3M encapsulators.