Chemical Feedback in Templated Reaction-Assembly of Polyelectrolyte Complex Micelles: A Molecular Simulation Study of the Kinetics and Clustering

1. Introduction

Polyelectrolyte (PE) complex micelles (PCM) which are formed in aqueous solutions by mixing charged ionic-neutral double hydrophilic diblock copolymers with oppositely charged homopolymers have recently attracted much attention [1,2,3,4,5,6] due to their potential use as nanocarriers and nanoreactors [7,8]. The oppositely charged moieties in the mixture form the PCM core, since the solvent conditions for these moieties worsen due to charge neutralization [9]. The PCM corona is formed of the neutral hydrophilic block. Although the existence of opposite charges on the PE moieties is the necessary condition for electrostatic complexation to occur, the contribution of electrostatics on the Gibbs free energy of mixing becomes significant only at high ionic strengths [10]. In sharp contrast, at low ionic strengths, the driving force of complexation is the entropy gain from the release of a large number of small counterions in the solution [10].

PCMs can be experimentally synthesized in two different ways. The most common is by using pre-synthesized building blocks and subsequent co-assembly to form PCMs. The other synthetic pathway is the polymerization induced electrostatic self-assembly (PIESA). In PIESA [11,12], the polymerization of the charged polyelectrolyte block of the copolymer is templated by the oppositely charged homopolymer chain with simultaneous co-assembly in one-pot reaction. The spatial and time colocalization of covalent and supramolecular electrostatic assembly, involving the same molecular compounds, results in chemical feedback between the different primary reactions, i.e., the polymerization in the solution and the polymerization on the template (Figure 1). Chemical feedback in the coupled reactions-assembly has profound effects on both the kinetics of polymerization and the final size and shape of PCMs.

Bos et al. [13] studied the effects of chemical feedback on the kinetics of the reversible addition-fragmentation chain-transfer (RAFT) polymerization of positive charged monomers of vinylbenzyltrimethylammonium cloride (VBTAC) with PEG chain-transfer agent containing 220 monomers, and poly(sulfopropylmethacrylate) (PSPMA) as the negative template with degree of polymerization (DP) equal to 47. The target length of the VBTAC block in the (PEG)-b-(VBTAC) copolymer was 50, with a monomer to template sites ratio of 1:1. The authors showed that the polymerization rate is strongly enhanced when the template is used at a solution concentration [Φ]≈0.03. This is because the binding of the charged monomer to the template increases the local monomer concentration near the template compared to the monomer concentration in the solution in the non-templated polymerization. The increase of [Φ] to 0.09 was shown to not further affect the polymerization rate. In the same study, the average micelle size in the templated polymerization was found to be smaller than the size of micelles formed by the addition of the template after the polymerization of the charged block of the copolymer.

In a similar experimental study, Ding et al. [14] reported the exact opposite results of Bos et al. but at higher [Φ]≈0.16 and 0.5. Ding et al. performed RAFT polymerization of 2-acrylamydo-2-methylpropanesulfonic (AMPS) acid monomers with targeted DP=50-150. Poly(2-hydroxypropyl methacrylamide) (PHPMA) was acting as CTA and polyethylimine (PEI) as a template with DP= 232 and 44 respectively. In that study, the templated polymerization proceeded slower than the homogeneous solution polymerization. It was hypothesized that the 17 times larger template concentration used in experiments by Ding et al., causes an increase in the viscosity, and thus, slows down the overall polymerization kinetics [13].

Simulations of RAFT polymerization accounting for the main reactions of the experimental process in synthesis of linear homopolymers were performed by Gavrilov et. al. [15] Using a Monte Carlo algorithm implemented with dissipative particle dynamics they found that if the RAFT/initiator ratio is large, a simplified model with no termination and intermediate radical formation can be used with good enough accuracy. Using the simplified model, they studied the heterogeneous polymerization -polymerization induced self-assembly (PISA). It was shown that the incompatibility between species results in different chain length distributions and polydispersity. The latter noticeably changes the phase behavior of the copolymer and the micelle size.

