Shape characteristics of the aggregates formed by amphiphilic stars in water: dissipative particle dynamics study

We study the effect of the molecular architecture of amphiphilic star polymers on the shape of aggregates they form in water. Both solute and solvent are considered at a coarse-grained level by means of dissipative particle dynamics simulations. Four different molecular architectures are considered: the miktoarm star, two different diblock stars and a group of linear diblock copolymers, all of the same composition and molecular weight. Aggregation is started from a closely packed bunch of $N_{\text a}$ molecules immersed into water. In most cases, a single aggregate is observed as a result of equilibration, and its shape characteristics are studied depending on the aggregation number $N_{\text a}$. Four types of aggregate shape are observed: spherical, rod-like and disc-like micelle and a spherical vesicle. We estimate"phase boundaries"between these shapes depending on the molecular architecture. Sharp transitions between aspherical micelle and a vesicle are found in most cases. The pretransition region shows large amplitude oscillations of the shape characteristics with the oscillation frequency strongly dependent on the molecular architecture.


Introduction
Star polymers represent one of the simplest branched polymeric architectures utilising several arms in the form of linear chains linked to a central core. Branching imposes intramolecular constraints on the monomers and modifies mechanical, viscoelastic, and solution properties of such molecules in bulk and in a solution. A relatively recent review of the synthesis, properties and application of star polymers is given in [1].
Recently, one observes a growing interest towards amphiphilic and polyphilic star polymers. There are two main types of such star polymers, namely, the miktoarm and the diblock stars. In the miktoarm star, all the arms are of a homopolymer type but their properties differ from arm to arm. In the diblock star, each arm is a linear diblock copolymer itself. Both cases generate interest in both synthetic methodologies and in the self-assembly of such molecules in solution [2][3][4]. In particular, depending on the details of the molecular architecture, a number of morphologies are observed, such as: Archimedean tiling patterns and cylindrical microdomains for miktoarm star copolymers as well as asymmetric lamellar microdomains for diblock stars, which have not been reported for linear block copolymers, for more details see [4].
Applications of star polymers range from nano-and micropatterning [1,4] in gel and bulk state to drug delivery via their micellisation in solution (see, e.g. [5]). In the latter case, both the size and the shape of a micelle plays an important role from the requirements of low toxicity and efficient transport of the drug. For example, it has been reported that smaller nanoparticles (∼ 25 nm) travel through the lymphatic more readily than the larger particles (∼ 100 nm) and accumulate in lymph node resident dendritic cells [6,7]. Besides that, the shape of the particles also affect their bio-application: e.g., wormshaped filamentous micelles show less phagocytosis as compared to spheres, rod-like and red blood cell discoidal particles boasted longer blood circulation time than spheres, which reduced their clearance from the bloodstream (see [5] and references therein).
Besides experimental studies, the aggregation of polyphilic star polymers has been a subject of theoretical investigations [8] as well as computer simulations using various approaches. In particular, Monte Carlo simulations [9,10] performed on a lattice revealed a spontaneous formation of roughly spherical aggregates of diblock stars in a selective solvent. Significant changes in the micellar properties, such as the critical micellar concentration and aggregation number are reported upon changing the fraction of the solvophobic part of the diblock while keeping the total arm length constant.
Coarse-grained approaches, such as the dissipative particle dynamics (DPD) [11,12], make it possible to consider both a variety of molecular architectures and chemical compatibility of the particular groups, combined with the computational efficiency. As a result, one can study microphase separation driven effects in amphi-and polyphilic molecules, self-assembly, adsorption, etc. (see, e.g. [13][14][15][16]). With respect to the aggregation of the polyphilic star molecules, there are a number of studies of polymerosomes [17], ABC star block copolymers [18][19][20], as well as related systems, e.g., a mixture of diblock and homopolymers [21]. A number of micellar shapes have been observed, such as spherical micelles, wormlike micelles, hamburgers and others.
Therefore, coarse-grained DPD simulations can be considered as a useful instrument to study the size and the shape of aggregates formed of amphi-and polyphilic molecules of star architecture and, therefore, are capable of sheding more light on their applicability for micelle-based drug delivery systems. The aim of this work is to study the effect of the molecular architecture of amphiphilic star polymers on the type of aggregates they form in water. We consider three types of stars: the miktoarm star and two different diblock stars. These are compared against the equivalent set of linear diblock copolymers. The molecular weight and the fraction of hydrophilic monomers are chosen to be the same in all cases. The outline of the study is as follows. The simulation approach and the properties of interest are described in section 2, the results are presented in section 3 followed by conclusions in section 4.

