Multicore Assemblies from Three-Component Linear Homo-Copolymer Systems: A Coarse-Grained Modeling Study

Multicore polymer micelles and aggregates are assemblies that contain several cores. The dual-length-scale compartmentalized solvophobic–solvophilic molecular environment makes them useful for, e.g., advanced drug delivery, high-precision synthesis platforms, confined catalysis, and sensor device applications. However, designing and regulating polymer systems that self-assemble to such morphologies remains a challenge. Using dissipative particle dynamics (DPD) simulations, we demonstrate how simple, three-component linear polymer systems consisting of free solvophilic and solvophobic homopolymers, and di-block copolymers, can self-assemble in solution to form well-defined multicore assemblies. We examine the polymer property range over which multicore assemblies can be expected and how the assemblies can be tuned both in terms of their morphology and structure. For a fixed degree of polymerization, a certain level of hydrophobicity is required for the solvophobic component to lead to formation of multicore assemblies. Additionally, the transition from single-core to multicore requires a relatively high solvophobicity difference between the solvophilic and solvophobic polymer components. Furthermore, if the solvophilic polymer is replaced by a solvophobic species, well-defined multicore–multicompartment aggregates can be obtained. The findings provide guidelines for multicore assemblies’ formation from simple three-component systems and how to control polymer particle morphology and structure.