Due to the lack of similar experiments studying for a wide range of concentrations, molecular simulation can be used to understand the underlying effects of chemical feedback on the kinetics and the size of the resulting PCMs. To this purpose, we perform coarse-grained molecular dynamics (CGMD) simulations using a stochastic reaction model (SRM) for the polymerization of the charged copolymer block. This approach has been successfully used to study radical and living polymerizations in solution, bulk, and flat surfaces [16,17]. Also, CGMD has been widely used to study the self-assembly behavior of copolymers [18,19,20]. Here, we compute the polymerization rate constants, the local monomer concentrations, and the polydispersity of the synthesized diblock copolymer for both the templated and non-templated polymerization to explain the experimental findings. The variation of these properties with the total solution concentration, the neutral block length, the targeted polymerization length, the template length, and the excluded volume interactions between the template and the monomers were also calculated. A new clustering algorithm in python based on graph theory was developed to compute the size and shape of micelles obtained from polydisperse diblock copolymers for the templated PIESA and the non-templated polymerization followed by the assembly with the addition of the template.

2. Model

2.1. Coarse-Grained Molecular Dynamics Simulation Details

All CGMD simulations were performed using the open-source Large-scale Atomic/Molecular Massive Parallel simulator (LAMMPS [21]). Periodic boundary conditions were imposed in all directions. The Murat-Grest bead-spring model [22,23,24,25] was used to describe homopolymer chains consisting of neutral beads (A type), templates of negatively charged beads (C type), positively charged monomer beads (B type), initiators I, and counterions. The van der Waals and electrostatic interactions were modeled by the Lennard-Jones (LJ) and Coulombic potentials respectively. The bonded interactions were modeled using the finitely extensible nonlinear elastic (FENE) potential [22]. The simulations were performed using dimensionless units. A bond stiffness k of 25ε/σ² (where ε and σ are the LJ parameters, both set equal to 1) and a maximum bond extension distance r of 1.5σ were used. The solvent was treated implicitly via the Langevin thermostat [20]. The long-range electrostatic interactions between the charged beads were handled using the particle-particle particle-mesh (PPPM) method [26] with Bjerrum length l_B =1. The real-space cutoff was set to 5.0σ. Different cut-off distances in the LJ potential were used [9,20] to describe the interactions between beads. C-C and monomer-monomer interactions were considered attractive with r_cij=2.5σ. All the other interactions were considered repulsive with r_cij=2^1/6σ.

The simulations were performed at a reduced temperature T*=k_BT/ε=2 corresponding to bad solvents conditions [22] for C-C and monomer-monomer interactions. The solvent conditions for the charged moieties are determined from the balance of hydrophobic attractions between beads and the electrostatic repulsions between charges; For non-neutralized charges, the electrostatic repulsions are predominant, and therefore, charged moieties are hydrophilic. Conversely, the neutralization of charges leads to the predominance of attractions between beads, making neutralized moieties hydrophobic. All types of beads were considered to have the same mass (m=1).

In the simulations of the non-templated polymerization, mixtures containing 500 homopolymer chains and positively charged monomers and counterions were used. The homopolymers consist of 25, 50, and 100 neutral A type beads (i.e., A₂₅, A₅₀, A₁₀₀). The number of monomers and counterions is determined by the target length of the polymerized copolymer block consisting of 20 type beads (B₂₀). In all simulations, 1500 I type beads were used. In the templated polymerization simulations, negatively charged templates consisting of 20, 40, 80 and 125 C type beads (C₂₀, C₄₀, C₈₀, C₁₂₅) were added to the mixtures. These are shown as red chains in Figure 1. The ratio of charged monomers to opposite charged template beads was set equal to 1:1 in most of the simulations since it was verified [9,13] that this leads to a high number of micelles with reasonable aggregation numbers (up to N=150). Simulations with a ratio of 1:2 were also performed to compare with experimental findings [13]. All solutions are electroneutral with the addition of the appropriate number of counterions. The total concentration of beads of all types in the simulation box was varied according to [Φ]=0.04, 0.12, 0.24 and 0.36.

In all simulations, initially, 1 million timesteps were performed with an integration step of $\Delta t=0.006\tau (\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\tau =\sqrt{\frac{m{\sigma }^2}{\epsilon }})$ setting all cutoff radii equal to r_cij=2^1/6σ to eliminate any biases introduced from the initial conformation. Then, the systems were allowed to equilibrate for half a million timesteps.