Simulation approach and properties of interest
We employ the mesoscopic method of DPD [11,12] simulations, which describes polymer molecules at the level of coarse-grained beads each representing a fragment of the chain. The same applies to the water solvent, in which case a single solvent bead is assumed to contain several molecules of water. All beads representing both a polymer and water are spheres of the same diameter which provides the lengthscale of the problem, whereas the energy scale is assumed to be * = k B T = 1, where k B is the Boltzmann constant, T is the temperature, and time is expressed in t * = 1 [22]. The monomers are connected via harmonic springs, which results in a force where F C i j is the conservative force resulting from the repulsion between i th and j th soft beads, F D i j is the dissipative force that occurs due to the friction between soft beads, and random force F R i j that works in pair with a dissipative force to thermostat the system. The expressions for all these three contributions

13802-2
Shape characteristics are given below [22] , v i is the velocity of i th bead, a is the amplitude for the conservative repulsive force. The dissipative force has an amplitude γ and decays with the distance according to the weight function w D (x i j ). The amplitude for the random force is σ and the respective weight function is w R (x i j ). θ i j is the Gaussian random variable and ∆t is the time-step of the simulations. As was shown by Español and Warren [12], to satisfy the detailed balance requirement, the amplitudes and weight functions for the dissipative and random forces should be interrelated: σ 2 = 2γ and w D ( We will denote the beads representing water as of type A. The compressibility of such coarse-grained solvent at a number density of beads equal to ρ = 3 matches that for water at normal condition if the repulsion parameter a in equation (2.3) is chosen equal to a AA = 25 [22]. Similarly, the level of hydrophilicity of a polymer fragment composed of beads of type P can be controlled by the value a PA for the repulsion interaction between beads P and A, where the difference a PA − a AA is related to the Flory-Huggins parameter [23]. In our study, we assume the hydrophilic fragments of a polymer chain to be composed of the beads A, the same as for water, whereas the hydrophobic fragments -of the beads B, characterised by the repulsion parameter a AB = 40, the value already used in a number of studies [14,15,24]. The other parameters are as follows: γ = 6.75, σ = 2γ = 3.67 and the time-step ∆t = 0.04. Another parametrisation of the repulsive potential in DPD studies aimed at polymerosomes can be found in [17].
We use the simulation box with the periodic boundary conditions. Its linear dimension L is chosen at least twice the linear dimension of the largest aggregate formed out of N mol molecules each containing n beads. This requirement appears from the need to restore the integrity of an aggregate if it is split across one or several walls of the simulation box. A rough estimate for L can be achieved by assuming the formation of a single spherical aggregate with the diameter D a . As far as each bead on average occupies the volume of 1/ρ, the volume of an aggregate is nN mol /ρ = πD 3 a /6, hence D a ∼ (6nN mol /πρ) 1/3 and, therefore, the condition for the linear dimension of the simulation box reads L > 2D a . For highly aspherical aggregates, the value of L should be accordingly increased.
To identify the existing aggregates at each time instance, we use the following algorithm: (1) for each (i mol , j mol ) pair of molecules, find the number n c of close bead-bead intermolecular contacts evaluated for hydrophobic beads (type B), and if n c > n min , then register the link between the molecules i mol and j mol ; (2) using the linking list for molecular pairs, build neighbour lists for linked molecules; (3) apply a "snowball" cluster identification algorithm, in which the neighboring molecules are added to the existing ones within each aggregate until no neighbours remain; (4) final information contains: the lists of molecules that belong to a particular kth aggregate, and the index array to determine the aggregate to which any molecule belongs.
A close bead-bead contact is registered when the distance x i j between them is less than 1.5, and the threshold number of close contacts n min is chosen to be equal to 10.
After the aggregates have been identified, the lists of molecules for each of it is used to rejoin the aggregates split by the periodic boundary conditions. Then, we proceed to the analysis of their size and shape properties. The radius of gyration and the shape characteristics of each aggregate are derived from the components of the gyration tensor Q defined as in [25,26]: (2.6)

13802-3
Here, N is the total number of beads in an aggregate, x α n denotes the coordinates of nth bead: n are the coordinates of the center of mass for the aggregate. Its eigenvectors define the axes of a local frame of the chain and the mass distribution of the latter along each axis is given by the respective eigenvalue λ i , i = 1, 2, 3, respectively. The trace of Q is an invariant with respect to rotations and is equal to an instantaneous squared gyration radius of the chain Here, the average over three eigenvalues,λ, is introduced to simplify the following expressions. The hydrodynamic radius R h is found according to the following expression where the averaging is performed over all pairs of beads that belong to the aggregate.
The shape properties of aggregates are characterised by the asphericity A (sometimes also referred to as "the relative shape anisotropy") and prolateness S [27][28][29][30], as well as the shape descriptor B introduced in this study: For a spherical shape: The convenience of the use of the shape descriptor B is that its range is symmetric spanning from −1 to 1 and it reverses its sign when the shape changes from the disc-like to the rod-like.
At each time instant t , the largest aggregate, or "the giant component" in percolation terminology, is identified and the set of its characteristics, R 2 Based on their behaviour with time t , we estimate the time t equ needed for the aggregate to forget its initial state and to start displaying stationary oscillatory behaviour in terms of its shape properties. The conservative estimate reads t * ini ∼ 8000 (200 · 10 3 DPD steps). Then, the average values 〈R 2 g 〉, 〈R −1 h 〉, 〈A〉, 〈S〉 and 〈B 〉 are evaluated by simple averaging of the instant values over the time interval t equ < t < t fin , where t * fin ∼ 10t * equ is the total duration of each run. One should note that the other type of averaging is also possible, where the nominators and denominators in equations (2.9) are averaged separately [27][28][29][30][31][32]. Besides the simple averages, we also consider the histograms of probability distributions for the shape characteristics, where the p(B ), for the shape descriptor B , is found to be the most informative.

Aggregates and their properties
We consider four types of molecular architectures, all shown in figure 1. All of them can be interpreted as four linear diblock copolymers, each of eight hydrophobic beads (type B) and of two hydrophilic beads (type A) bonded in a different way. In particular, architecture (a) is just a set of four non-bonded diblocks and serves as some kind of a reference system. Architecture (b) represents the bonding of four  The aggregation transition in a solution of polymers is of much interest and is the subject of numerous recent works [10,[33][34][35]. In this study, we consider the setup that mimics a drop of a highly concentrated polymer solution immersed into water. The drop is modelled via a bunch of N mol molecules closely packed inside a central part of the simulation box. Namely, a central bead of each star is positioned randomly with all its coordinates falling into the interval of [L/3; 2L/3]. Each arm is then generated as a random walk. This algorithm results in inevitable overlaps of beads in the initial configuration. The overlaps, however, do not produce forces of large magnitude, as it would occur in atomic simulations, due to the soft nature of the interaction potential (2.3). In the course of simulation, this bunch of molecules looses the memory concerning its initial state during the time t equ . After that it starts to display stationary oscillations in its shape properties. The values for t ini and for the simulation duration t fin have already been provided in section 2.
In most cases being studied, the initial setup equilibrates into the single aggregate state. This, however, might be the result of a finite system size and also a finite simulation time. In a case of a single aggregate, the aggregation number N a coincides with the total number of star molecules N mol in the solution.
The increase of the latter results in a growth of the average aggregate size characterised by the estimates for the gyration and hydrodynamic radii, as shown in frames (i) and (ii) of figure 2, respectively. We find the respective data points obtained for the molecular architectures (a)-(c) to be very close for both size properties being considered, especially at larger N a . Therefore, the difference in intramolecular binding in these cases has a minor impact on the average aggregate size. The aggregate size is dominated by a collapse of the hydrophobic part of stars characterised by the scaling laws 〈R 2 g 〉 1/2 , 〈R −1 h 〉 −1 ∼ (n h N a ) 1/3 , where n h is the number of hydrophobic beads in each star. Data points obtained for the case (d) are found to be higher than the respective values for the cases (a)-(c), especially in frame (i). This deviation, as well as local deviations from the scaling law observed at specific values of N a in cases (a)-(c) are related to the peculiarities of the shape oscillations discussed below.
On a classification note, we found four types of aggregate shape: spherical, rod-like and disc-like micelles, as well as a spherical vesicle, all visualised in figure 3. For most molecular architectures considered here, the increase of N a results in the following sequence of shapes: spherical micelle → aspherical micelle (rod-like or oscillating between the rod-and disc-like) → spherical vesicle.
This sequence can be understood in terms of the competition between the enthalpic and entropic contributions to the free energy related to the presence of the aggregate. The enthalpic contribution is minimal when the number of contacts between the hydrophobic and hydrophilic beads are minimised, which is achieved in the case of a spherical shape. The entropic contribution is minimised when the hydrophobic chains within the aggregate display coil-like conformations. Both factors do not compete in the case of a spherical micelle until its radius R a is smaller than the average end-to-end distance r e of the hydrophobic chains in a coiled state, see frame (i) in figure 4. With the growth of N a , there are two options: sperical micelles with R a > r e shown in frame (ii); or aspherical micelles shown in frame (iii) of the same figure. For the case (ii), the need to fill-in the center of the micelles inevitably leads to stretching some of the hydrophobic chains beyond the r e distance (these chains are shown in red), thus reducing the entropy (i.e., increasing the free energy). The effect strengthens with the increase of N a and, hence, of R a . Therefore, at a certain value of N a , the aspherical aggregate becomes more favourable, when the penalty of aspherical shape is compensated by the gain from a coiled state for all chains, as shown in frame (iii).
A further increase of N a forces a higher asphericity of this shape and a spherical vesicle with the internal void, shown in frame (iv), becomes more favourable. It combines both overall spherical shape and coiled conformations for the chains within a hydrophobic shell, which is about r a for a single layer vesicle.
Let us consider now how the details of the molecular architecture affect the "phase boundaries" for  showing rather sharp peaks centered around zero for N a = 20 with the distribution getting broader and shifting towards larger values of B for larger N a < N * a . In this interval both disc-like (B < 0) and rod-like (B > 0) shapes are found, albeit the disc-like tendencies are rather weak not reaching the values below −0.25. Above the transition, N a > N * a , the distribution is back to a single sharp peak centered around zero indicating a spherical vesicle. Its shape practically overlaps with the one for the spherical micelle (N a = 20). The case (b) shows similar tendencies for the probability distribution p(B ), but the effect is much weaker. The case (d) is rather different. At N a = 20, the distribution p(B ) shows a series of peaks at relatively large values of B indicating that the spherical shape is not the single favourable one, presumably due to specific intramolecular constraints of the case (d). The tail of the distribution which spreads towards the larger values of B persists even for the largest value of N a = 120 being considered. It lifts the simple average 〈B 〉 to the nonzero values, as seen earlier in figure 5.
The The case (d) is characterised by fast oscillations between the weakly disc-like and highly rod-like shapes, resulting in essentially non-zero simple averages for this case in figure 5.
One should note that the aspherical shapes of aggregates observed for the amphiphilic stars in this study, are reminiscent of the worm-like prolate and burger-like oblate multicompartment micelles observed for the case of the ABC miktoarm stars [18,20,36]. In the latter case, however, these shapes are stabilised by adding the third component, whereas for the two-component amphiphilic stars they are prone to large scale oscillations observed in our work. It is plausible to suggest that these oscillations occur due to a delicate balance between the enthalpic and entropic contributions, whereas the addition of the hydrophobic drug agent [19] may stabilise the aggregate shape. This will be a subject of subsequent studies.

Conclusions
The simulation study presented in this work reveals strong dependencies of the shapes of aggregates formed by amphiphilic stars in water on the details of their molecular architecture. It has a prospect of the application in the drug delivery systems, where both the size and the shape of the aggregate is known to play an important role for the flow of the enveloped agent through the vessels.
Four molecular architectures have been examined: (a) four disjoint linear diblocks, (b) asymmetric miktoarm polymer, (c) diblock star 1 (hydrophilic parts pointing outwards) and (c) diblock star 2 (hydrophilic parts next to a central bead). For all cases, the same general sequence of shapes is found with an increase of the aggregation number, namely: spherical micelle, aspherical micelle and a spherical vesicle. The "phase boundaries" between these are found to depend on the details of the molecular architecture. For the case (a)-(c), the transformation between a spherical and aspherical micelle occurs gradually, whereas the transition from an aspherical micelle into a spherical vesicle is in a form of a sharp transition. In the case (b), aspherical micelle is less stable and transition to a vesicle occurs at a lower aggregation number. The case (d) is characterised by gradual transitions between all the shapes.