Introduction
Multicore polymer assemblies, i.e., polymer aggregates consisting of multiple microphasesegregated cores and in the case of micelles, a solvophilic corona, are functional selfassembling materials with important applications in, e.g., water purification, catalysis, pharmacology, electronics, and oil recovery [1][2][3]. The high tunability of the chemistry, size, and distribution of the solvophobic cores of the multicore assemblies make them unique nanoscale capsules that can be loaded with chemical reagents or catalysts, making them nanoreactors [4]. Such polymeric capsule nanoreactors are ubiquitous, for example, in emulsion polymerization [4,5] and confined catalysis [6,7], but also as artificial organelles in synthetic biology and therapeutic applications [8][9][10].
Furthermore, multicore-multicompartment assemblies and micelles, which differ from the basic multicore micelles by having multiple levels of compartmentalization, are linear ABC tri-block ter-polymers containing a solvophilic midblock and two mutually incompatible solvophobic head groups. In another study, Wengenmayr et al. [51], via a mean-field model, observed the splitting of unimolecular single-core micelles of dendriticlinear copolymers to a multicore structure, with increasing dendrimer generation and decreasing solvent selectivity. On the computational side, DPD simulations which employ particle-based, mesoscale level modeling of the polymer chains and their assembly and dynamics, have turned out to be an adaptable and reasonably successful tool for examining the formation and assembly responses of polymeric micelles [46,[52][53][54]. In particular, multicore and multicompartment micelles have been studied via DPD by Zhao et al. [3], who investigated the formation of multicompartment micelles from a dilute aqueous solution of the mixture of di-block copolymer poly (ethyl ethylene)-block-poly (ethylene oxide) (PEE-b-PEO) and homopolymer poly(propylene oxide) (PPO). Chen et al. [54] also employed DPD simulations to examine the multicomponent-multicore micelle formation from two kinds of star-shaped copolymers and the degradation of multicomponent-multicore micelles. Chen et al. [55], via DPD, examined the self-assembly of π-shaped copolymers that had rigid backbones and two flexible side chains in a selective solvent. By varying the model parameters, such as the concentration of the copolymer, the nature of the solvent, and the repulsive interactions between the two coiled side chains and the rigid backbone block, they observed a variety of nanostructures, including interacting micelles, multicore micelles, nanocages, and nano-tires. Callaway et al. [2] studied the self-assembly of a lipophilic−hydrophilic−fluorophilic triblock copolymer system that formed non-spheroidal micelles with multicores of both lipophilic and fluorophilic species.
To find a more simple and economical approach of forming multicore assemblies, we present here a study where an amphiphilic di-block copolymer and two linear homopolymers, one weakly and the other strongly solvophobic, are mixed in a solvent. Since the two homopolymers have a significant difference in their level of solvophobicity and the weakly solvophobic species is relatively compatible with the solvent, we denote the weakly and strongly solvophobic homopolymers as solvophilic (solvophile) and solvophobic (solvophobe), respectively. Such a setup is easily attainable in industrial scale application processes and provides many readily accessible tuning parameters, for example, choice of the polymers, their concentration, degrees of polymerization, and system composition. To reduce the number of variables in the system, the di-block components are considered identical in chemistry to the free solvophile and solvophobe. Here, we show the polymer property range that leads to multicore assemblies, via DPD simulations of this model system. We map the aggregation dependencies and report the range of polymer system parameters over which multicore assemblies form for a fixed degree of polymerization. Furthermore, well-defined multicompartment aggregates can form if the solvophilic polymer is replaced by a more solvophobic species.

Simulation Details
Polymer particle formation and aggregation was examined via DPD [56]. It is a mesoscale simulation method in which a group of atoms comprise beads that capture the collective interactions of the encompassed atoms and groups. In the DPD approach, the total force exerted on the ith bead by all the other beads in the system is provided as: Here, the four force components correspond to a conservative force distribution, and a harmonic spring force → F ij s = C → r ij acts between the polymer beads. The spring force is specific to the DPD setup for a polymer system and exists only for neighboring beads in the chain.
In Equation (1), the total force is summed over all particles i = j, within a cut-off radius R c . Here, ν i are the velocities of beads i and j, respectively. The constant a ij describes the maximum repulsion between interacting beads. γ and σ are the friction coefficient and noise amplitude, respectively. ξ ij is a random number with zero mean and unit variance, chosen independently for each interacting pair of beads at each time step ∆t, while ω C , ω D , and ω R are three weight functions for conservative, dissipative, and random forces, respectively. For the conservative force, we chose ω C r ij = 1 − r ij /R c for r ij < R c and ω C r ij = 0 for r ij ≥ R c . The choice of ω D r ij and ω R r ij is connected via the relation ω D = ω R r ij 2 , σ 2 = 2γk B T, which comes from the fluctuation-dissipation theorem. Here, we chose a simple form of ω D and ω R following Groot and Warren [57]: The equations of motion in the DPD simulations are integrated using a modified version of the velocity-Verlet algorithm. For simplicity, the DPD simulations employ reduced units such that the cut-off radius R c , the bead mass m, and k B T specify the unit magnitude for the distance, mass, and energy units, i.e., R c = m = k B T = 1.
For the like-like-repulsion parameters, a ii , in Equation (1), Groot and Warren derived an expression a ii ≡ a jj = 75 ρ , based on the equation of state for the soft repulsive DPD fluid and the compressibility of bulk water. Here, ρ is the number density in the reduced units. For p = 3, a ii = 25, which corresponds to the compressibility of water. They also suggested the following linear relation between a ij and the Flory-Huggins χ parameter at p = 3: The harmonic spring force constant C was set to be 4.0 [3] to keep the polymer beads connected.
The bead interaction parameters, a ij , are presented in Table 1. Here, we take the interaction parameters between the same types of A, S, and W beads to be a WW = a SS = a AA = 25, reflecting the compressibility of water. This value also describes unbiased miscibility between the components, i.e., χ = 0 in the Flory-Huggins model. We set the interaction parameter between the solvophilic block S and the solvent W to a SW = 30 (χ = 1.53), as this resides between a ij = 26.64 that emulates the θ-state (χ = 0.5) and the critical value for hydrophobicity a ij = 32 [58]. The value of the interaction parameter between A and S is selected as a AS = 72 (χ = 14.37). This implies that A and S components are mutually incompatible. The interaction parameter between A and solvent bead W is selected as a AW = 115. This value models extremely unfavorable interactions corresponding to the Flory-Huggins parameter χ = 27.72. Throughout, we consider a system with 20% of polymeric phase concentration. The basic system is composed of A 19 (4%), A 1 -b-S 6 (13%), and S 6 (3%) in solvent medium. As our target is to provide a general perspective to the aggregation behavior of the system chosen, the parameters employed do not rise from a bottom-up coarse-graining of a Polymers 2021, 13, 2193 5 of 17 particular chemical system. However, the bead-spring chain models were chosen to be in the range of oligomers; generally, in a bottom-up coarse-graining study with DPD methodology, each bead-spring model chain can be mapped on a real polymer chain with a degree of polymerization in a range of oligomers to polymers, depending on the degree of coarse-graining.
The DPD simulations were carried out in a fully periodic cubic simulation box of size 40R c × 40R c × 40R c . Assessment of the simulation box size was performed by reproducing the multicore morphologies in a significantly larger system of 100R c × 100R c × 100R c , see Supplementary Materials for details. All simulations were carried out using the LAMMPS (large-scale atomic-molecular massively parallel simulator) package [59]. The DPD simulations were performed in the NVT ensemble (k B T = 1). The equations of motion were integrated using the modified velocity-Verlet algorithm, with λ = 0.65 and a time step ∆t = 0.05τ. The total simulation times varied between 5 × 10 5 and 5 × 10 6 steps.

Characterization of Self-Assembled Structures
Classification to aggregates was performed based on a separation distance criterion [60,61]: any bead, or set of beads, within a cut-off distance 1.5R c from each other are part of the same aggregate. The cut-off distance was determined via visual analysis. The classification scheme is not particularly sensitive to the precise value of the cut-off distance, and the same value was used for all systems.
The aggregate sizes were characterized by calculating the average particle number (X n ) and the weight average (X w ), as: where N i is the number of aggregates containing X i beads. Equilibration of the simulation system was assessed via monitoring the aggregate size distributions, and equilibrium state was deemed when the average aggregate size stabilizes at a constant level. Polydispersity was evaluated via the polydispersity index: Notably, P = 1 corresponds to a monodisperse system with identical core sizes. The larger the deviation from unity, the more polydisperse the cores in the assemblies are.
The internal structure of the aggregates and phase separation in them was quantified via defining an average order parameter, ϕ, as: where the summation is over N slices of equal thickness along the xy-plane, and in each slice, the volume fractions of beads A and B are ϕ Ai and ϕ Bi , respectively. Shape of the aggregates was assessed by determining their deviation from the spherical shape via asphericity, δ, calculated as: where λ 1 2 , λ 2 2 , and λ 3 2 are the principal moments of the gyration tensor and λ l > λ 2 > λ 3 . The value of δ is positive, except for a perfect sphere, where it is precisely zero.

Results
The DPD simulations comprise a mixture of a linear solvophobic polymer, a linear solvophilic polymer, and a linear amphiphilic di-block copolymer in solvent medium, referred to as A x , S y , and A c -b-S d , respectively. The model of the system components is shown in Figure 1. The subscripts in each abbreviation refer to the chain length. In all simulation runs, the initially randomly dispersed solvophobic-free chains self-assembled into small micelles very quickly. After this, small micelles gradually formed, and coalesced to form large multicore assemblies. Figure 2 shows the time evolution of the aggregation process corresponding to this self-assembly for a 20% polymer concentration of A 19 (4%), A 1 -b-S 6 (13%), and S 6 (3%) in aqueous solutions by simulation snapshots, and the corresponding average sizes of aggregates measured by the average bead number and its weighted average in the aggregates of the system. As an initial step, we extensively tuned the polymer system composition. The final self-assembling structure is highly dependent on the composition of the polymer mixture. This composition was selected as the focus of the work since multicore assemblies readily and stably assembled in this composition and other close-by compositions tested. erage in the aggregates of the system. As an initial step, we extensively tuned the polymer system composition. The final self-assembling structure is highly dependent on the composition of the polymer mixture. This composition was selected as the focus of the work since multicore assemblies readily and stably assembled in this composition and other close-by compositions tested.
The initially formed small micelles are single-core aggregates with double-core corona structure. They are composed of a solvophobic spherical core of A19, a layer of solvophobic block of the copolymer as the first corona shell, and both a solvophilic block of copolymer and S6 as the second shell. Finally, the small micelles merge with the neighboring micelles to a large spherical aggregate. As shown in Figure 2a-d, the small micelles agglomerate, which is followed by partial merging of their corona regions. However, the small hydrophobic cores remain as separate hydrophobic compartments in the final aggregate. Figure 2e shows a zoomed-in cross-section of the final multicore configuration. Additional time evolution snapshots and a visualization of the aggregate internal structure are provided in the Supplementary Materials Figure S1. Moreover, comparison with the assembly response in the system with the box size of 100Rc × 100Rc × 100Rc (see the Supplementary Materials Figure S2) shows that the simulation box size influences only the aggregate size but not the multicore structure.
Due to the simulation box size dependence of the particle-like assembly size in the simulations, we note that a more general interpretation of the multicore assembly response is a concentrated phase with phase-separation into multiple droplets. In the case of a finite sized aggregate, this would refer to the formation of a concentrated phase droplet containing several small droplets that differ in their polymer composition from the surroundings. The initially formed small micelles are single-core aggregates with double-core corona structure. They are composed of a solvophobic spherical core of A 19 , a layer of solvophobic block of the copolymer as the first corona shell, and both a solvophilic block of copolymer and S 6 as the second shell. Finally, the small micelles merge with the neighboring micelles to a large spherical aggregate. As shown in Figure 2a-d, the small micelles agglomerate, which is followed by partial merging of their corona regions. However, the small hydrophobic cores remain as separate hydrophobic compartments in the final aggregate. Figure 2e shows a zoomed-in cross-section of the final multicore configuration. Additional time evolution snapshots and a visualization of the aggregate internal structure are provided in the Supplementary Materials Figure S1. Moreover, comparison with the assembly response in the system with the box size of 100R c × 100R c × 100R c (see the Supplementary Materials Figure S2) shows that the simulation box size influences only the aggregate size but not the multicore structure. Schematic illustration of the solvophobic and solvophilic homopolymer, block copolymer, and solvent. Different color codes for beads A and S in the copolymer (white and green) and homopolymer (red and purple) chains are to clarify the chain locations in the assemblies in the simulations. Despite the color coding, beads A and S are chemically identical in both polymers. The time evolution of the aggregate sizes in terms of and for both the solvophobic cores and total aggregates in Figure 2f shows three stages of evolution. In the first stage (I), small single-core micelles, with the cores and double-corona composed A19 and A1-b-S6 respectively, are formed, in stage (II), the small micelles aggregate to form the multicore aggregate, and stage (III) corresponds to equilibrium. In stage (I), and increase for both the A19 solvophobic cores and the whole aggregate, and the nucleation of several small aggregates takes place. At the end of stage (I), the small micelles reach their approximate equilibrium morphology and size, but the total aggregate is still evolving both in shape and size. In the second stage (II), and for the solvophobic cores saturate, but the whole aggregate still increases in size. This stage occurs over a longer time period than the initial micelle formation as the system changes slowly by the small micelles coalescing their solvophilic outer coronas. Such events show a sharp increase in and of the whole aggregate in Figure 2. This process may also have started in stage (I). The final stage (III) is equilibrium, in which and for the whole aggregate saturate, too. The different stages are visualized by the snapshots in Figure 2a The different formation stages were also assessed via calculating the order parameter both for the solvophobic chains' phase separation to form the cores and for the polymer Due to the simulation box size dependence of the particle-like assembly size in the simulations, we note that a more general interpretation of the multicore assembly response is a concentrated phase with phase-separation into multiple droplets. In the case of a finite sized aggregate, this would refer to the formation of a concentrated phase droplet containing several small droplets that differ in their polymer composition from the surroundings.
The time evolution of the aggregate sizes in terms of X n and X w for both the solvophobic cores and total aggregates in Figure 2f shows three stages of evolution. In the first stage (I), small single-core micelles, with the cores and double-corona composed A 19 and A 1 -b-S 6 respectively, are formed, in stage (II), the small micelles aggregate to form the multicore aggregate, and stage (III) corresponds to equilibrium. In stage (I), X n and X w increase for both the A 19 solvophobic cores and the whole aggregate, and the nucleation of several small aggregates takes place. At the end of stage (I), the small micelles reach their approximate equilibrium morphology and size, but the total aggregate is still evolving both in shape and size. In the second stage (II), X n and X w for the solvophobic cores saturate, but the whole aggregate still increases in size. This stage occurs over a longer time period than the initial micelle formation as the system changes slowly by the small micelles coalescing their solvophilic outer coronas. Such events show a sharp increase in X n and X w of the whole aggregate in Figure 2. This process may also have started in stage (I). The final stage (III) is equilibrium, in which X n and X w for the whole aggregate saturate, too. The different stages are visualized by the snapshots in Figure 2a-d.
The different formation stages were also assessed via calculating the order parameter both for the solvophobic chains' phase separation to form the cores and for the polymer components (solvophobic, di-block copolymer, and solvophilic chains) to form overall aggregates. A larger value of order parameter indicates a higher degree of order, or in this case, component segregation [62]. For the overall aggregate, order parameter values of 0.514 ± 0.001 and 0.134 ± 0.001 were obtained for the final aggregates and the solvophobic cores, respectively. The time evolution of the order parameter is presented in the Supplementary Materials Figure S3. The data show the same response as demonstrated via the X n and X w datasets in Figure 2, including the different stages of aggregation. To examine the dispersity of the cores and aggregation, we also calculated the polydispersity, P, from Equation (5). The polydispersity data presented in Supplementary Material Figure S4 show that the initially forming polydisperse cores become more monodisperse in size during stage (I) and early in stage (II), i.e., the hydrophobic cores still evolve in stage (II). The practical significance of this is that controlling the growth rate in stage (II) provides a key handle for controlling the size polydispersity of the small cores in applications. For this purpose, e.g., the copolymer block lengths, the chemical difference between the blocks and copolymer concentration provide easily accessible means.
To assess the shapes of the aggregates and their fluctuations, the time evolution of the principal components of the gyration tensor (λ l , λ 2 , and λ 3 ) for the largest aggregate in the system are presented in Figure 3. At the initial stages of the aggregation, as shown by the inset of Figure 3a, λ l , λ 2 , and λ 3 are almost equal, indicating the formation of small, separated spherical aggregates (Figure 3b, aggregate A). These small micelles rapidly find each other and form small aggregates with multiple cores (aggregate B). The aggregation continues by formation of anisotropic supramolecular structures and a transition from spherical particle to worm-like multicore assemblies, as shown by λ l increasing, while λ 2 and λ 3 remain nearly constant, see aggregates C-E in Figure 3b. The shapes of these worm-like aggregates fluctuate, see aggregates D, E, and finally, the system approaches (aggregate G) and reaches (aggregate H) a final spherical structure. Deviation of the principal component values here corresponds to structural fluctuations, and the overall structure is a sphere. Figure 3a also presents an inset showing the corresponding asphericity, δ. In the beginning of the assembly, δ is close to unity, which matches the corresponding morphologies in Figure 3b. The final equilibrated aggregate is spherical, as shown by δ smaller than 0.01.
To study the range of conditions where the observed multicore micellization occurs and to map out its polymer species' sensitivity, we varied the solvent selectivity for the solvophobic and solvophilic homopolymers. This was performed by simply changing the corresponding bead-solvent interaction parameters between 35 and 115 for a AW and between 30 and 115 for a SW, respectively. We first discuss the results shown in Figure 4 when the solvophobic homopolymer solvent selectivity varies by reducing the interaction between beads A and W, i.e., a AW . Decreasing a AW from the reference value of 115 (making the bead type A less solvophobic) decreases the number of solvophobic cores. This is a direct consequence of the solvophobic polymer becoming less solvophobic or more soluble. Multicore micellization still persists for a AW , that is about 65. For a AW less than that, there is a gradual transition to a regime where single-core aggregates surrounded by the solvophilic polymer form. These singlecore assemblies are structurally similar to the small micelles formed in multicore aggregate formation step I, see Figure 4, visualization of a system corresponding to a AW = 55. On the other hand, when the homopolymers become more and more soluble, the system no longer forms clear aggregates with internal structure. Instead, the polymer chains form partially soluble spherical particles in the solvent due to their attractive interactions, see Figure 4, visualization of a system corresponding to a AW = 35. The mean aggregation number, , and polydispersity, P, of the solvophobic cores versus the more solvophilic polymer-solvent interaction parameter aAW are presented in Figure 4. The data show the general trend of a decreasing aAW, leading to a significant increase in . There are two sharp steps in the aggregation number vs. aAW. The increase taking place with aAW between 75 and 65 corresponds to a transition from multicore to single-core structures. The second transition from single-core aggregate to soluble particle takes place with aAW decreasing from 45 to 35. It occurs because in this range of aAW, the solvophobic A19 chains are relatively solvophilic, but still assemble as spherical particles that have an interface with the solvent. This results from the stronger like-like interaction between the A beads (aAA = 25) than between the A-W beads (aAW = 35). Block copolymers are mostly solvated, with only a few copolymer chains remaining connected to the particle surfaces by their A block head, see Figure 4, panel corresponding to aAW = 35 for visualization. The mean aggregation number, X n , and polydispersity, P, of the solvophobic cores versus the more solvophilic polymer-solvent interaction parameter a AW are presented in Figure 4. The data show the general trend of a decreasing a AW , leading to a significant increase in X n . There are two sharp steps in the aggregation number vs. a AW . The increase taking place with a AW between 75 and 65 corresponds to a transition from multicore to single-core structures. The second transition from single-core aggregate to soluble particle takes place with a AW decreasing from 45 to 35. It occurs because in this range of a AW , the solvophobic A 19 chains are relatively solvophilic, but still assemble as spherical particles that have an interface with the solvent. This results from the stronger like-like interaction between the A beads (a AA = 25) than between the A-W beads (a AW = 35). Block copolymers are mostly solvated, with only a few copolymer chains remaining connected to the particle surfaces by their A block head, see Figure 4, panel corresponding to a AW = 35 for visualization. s 2021, 13, x FOR PEER REVIEW 10 of 18 The color scheme is the same as that of Figure 1. The data graphs show mean aggregation number , polydispersity index (P) of the solvophobic cores, and the order parameter (φ) of the overall aggregate versus the more solvophobic polymer-solvent interaction parameter aAW. Error bars are smaller than the symbol size. The scatter in the data is due to finite size effects in the simulations.
We note that the polydispersity index of the aggregates is relatively insensitive to aAW. Nevertheless, the snapshots in Figure 4 indicate that decreasing aAW eventually leads to less polydisperse cores when the multicore structures are replaced by single-core aggregates. Furthermore, the variation of aAW, i.e., the solvophobicity of the solvophobic polymer, provides a means to tune the solvophobic core size for single-core aggregates.
Additionally, the dependence of the aggregate order parameter on aAW is shown in Figure 4. The data show a decreasing trend with decreasing aAW. This is a direct consequence of the increasing compatibility between the solvophobic and solvent beads, which leads to a smaller degree of molecular separation.
The second important parameter in the system is the solvophilic homopolymer solvophilicity parameter aSW, whose reference value was set to 30. Figure 5 shows simulation snapshots with a larger value of aSW = 65. Increasing the solvophobicity of the S6 homopolymer leads to the chains entering the aggregate to reduce their solvent contact, which promotes the formation of multicore-multicompartment assemblies. Similar to the case where aSW = 30, the more solvophobic A19 polymer aggregates to form cores, but now also the less solvophobic polymer S6 is sufficiently solvophobic to aggregate and form cores instead of just filling up the aggregate. Both types of cores are covered by the A1-b-S6 diblock copolymer to form small micelles. Figure 6 shows the effect of varying aSW on the aggregate structure via snapshots and analysis of the aggregate structural changes as the function of aSW. The snapshots show that increasing aSW promotes the aggregation of the solvophobic irregular cores to the center of the total aggregate and formation of a more united core. It also leads to the initial micelles having more elongated core structures. This elongation is promoted by the increasing surface contact area between the cores and the solvophilic blocks. Similar to the multicore aggregate formation, small micelles assemble together, but as there are now two types of cores, the aggregate structure is multicore-multicompartment. Notably, the larger aSW values lead to the solvophilic S6 chains being present also in the center of the aggregate, surrounded by the solvophobic cores. Simulation snapshots and analysis of the aggregate structural changes when a AW is varied in the system containing A 19 (4%), A 1 -b-S 6 (13%), and S 6 (3%). Solvent beads are omitted for clarity. The color scheme is the same as that of Figure 1. The data graphs show mean aggregation number X n , polydispersity index (P) of the solvophobic cores, and the order parameter (ϕ) of the overall aggregate versus the more solvophobic polymer-solvent interaction parameter a AW . Error bars are smaller than the symbol size. The scatter in the data is due to finite size effects in the simulations.
We note that the polydispersity index of the aggregates is relatively insensitive to a AW . Nevertheless, the snapshots in Figure 4 indicate that decreasing a AW eventually leads to less polydisperse cores when the multicore structures are replaced by single-core aggregates. Furthermore, the variation of a AW , i.e., the solvophobicity of the solvophobic polymer, provides a means to tune the solvophobic core size for single-core aggregates.
Additionally, the dependence of the aggregate order parameter on a AW is shown in Figure 4. The data show a decreasing trend with decreasing a AW . This is a direct consequence of the increasing compatibility between the solvophobic and solvent beads, which leads to a smaller degree of molecular separation.
The second important parameter in the system is the solvophilic homopolymer solvophilicity parameter a SW , whose reference value was set to 30. Figure 5 shows simulation snapshots with a larger value of a SW = 65. Increasing the solvophobicity of the S 6 homopolymer leads to the chains entering the aggregate to reduce their solvent contact, which promotes the formation of multicore-multicompartment assemblies. Similar to the case where a SW = 30, the more solvophobic A 19 polymer aggregates to form cores, but now also the less solvophobic polymer S 6 is sufficiently solvophobic to aggregate and form cores instead of just filling up the aggregate. Both types of cores are covered by the A 1 -b-S 6 di-block copolymer to form small micelles.   Figure 6 shows the effect of varying a SW on the aggregate structure via snapshots and analysis of the aggregate structural changes as the function of a SW . The snapshots show that increasing a SW promotes the aggregation of the solvophobic irregular cores to the center of the total aggregate and formation of a more united core. It also leads to the initial micelles having more elongated core structures. This elongation is promoted by the increasing surface contact area between the cores and the solvophilic blocks. Similar to the multicore aggregate formation, small micelles assemble together, but as there are now two types of cores, the aggregate structure is multicore-multicompartment. Notably, the larger a SW values lead to the solvophilic S 6 chains being present also in the center of the aggregate, surrounded by the solvophobic cores.   Figure 1. At the bottom, the mean aggregation number, , polydispersity index (P) of solvophobic cores, and order parameter (φ) of the aggregates are shown as a function of the interaction parameter between the more solvophilic polymer and solvent aSW. Error bars are smaller than the symbol size. The scatter in the data is due to finite size effects in the simulations. Furthermore, Figure 6 presents and polydispersity of the solvophobic cores vs. aSW. The data show that the aggregation number does not have a simple relation to increasing aSW. While relatively large aggregates form for small aSW, a modest increase leads to a significant decrease in aggregate size, and further increasing aSW shows that the aggregate size increases. This can be understood if one considers that a small value of aSW = 30, i.e., good solvent solubility of the more solvophilic S6 polymer, leads to a minority of them assembling in the aggregate. Increasing aSW makes the solvophilic polymer less soluble in the solvent and thus more of it assembles in the center of the aggregate. This may lead to promoting the merging of the solvophobic cores. This is shown as the increase in in Figure 6. The polydispersity index of the aggregates in Figure 6 remains relatively insensitive to variations of aSW. For the larger aSW values, the increase in the polydispersity index results from the solvophilic homopolymers that remain incompatible with the solvophobic homopolymer aggregating into the aggregate. The presence of the solvophilic polymer in Figure 6. Aggregate snapshots and analysis of the aggregate structural changes when a SW is varied. At the top, the cross-sections of the final morphology of the system containing A 19 (4%), A 1 -b-S 6 (13%), and S 6 (3%) are presented. Solvent beads are omitted in the visualizations for clarity. The color scheme is the same as in Figure 1. At the bottom, the mean aggregation number, X n , polydispersity index (P) of solvophobic cores, and order parameter (ϕ) of the aggregates are shown as a function of the interaction parameter between the more solvophilic polymer and solvent a SW . Error bars are smaller than the symbol size. The scatter in the data is due to finite size effects in the simulations. Furthermore, Figure 6 presents X n and polydispersity of the solvophobic cores vs. a SW . The data show that the aggregation number does not have a simple relation to increasing a SW . While relatively large aggregates form for small a SW , a modest increase leads to a significant decrease in aggregate size, and further increasing a SW shows that the aggregate size increases. This can be understood if one considers that a small value of a SW = 30, i.e., good solvent solubility of the more solvophilic S 6 polymer, leads to a minority of them assembling in the aggregate. Increasing a SW makes the solvophilic polymer less soluble in the solvent and thus more of it assembles in the center of the aggregate. This may lead to promoting the merging of the solvophobic cores. This is shown as the increase in X n in Figure 6.
The polydispersity index of the aggregates in Figure 6 remains relatively insensitive to variations of a SW . For the larger a SW values, the increase in the polydispersity index results from the solvophilic homopolymers that remain incompatible with the solvophobic homopolymer aggregating into the aggregate. The presence of the solvophilic polymer in the aggregate imposes constraints, i.e., barriers for solvophobic beads to move in the system randomly and to find other solvophobic beads for favorable like-like interactions. This leads to more variation in the aggregate size.
The a SW dependence of the aggregate order parameter curve presented in Figure 6 shows an increasing trend for the order parameter with increasing a SW , i.e., increasing solvophobicity of the solvophilic polymer. This is a direct consequence of the decreasing compatibility between the solvophilic beads, solvophilic homopolymer, and solvent beads, which promotes phase separation.

Theoretical Arguments
Formation of micelles by amphiphilic di-block copolymers in aqueous solutions is wellunderstood [63][64][65]. For example, for amphiphilic macromolecules in water, micellization at low temperatures is driven by a reduction in free energy by the entropy associated with the water molecules. In general, the total Gibbs free energy can be written as a sum of all the chemical potential energy terms in the solution, including those of the solvent, singly dispersed di-block copolymer, and aggregates of different sizes. In the simplest picture, micellization is driven by the free energy change, ∆G g , to transfer a single amphiphilic di-block copolymer molecule from the solvent to an aggregate of size g: k B T∆G g = ∆µ g = ∆µ g,SO + ∆µ g,I + ∆µ g,SI where ∆µ g,SO is a negative contribution arising from removing the solvophobic polymer segments from having contact with the solvent (including the intramolecular/intermolecular solvophobic-solvophobic interactions), ∆µ g,I is a (positive) interface contribution, and ∆µ g,SI is a (negative) contribution from the intramolecular/intermolecular solvophilicsolvophilic interactions. We note that in this model, the interface free energy contribution comes from the change in free energy per unit area when the aggregate molecules are transferred from their bulk liquid phase to either a core-corona interface or an interface with the solvent. In contrast to the simple core-shell aggregate formation in single species di-block copolymer systems, multicore and multicompartment aggregation that emerges in multicomponent systems, such as that considered in the present DPD simulations, does not have a solid theoretical underpinning. At the microscopic level, the problem stems from the appearance of multiple competing length scales due to the many different interactions in the system. Although general formulation is lacking, the case of amphiphilic ABC copolymers has been treated with self-consistent field theory [34]. It allows one to write the total free energy as a sum of Flory-Huggins-type internal energy contribution and configurational entropy terms. In such systems, multicore assemblies can appear with a range of interaction parameters and chain lengths.
Within the thermodynamic approach, for multicore micellization, we can write the total free energy of the multicomponent system, see Equation (8), as: where the sum over i is over all the molecular species (solvophobic/solvophilic homopolymers and di-block copolymer) with aggregation numbers g. Detailed evaluation of the role of various configurational and entropic contributions to ∆G g,tot is complicated and outside the scope of the present work. Instead, we will focus on the role of the various dominant terms from the non-bonded interactions in Equation (9) between the different species to qualitatively explain the DPD simulation data.
To this end, we focus on the multicore aggregate case and divide it into six parts, as schematically shown in Figure 7. The multicore aggregates, such as those formed in this work, are composed of the following regions: • (i) Cores filled with solvophobic-free homopolymer A 19 with a negative free energy contribution, ∆µ g,SO (A 19 ), which dominates here, favoring aggregation of A 19 into the core. • (ii)-(iii) An interface between the core surface and the solvophobic A blocks of the A 1 -b-S 6 copolymers that surround the cores. Here, the free energy has a positive contribution, ∆µ g,I (A − W), for the interface between the solvophobic core and the solvent, and a negative contribution, ∆µ g,SO (A), resulting from removing core-solvent contacts by the A 1 -b-S 6 copolymers. • (iv)-(v) A corona composed mainly of the solvophilic copolymer blocks and the S 6 homopolymer. The free energy contribution is negative and comprises of ∆µ g,SI (A 1 − b − S 6 ) and ∆µ g,SI (S 6 ). This allows the solvophilic S beads to form the large corona. • (vi) The external interface between the aggregate and solution mainly composed of solvophilic homopolymers with interface energy ∆µ g,I (S − W). • (i) Cores filled with solvophobic-free homopolymer A19 with a negative free energy contribution, ∆ , (A ), which dominates here, favoring aggregation of A19 into the core. • (ii)-(iii) An interface between the core surface and the solvophobic A blocks of the A1-b-S6 copolymers that surround the cores. Here, the free energy has a positive contribution, ∆ , (A − W), for the interface between the solvophobic core and the solvent, and a negative contribution, ∆ , (A), resulting from removing core-solvent contacts by the A1-b-S6 copolymers.  First, we will qualitatively describe the multicore aggregate formation, as seen in the DPD simulations data. The multicore assemblies form because of the dominance of the term ∆ , (A ), resulting from the significant hydrophobicity of the A beads in A19. This leads to the A19 forming the inner core (i) in Figure 7. Following this, due to the energy cost, ∆ , (A-W), from the interaction between the hydrophobic A beads and the solvent W, the A1-b-S6 assemble on the surface of the core to shield the cores. The copolymer assembly on the core surface results in the formation of one layer of A1 blocks on the core surface. This forms the interface region (ii), (iii) due to the negative contribution from ∆ , (A). The S6 blocks and S6 homopolymer chains thus extend outwards, forming the solvophilic corona (iv), (v) and the small micelles. Eventually, the repulsion between S and W (S-S interaction is stronger than S-W) drives the small micelles to merge because of the large negative ∆ , (S). This finally drives the system to self-assemble to a stable multicore aggregate with an interface with the solvent (vi).
Next, we consider the influence of changing the solvophobicity interaction. By decreasing aAW, the solvophobe-solvent interfacial tension decreases, and the assemblies can grow in size, see Figure 5. The multicore aggregate structures are still observed for aAW between 115 and 65. When the A-W repulsion becomes weak enough, the A19 cores First, we will qualitatively describe the multicore aggregate formation, as seen in the DPD simulations data. The multicore assemblies form because of the dominance of the term ∆µ g,SO (A 19 ), resulting from the significant hydrophobicity of the A beads in A 19 . This leads to the A 19 forming the inner core (i) in Figure 7. Following this, due to the energy cost, ∆µ g,I (A-W), from the interaction between the hydrophobic A beads and the solvent W, the A 1 -b-S 6 assemble on the surface of the core to shield the cores. The copolymer assembly on the core surface results in the formation of one layer of A 1 blocks on the core surface. This forms the interface region (ii), (iii) due to the negative contribution from ∆µ g,SO (A). The S 6 blocks and S 6 homopolymer chains thus extend outwards, forming the solvophilic corona (iv), (v) and the small micelles. Eventually, the repulsion between S and W (S-S interaction is stronger than S-W) drives the small micelles to merge because of the large negative ∆µ g,SI (S). This finally drives the system to self-assemble to a stable multicore aggregate with an interface with the solvent (vi).
Next, we consider the influence of changing the solvophobicity interaction. By decreasing a AW , the solvophobe-solvent interfacial tension decreases, and the assemblies can grow in size, see Figure 5. The multicore aggregate structures are still observed for a AW between 115 and 65. When the A-W repulsion becomes weak enough, the A 19 cores shielded by the copolymer become unstable as the (positive) surface energy term ∆µ g,I (A-W) (not shown in Figure 7) decreases and destabilizes the multicore configurations. This leads to the formation of single-core assemblies of one core of A 19 and a double-layer corona made first by A 1 and second by S 6 blocks of A 1 -b-S 6 (see Figure 4) or soluble particles of A 19 with a sparse coating of A 1 -b-S 6 (see Figure 4, configurations corresponding to a AW = 55 and 35).
Finally, we consider the case where the solvophobicity of the homopolymer S 6 is increased by increasing the parameter a SW from its reference value of 30 (data of Figure 6). As expected, increasing a SW leads to aggregation of the S 6 chains in the emerging aggregates that now contain cores made of both the S 6 and A 19 molecules, as the large repulsive a SW prevents them from mixing within the aggregates. The negative contributions from both ∆µ g,SO (A 19 ) and ∆µ g,SI (S 6 ) are responsible for the formation of two types of solvophobic cores. Due to the large positive contributions, ∆µ g,I (A-W) and ∆µ g,I (S-W), the cores of A 19 and S 6 are shielded by the copolymer via blocks A 1 and S 6 , respectively. The former leads to small-core-double-shell aggregates, while only irregularly shaped aggregates of S 6 appear due to the short length of the A 1 blocks. Changing the di-block copolymer block length ratio and increasing its degree of polymerization, the cores of S 6 can be stabilized through small-core-double-shell aggregate formation in the overall multicoremulticompartment assembly.

Conclusions
This work has investigated self-assembly in a three-component mixture of linear polymers composed of a solvophobic and a solvophilic homopolymer species, and an amphiphilic block copolymer species in solvent using DPD simulations. The study conditions considered fixed chain length and concentration. We showed that the linear threecomponent mixture self-assembles to multicore and multicore-multicompartment aggregates, in addition to standard core-shell particles. The assembly response and transitions elucidated the role of mutually competing interactions when the interactions in terms of the solvent selectivity of the components were varied. As our main result, we found that the solvophobic homopolymer assembled to finite size aggregates, that due to their high interfacial tension with the solvent, were covered and stabilized by the copolymer chains and the solvophilic homopolymer, leading to small core-shell micelles. Depending on the solvophobic polymer solvent selectivity, the small micelles either aggregate further to form multicore assemblies or fuse to single-core aggregates (in finite size, standard core-shell polymer micelles). By variation of the polymer-solvent interaction for both the solvophobic and solvophilic polymer segments, we mapped the range over which multicore aggregates assemble. Finally, we discussed the conditions for the formation and stability of multicore-multicompartment aggregates and presented qualitative free-energy arguments governing the formation of such assemblies.
The work revealed some important practical guidelines for multicore assembly formation in the present system. The main requirements are a sufficient solvophobicity of the solvophobic component and a relatively high solvophobicity difference between the solvophilic and solvophobic polymer components. Additionally, there should be two different types of immiscible solvophobic homopolymers and commensurability of each block of the copolymer with one of the homopolymers. Furthermore, decreasing the solvophobicity of the solvophobic polymer segments leads to two transitions in the aggregate structures: first, the multicore aggregates transition to phase separation of the solvophobic and the solvophilic polymer segments (single-core assemblies), and then into mixing of the components (partially soluble particles). We also demonstrated that increasing the solvophobicity of the solvophilic homopolymer changes the aggregate structure from multicore to multicore-multicompartment assemblies. In particular, in our model system, spherical multicore assemblies form for a AW above 95 and multicore-multicompartment assemblies should have a SW larger than 45.
The structural transitions presented here, and their relation to direct polymer chemistries via the interaction parameters, enable extracting guidelines for the multicore and the multicore-multicompartment assembly. Naturally, for real polymeric systems, practical factors such as polydispersity, the effect of solution impurities, and temperature need to be considered. However, the presented study provides a strong guideline for directing practical polymer system realizations for multicore and multicore-multicompartment assemblies by revealing trends and control features driving the assembly morphologies and their transitions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/polym13132193/s1, Figure S1: Time evolution of the aggregate formation, Figure S2: Final simulation configuration from the simulation system (100R c ) 3 in size, Figure S3: Time evolution of the dimensionless order parameter ϕ, Figure S4: Time evolution of the polydispersity index versus time, Table S1: Morphological characteristics of the aggregates formed in the (100R c ) 3

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.