2.2. Modeling the Polymerization

For the stochastic polymerization of the positive charged monomers for the synthesis of linear AB diblock copolymers shown in Figure 2, the “bond /create” functionality of LAMMPS was used. This functionality is based on a Monte Carlo algorithm that creates new bonds between atoms according to specific criteria. Possible bond pairs are identified when two non-bonded beads i and j are within a set distance R_cutoff of each other, given that the maximum number of bonds allowed per bead is not reached (i.e., 2 for linear chains). If multiple neighbors are within R_cutoff of a bead, the closest one is chosen as the sole bond partner. This bond can be created based on a predetermined reaction probability (RP). The number of maximum bonds and the types of the beads can be changed after a successful bond creation. A check for possible new bonds is performed every N_every time steps during the simulation. A schematic representation of the “bond/create” scheme is shown in Figure 2.

Three different polymerization steps were performed here: (a) the activation of the neutral homopolymer chain end bead by the addition of initiator beads (Figure 2b), (b) the bond creation between the homopolymer active center and the monomer (Figure 2c), and (c) the propagation of polymerization for the synthesis of the charged B block of the diblock copolymers. The addition of the neutral initiator alters the length of the neutral block by one bead e.g., A₅₀ to A₅₁. Since the “bond/ create” is a Monte Carlo algorithm, only one of the polymerization steps can be executed at a time. This means that only one bond can be created per timestep. R_cutoff was set equal to 1σ in all simulations. This stochastic algorithm approaches the achievements of RAFT polymerization in the sense that the main target of RAFT polymerization is to find the suitable stoichiometry (and the appropriate probabilities for the different steps in the simulation) in order the polymerization of the diblock copolymer block to take place in the vast majority from the end of the neutral block [15]. Our algorithm allows all monomers to polymerize exclusively from the end of the neutral block. As we mentioned in the introduction if the RAFT/initiator ratio is large, a simplified model with no termination and intermediate radical formation, as in our simulations, can be used with good enough accuracy [15]. Thus, our model is a coarse-grained model of the full RAFT mechanism.

2.3. Clustering Analysis

In the templated reaction assembly, after the polymerization phase is complete, the simulation was carried out for 60 million timesteps with integration step $\Delta t=0.006\tau$. The duration of the simulation was determined by the relaxation time of the tracer autocorrelation function [27,28] (Equation (S1) in the Supporting Information) of the instantaneous chains involved in a micelle. In the non-templated reaction, after the completion of the polymerization, simulations with the newly formed AB diblock copolymers along with template chains and counterions are performed to study the micellization at the desired concentrations. These simulations were performed for 15 million timesteps with $\Delta t=0.006\tau$. The properties of interest are calculated from 2000 to 4000 snapshots using the block average method with ten blocks. Following the Stillinger criterion [29] a diblock copolymer and a template chain were assumed to reside to the same micelle if any two oppositely charged beads B and C, were found within 1.5σ. In our previous study [9] on the micellization through complexation of oppositely charged diblock copolymers, we have shown that two beads of identical charge to be within 1.5σ is highly unlike. Thus, clustering algorithms that do not distinguish the charge types of the beads can safely be used for the micellization study of oppositely charged polymers.

We used graph clustering analysis to analyze the simulation data. We first identified the micelles with the data clustering algorithm DBSCAN implemented in the Python library Sklearn [30] with a maximum allowable neighborhood radius of 1.5σ. For a point (bead) to be considered as a core point, at least two points (including the point itself) must be in the neighborhood. We used a precomputed neighbor sparse array as input to the DBSCAN algorithm. To compile this array, we used the KDTree neighbor data structure from the Python library SciPy [31], and in particular the Sparse Distance Matrix algorithm with the max distance between two points of 1.5σ (note that the distance matrix algorithm ignores points with distance greater than the max distance parameter). The KDTree neighbor data structure takes periodic boundary conditions into account, and thus, the clustering analysis includes the periodic images.

To identify polymer chains, we used the Python NetworkX library [32]. Accordingly, beads were represented as nodes and bonds were represented as edges. From the graph created by this library, we could extract the polymer chains using the algorithm “Connected/Components”. This algorithm generates connected components from a graph, i.e., bead spring chains in our case. Then, the polymer chains were assigned to micelles based on the previous steps. To compute properties such as the radius of gyration of the micelles (core, corona, and total) and the shape anisotropy parameter κ² (Equation (S2) in the Supporting Information) the outbox coordinates were used. To this purpose, micelles split due to the periodic conditions (inbox coordinates) were determined and unified using the data clustering algorithm DBSCAN, KDTree neighbor data structure and Sparse Distance Matrix algorithm without periodic conditions.

3. Results and Discussion

3.1. Polymerization Rate

The kinetics of polymerization of the cationic monomers for the synthesis of AB diblock copolymers can be described by a pseudo-first-order reaction [13] according to

$$ln\left(\frac{{\left[B\right]}_0}{\left[B\right]}\right)=kt$$

(1)

where [B]₀ is the initial monomer concentration, [B] is the monomer concentration, k is the polymerization rate constant, and t is the time. To determine the effect of chemical feedback on the polymerization rate in the templated reaction when the concentration of monomers, the length of the neutral block, and other parameters vary, we compare with the respective polymerization rate of the non-templated polymerization, in which there is no assembly, and hence, no chemical feedback. To investigate the influence of concentration on chemical feedback, simulations of mixtures containing 500 neutral homopolymer chains A₅₀, 10000 positively charged B type monomers, 1500 initiator beads, and 10000 counterions were performed. For the templated polymerization, 500 additional templates C₂₀ and 10000 counterions were considered. The reaction probability was set to 0.125. Using different simulation box sizes, concentrations of [Φ]=0.04, 0.12, 0.24 and 0.36 were achieved. The target length of the polymerized positively charged block of the copolymer was set to 20 beads (A₅₁B₂₀). The kinetic plots obtained from the simulation are illustrated in Figure 3. The resulting polymerization rate constants are presented in Table S2.

As can be observed the polymerization rate in the templated reaction at [Φ]=0.04 is much faster than in the non-templated polymerization (Figure 3a,c). This is in full agreement with the experimental results of Bos et al. [13] reported for similar concentrations. The increase to [Φ]=0.12 increases both the reaction rates. However, the ratio of the rates of the two polymerization types becomes much smaller. Furter increase to [Φ]=0.24 and 0.36 (Figure 3b,c) leads to the opposite result, i.e., the non-templated polymerization rate becomes higher than the templated polymerization. This is in full agreement with the experimental results by Ding et al. [14] for concentrations of 16% and 50% w/w (coacervate in water). Since the reaction probability in the simulations is constant in both cases, the local monomer concentration before the polymerization takes place may is clear measure of the complex dependency of the polymerization rate on the concentration. To quantify the local monomer concentration in the templated polymerization, the radial distribution functions g(r_CB) of the template and the monomer beads were calculated from 200 simulation snapshots obtained by a trajectory of one million timesteps after the initial equilibration (Figure 4). Three independent simulations of the mixture with [Φ]=0.36 each one starting from a different initial configuration were used for the calculation of the standard deviation of g(r_CB). The integration of g(r_CB) up to a radius of ca. 3.5σ yields the number of B type beads that surround each C type bead. From this number the local concentration of monomers within the spherical volume with this radius can be computed. For the non-templated polymerization, the calculation of local concentration of monomers is straightforward, since the system is homogenous without template beads, and equals to the total monomer density in the simulation box.

As shown in Figure 5 and Table S2, the ratio of local monomer concentrations of the templated to the non-templated polymerizations decreases with the total concentration in a similar way as the ratio of polymerization rate constants decreases with [Φ] (Figure 3c). This shows that the local concentration may is the key parameter for understanding the variation of polymerization rate with the concentration. At [Φ]=0.04, the local concentration of the monomers around the template is two times the respective of the non-templated solution accelerating this way the templated polymerization. The difference in local monomer concentrations for the two types of polymerizations decreases at [Φ]=0.12. Nevertheless, the local monomer density in the templated polymerization is always higher, leading in a higher polymerization rate. At the most concentrated systems ([Φ]=0.24 and 0.36), the LJ interactions have a strong impact on the templated polymerization rate. The repulsion between the template, and the monomers prevents them from coming close, and thus, both the local density and the polymerization rate in the templated reaction become smaller than in the non-templated reaction.

To study the effect of LJ interactions, separate simulations with the LJ interaction parameters set to ε_BC=1, ε_BB=1.5, and ε_CC=1.5, at T*=3.0 were performed at [Φ]=0.24. In this way, the B-B and C-C type interactions remain in bad solvent conditions as previously, but the B-C interactions correspond to theta solvent for the neutral chains [9]. As shown in Figure 6, the templated polymerization rate becomes much higher than in the non-templated case since both the neutralized template and the monomers prefer to be close to each other, increasing this way, the local concentration around the template.

To study the effect of the template length on the polymerization rate, simulations of mixtures of 500 A₅₀ neutral chains and a varying number (i.e., 500, 250, 125 and 80) of template chains consisting of 20, 40, 80 and 125 C type beads, respectively, were performed at [Φ]=0.24. The target length of the copolymer chains was set to A₅₁B₂₀. B-B and C-C interactions are kept attractive, corresponding to bad solvent conditions (T*=2.0), while all other interactions are considered repulsive. As shown in Figure 7 the polymerization rate increases non-linearly with the increase of template length. In the mixtures with the longest template C₁₂₅, the polymerization rate approaches the rate of non-templated polymerization, which is in general higher for [Φ]>0.18. Long template chains may conform to an elongated shape which favors the attractions with the oppositely charged monomers. This increases the monomer local concentration and the polymerization rate.

The effects of chemical feedback on the polymerization rate when the neutral chain length varies from 25 to 50 and 100 A type beads for both templated and non-templated reactions are studied for the lowest concentration ([Φ]=0.04). In all cases, the target length of the B type block for the copolymers was fixed, i.e., A₂₆B₂₀, A₅₁B₂₀ and A₁₀₁B₂₀.

As shown in Figure 8, both the templated and the non-templated polymerization, the increase in the length of the neutral block leads to the linear decrease of the polymerization rate. This is because the excluded volume interactions between the A and B hinder monomers from approaching the active B type end beads (reduced monomer concentration around the active centers). The decrease in the non-templated polymerization rate is enhanced compared to the templated polymerization. As discussed earlier, in the templated polymerization a lot of monomers are stuck on the template, even before the polymerization starts. Thus, the excluded volume interactions with the A type blocks concern a smaller number of free monomers. In contrast, in homogeneous non-templated polymerization, all monomers experience excluded volume interactions with A type blocks. This hinders monomers from approaching the active end beads.

3.2. Molecular Weights and Polydispersity

The mass distributions of the B block of the A₅₁B₂₀ copolymers obtained from both types of polymerizations are presented in Figure 9 for [Φ]=0.04, 0.12, and 0.24. A C₂₀ template is used and RP is set to 0.125. From these distributions, the number (Mn) and weight (Mw) average molecular weights, and the polydispersity index (PDI=Mw/Mn) can be calculated. The results are listed in Table S1 in the Supporting Information. From the simulations, we have verified that copolymer chains consisting of one or two B type beads (i.e., A₅₁B₁, A₅₁B₂) do not participate in the micellization. Thus, these copolymers are considered as impurities and are not counted in the calculation of Mw, Mn and PDI. As can be seen in Table S1, the Mn, Mw and PDI of the B type block obtained from the templated polymerization are higher than the respective quantities in the non-templated polymerization, especially at low and moderate concentrations [Φ]=0.04 and 0.12. This is because in the templated polymerization where half of the monomers are lying on the oppositely charged templates, many neutral chains remain without or with only 1 or 2 B type beads. In the non-templated polymerization, where the monomers are homogeneously distributed in the solution, almost all A type chains participate in the polymerization. In this case, the chains have narrow molecular weights and low PDI. Experimental results of the PDI values of the whole diblock copolymer chains (PDI_diblock) and the pre-synthesized A type precursors (PDI_B) are reported [14]. To extract PDI_B and to compare with the simulation results, the following equation is used [33]

$${PDI}_{diblock}=w_A^2\left({PDI}_A-1\right)+w_B^2\left({PDI}_B-1\right)+1$$

(2)

where w_A and w_B are the ratios of the number of beads of the A and B block to the total number of beads in the diblock copolymer, respectively. Figure 10a and 10d show the PDI_B and Mw as a function of [Φ], respectively. As can be observed, in the templated polymerization, the increase in the concentration decreases both the PDI_B and Mw. In the non-templated homogeneous polymerization, the trend is the opposite; Both PDI_B and Mw increases with [Φ]. At [Φ]=0.24, the difference between the PDI_B values for the two types of polymerizations become very small. Using the experimental values ( ${PDI}_{diblock}$=1.10 and ${PDI}_A$=1.2) of Ding et al. [14], which are the same for both types of polymerizations at [Φ]=0.5 in Equation (2), we predicted the ${PDI}_B$=1.2. This value is close to the simulation mean value of 1.26 at [Φ]=0.24.

Our simulation results show that the difference in the PDI_B values between templated and non-templated polymerization at [Φ]=0.04 is significant in full agrement with the simulation results of Gavrilov et. al. [15] for PISA polymerization. They found that the polydisperity of the diblock copolymer varies from 1.07 to 2.15 as the strength of the interactions between moieties is changed. However, no direct comparison with the experimental results can be done, since Bos et al. [13] do not report the PDI_B values for the non-templated polymerization. PDI and Mw as functions of the A type block length are presented in Figure 10b, and e respectively, for [Φ]=0.04. It can be seen that the increase from A₂₆ to A₅₁ leads to a decrease in PDI. This is because the longer neutral chain covers the active end bead of the polymerized B block preventing the monomers to approch. The effect of RP on PDI and Mw are presented in Figure 10c,f for [Φ]=0.04. In the non-templated polymerization, the increase of the RP from 0.125 to 0.25 leads to the increase of Mw and PDI values. However, further increase to 0.5 has no extra effect on them because the local concentration of monomers around the active polymerization center cannot further increase.

3.3. Micelle Size and Shape

The effects of the chemical feedback on the size and shape of the micelles were studied at the lowest concentrations, i.e., [Φ]=0.04. The simulations of the PIESA one-step micellization were performed for the following mixtures: (a) A₂₅ + C₂₀, (b) A₅₀ + C₂₀, and (c) A₁₀₀ + C₂₀. The target diblock copolymers were A₂₆B₂₀, A₅₁B₂₀, and A₁₀₁B₂₀, respectively. To model the two-step micellization, after the end of polymerization, the C₂₀ templates were added to the simulation box. The simulation scheme is described in detail in the model section. The mass distribution functions of the micelles computed from the molecular dynamics trajectories using the new python code are shown in Figure 11 as a function of the aggregation number N.

Figure 11 clearly shows that in the templated polymerization micelles with higher aggregatiom numbers are formed. Regardless of the micellization method, the increase of the neutral block in the diblock copolymers from A₂₆B₂₀ to A₅₁B₂₀ and further to A₁₀₁B₂₀ leads to smaller aggregates. This is expected because the increase in the hydrophilic block makes the corona of the micelle bulkier, better protecing the hydrophobic core formed by the complexation of the oppositely charged B and C beads. However, the mass distibution profiles obtained from the PIESA differ significantly from the two-step micellization. In PIESA the mass distribution profile is a Gaussian- like function (micelles with preferential aggregatiom number are formed as shown in Figure 11b). The deviations from the Gaussian function Figure 11a,c can be attributted to the difficulties in the equilibration arising from the smaller neutral block and the higher Mw of PIESA chains respectively. In sharp contrast, the mass distribution profiles obtained from the two-step micellization do not reveal preferential aggregation [9] (decaying function with N, Figure 11d–f). The evolution of micelle mass distribution with the time in the templated polymerimerization is presented in Figure S2 for the mixture of A₅₁B₂₀ + C₂₀. Initially small micelles are formed then progresivelly the aggregation number of micelles is increasing and after completion of the polymerization the size reareangements leads in the gaussian type distribution presented in Figure 11b. Snapshots of the simulation box are presented in Figure S3 for the same mixture at concentrations [Φ]= 0.04 and 0.36 and different simulation times τ.

The mean squared radii of gyration <S²>_PIESA and <S²>_two-step, describing the size of the micelles for the PIESA and the two-step micellization, respectively, are shown in Figure 12 as a function of the aggregation number. <S²>_two-step is always higher than <S²>_PIESA. This finding is in full

aggrement with the experimental results of Bos et. al. [13]. The deviation between <S²>_two-step and <S²>_PIESA decreases as the length of the neutral block forming the corona becomes much larger than the B type block (Figure 12b,c). The lower values of <S²>_PIESA may be due to the higher PDI_B which is 1.3 for PIESA, compared to 1.1 for the two-step micellization. Van der Kooij et. al. [3] and Gavrilov et. al. [15] have shown that the diblock copolymer chains with higher PDI resulted in denser packing in the micelle core. Thus, the overall mean squared radius of gyration value of the polyelectrolyte complex micelles was much lower than the respective of micelles formed by the low PDI copolymers. The shape anisotropy parameter [9,20] κ² (Equation (S2) in the Supporting Information) is shown in Figure 13 as a function of the aggregation number of the micelles for the two micellization schemes. From Figure 13, it is evident that micelles with very small aggregatiom numbers (N<10) are elongated (κ² > 0.1). The micelles are shperical (κ² < 0.1) at moderate aggregation numbers, and again elongated for higher aggreegation numbers (N>90). The κ² values for micelles with high aggregation numbers are scattered since this calculation suffer from bad statistics because such big micelles rarely form in the simulation. In general, micelles with aggregation numbers 20< N <80 formed by the PIESA micellization are more spherical than the ones in the two-step micellization for the A₂₆B₂₀ + C₂₀, and A₅₁B₂₀ + C₂₀ mixtures. In contrast, the micelles obtained from the A₁₀₁B₂₀ + C₂₀ mixture (i.e., the system with the longest neutral A type block) have similar shape in both micellization schemes. The reason is that the shape of the large corona consisting of the A₁₀₁ blocks predominaly determines the overall shape of the micelles.

So far in this study, we have focused on micelles formed by mixtures in which the ratio of the charged template beads to the oppositely carged monomers is 1:1. To study the effect of chemical feedback in mixtures with an excess of templated negative beads, simulations are performed for a 2:1 ratio. The mass distributions of the micelles obtained from A₁₀₁B₂₀ + C₄₀ mixtures are presented in Figure 14 as a function of the aggregation number for both PIESA and two-step schemes. Our results show that, regardless of the micellization method, only very small aggregates are formed (N ≤ 7). This finding is in full agreement with the experimental results of Boss et al. [13], and also agrees with previous theoretical predictions of PCMs being formed only at approximately equal charge stoichiometries [9].

4. Conclusions

To understand the effects of chemical feedback in the templated polymerization of a charged copolymer block and the simultaneous self-assembly with the opposite charged template on the kinetics and the size of the resulting complex micelles we performed CGMD simulations. The polymerization was based on a Monte Carlo stochastic reaction model. We show that chemical feedback fundamentally changes both the polymerization and self-assembly. At low solution concentrations [Φ], the location of monomers on the templates significantly accelerates the polymerization, while at high [Φ], the local monomer concentration on the templates becomes lower than that of the solution, thus, slowing down the polymerization. This is in full agreement with the experimental data. Both the increase in the template length and the decrease in the neutral block length led to the linear increase of the polymerization rate.

In the templated polymerization, the Mn, Mw, and PDI values of the B type block were computed to be higher than the respective quantities in the non-templated polymerization, especially at the low and moderate concentrations [Φ]=0.04 and 0.12. The micelles formed by the templated PIESA method are have higher aggregatiom numbers than those formed by the two-step micellization. However, the <S²>_two-step is always higher than <S²>_PIESA, in aggrement with the experimental results. This is most probably a result of the higher PDI of the templated polymerization diblock copolymer block and the difference in the shape of micelles. The micelles with moderate aggregation numbers formed by the PIESA method are more spherical than the ones formed in the two-step micellization in full agreement with experimental findings. This useful insight into the templated reaction assembly process in polymers obtained from molecular simulation is necessary for the rational design of new synthetic supramolecular materials.