Modelling realistic microgels in an explicit solvent

Modelling realistic microgels in an explicit solvent

F. Camerin N. Gnan L. Rovigatti E. Zaccarelli

Thermoresponsive microgels are polymeric colloidal networks that can change their size in response to a temperature variation. This peculiar feature is driven by the nature of the solvent-polymer interactions, which triggers the so-called volume phase transition from a swollen to a collapsed state above a characteristic temperature. Recently, an advanced modelling protocol to assemble realistic, disordered microgels has been shown to reproduce experimental swelling behavior and form factors. In the original framework, the solvent was taken into account in an implicit way, condensing solvent-polymer interactions in an effective attraction between monomers. To go one step further, in this work we perform simulations of realistic microgels in an explicit solvent. We identify a suitable model which fully captures the main features of the implicit model and further provides information on the solvent uptake by the interior of the microgel network and on its role in the collapse kinetics. These results pave the way for addressing problems where solvent effects are dominant, such as the case of microgels at liquid-liquid interfaces.

microgels, molecular dynamics, dissipative particle dynamics, polymers,explicit solvent


In recent years microgels — colloidal-scale polymer networks — have emerged as a popular model system in condensed matter physics[1] thanks to their colloid/polymer duality[2]. The combination of colloidal properties and responsiveness to external stimuli is the key for their appeal for both applications and fundamental science[3]. Among microgels, the most widely studied are those based on Poly(N-isopropylacrylamide) (PNIPAM), a thermoresponsive polymer able to swell and deswell reversibly as a result of temperature changes. When PNIPAM chains are crosslinked with bisacrylamide (BIS), microgel particles can be prepared in a range of sizes of  nm by standard synthesis methods [4], and even reach much larger scale (up to m) with microfluidic techniques [5]. These particles undergo a Volume Phase Transition (VPT) in water at a temperature of , from a swollen state at low temperatures to a collapsed one at high temperatures. This swelling-deswelling transition is fully reversible and can be exploited to tune the size of the particles in situ. The VPT is completely controlled by the polymer-solvent interactions, echoing the coil-to-globule transition of linear PNIPAM chains in water[6]. As a matter of fact, the role of water is highly relevant, as the VPT originates from changes in the hydrophilic/hydrophobic character of the interactions of the polymer with the solvent upon temperature variations.

Experimental work on microgels has enormously increased in the last couple of decades, and a comparison of experimental data with theory has been possible thanks to the use of the classical Flory-Rehner theory of swelling[7]. On the other hand, microgel simulations have been less abundant due to the complex, multi-scale nature of the particles. So far, most efforts have relied on the use of unrealistic networks, often based on ordered, diamond-like topologies, in which all strands are of equal length[8, 9, 10, 11, 12]. Only a few of these approaches have explicitly considered the role of the solvent[13, 14, 12].

Recently, we have introduced a novel method to synthesize realistic microgel particles in silico through the assembly of fully-bonded, disordered networks with arbitrary topology[15, 16]. In this approach we initially consider the self-assembly of a mixture of patchy particles, respectively bivalent and tetravalent, to mimic monomers and crosslinkers. To retain a spherically-shaped network, the mixture is confined within a sphere of a given radius. Fully-bonded configurations are obtained by introducing a swapping mechanism that makes it possible to equilibrate the system even at the strong attractions required to maximize the bonding. In this protocol there are two parameters controlling the topology of the resulting network: the concentration of crosslinkers and the confinement radius. Thus, more compact and homogeneous networks are obtained in presence of a large number of crosslinkers and/or for a tight confinement, while looser and more heterogeneous microgels can be produced with a smaller amount of crosslinkers and a very weak confinement. A thorough discussion on how the internal structure of the microgels depends on these parameters can be found in Refs.[15, 16].

Once the network is assembled, we replace the patchy (reversible) interactions with permanent bonds by adopting the classical Kremer-Grest bead-spring model for polymers[17] to preserve the network topology throughout the course of the simulation. In order to reproduce the swelling behavior, it is possible to incorporate in the model an attractive potential that has been shown to capture the variation in polymer-water interactions upon changing temperature. With this approach, the solvent is implicit and the solvophobic potential accounts for it within the thermodynamic properties of the system in an effective way. This implicit solvent model was shown to be able to faithfully reproduce swelling data of individual microgels measured with Dynamic Light Scattering experiments [15]. Even though the use of an explicit versus an implicit solvent model[18, 19] should give identical results in terms of equilibrium properties, there are a number of features that cannot be correctly captured and/or described by an implicit model. In particular, the kinetics of swelling and deswelling will depend on the presence of the solvent and on how it is modelled. Besides that, there are situations of fundamental and practical interest in which an explicit solvent will dramatically affect the picture. For instance, to model a system at a liquid-liquid interface, it is necessary to take into account the presence of the two different media in order to capture effects related to the surface tension [20, 21].

In order to be able to handle these situations, here we take the implicit-solvent model of Refs.[15, 16] and extend it by developing an explicit solvent description that accurately predicts the swelling behavior of microgel particles. We use the swelling properties exhibited by the implicit solvent model, which has been shown to faithfully reproduce the experimental results, as reference data to calibrate the explicit-solvent parametrisation. By comparing the swelling ratio as a function of temperature and the microgel density profile and form factor with and without solvent, we are able to discriminate among different solvent models and choose the explicit description that works best. In particular, we intend to model a generic solvent that ensures that the key properties of microgel colloids are accurately reproduced rather than to provide a systematic and exhaustive study on the influence of the system parameters on the properties of the particle. We further test the robustness of our approach by repeating the analysis for microgels generated with different topologies and confinement radii. Once established our explicit model, we first look at the arrangement of the solvent inside the microgel across the volume phase transition, and then study the kinetics of the deswelling. Overall, our results open up the possibility to obtain more and more realistic descriptions of microgels, thanks to which it will be possible to tackle exciting problems in which the explicit role of the solvent plays a crucial role [21, 22, 23, 24].


Swelling behavior

Figure 1: Microgel swelling curves. Radius of gyration across the VPT transition for (a) the implicit model, ; (b) the explicit LJ solvent with LJ monomer-solvent interactions at a solvent density ; (c,d) explicit solvent with monomer-solvent interactions at and , respectively; (e) DPD simulations where the microgel is modeled as a bead-spring polymer network. All curves report the gyration radius as a function of the parameter controlling the solvophobic interactions in each model: (a) , (b) , (c-d) and (e) .
Figure 2: Effect of microgel topology and solvent arrangement. Swelling curves for the implicit- (full line) and explicit-solvent models that best reproduce the swelling behavior, namely MD simulations with at (dashed lines) and DPD simulations (dotted lines) for (a) a loose microgel () and (b) a more compact microgel (). Corresponding microgel snaphots are also shown. Symbols refer to state points in explicit solvent simulations (MD: circles, DPD: triangles) for which further analysis is provided in the next sections, whereas similar colors/shapes refer to similar swelling degrees between the two explicit solvent models. Panels (c.I-c.III) display a central slab of the simulation box for three different values of , respectively corresponding to the swollen state (c.I), a state very close to the VPT (c.II) and the collapsed state (c.III). The arrangement of solvent (blue spheres) within/around the polymer network (red spheres) depends on . For visual clarity, only half of the solvent particles are shown.

We start by discussing the swelling behavior of microgels in the presence of an explicit solvent as compared to the reference case of the implicit model , Fig. 1(a) (see Methods), discussed in Ref.[15]. To this aim, we perform simulations of an individual microgel assembled with a rather loose topology (using a confining radius ) in different solvents. In particular, we make a comparison between “atomistic” and coarse-grained solvent representations by employing Molecular Dynamics (MD) and Dissipative Particle Dynamics (DPD) simulations, respectively (see Methods). In the former type of approach, we first need to adjust the solvent-solvent interactions, for which the most natural choice is to use a Lennard-Jones (LJ) potential. Next, we address the choice of the monomer-solvent (ms) interactions: these are crucial to describe the swelling transition, because they control the contraction or extension of the polymers chains in the solvent environment. To discriminate between different models and identify the best possible candidate, we explore multiple ms potentials and compare them to the implicit solvent case. The choice of the solvent density allows to tune the pressure exerted by the solvent on the polymer network thus determining the swelling range of the microgel particle, as discussed below.

Similarly to solvent-solvent interactions, a straightforward choice for the monomer-solvent ones is the LJ potential[25] where, by varying the energy minimum , we control the polymer-solvent affinity. In this way, we obtain the swelling curve reported in Figure 1(b), where the radius of gyration of the microgel is shown as a function of : by decreasing this parameter (with respect to solvent-solvent interaction, which sets the energy scale), the polymer-solvent interactions are less favoured than solvent-solvent ones, giving rise to a reduction of the microgel size. However, an unphysical increase of is observed for : under this condition, both terms in the LJ potential go to zero, i.e. the microgel feels neither attraction nor repulsion with the solvent. Consequently, the network relaxes as the external pressure on the polymer network vanishes, and the microgel swells again, maximizing its configurational entropy.

Such behaviour clearly indicates the unsuitability of the LJ potential to mimic the solvent-monomer interactions. Consequently, the next step is to use a potential in which the attractive term can be tuned arbitrarily without affecting the short-range repulsion. To this aim, we adopt the model, defined in Eq. (3), where the repulsion remains unchanged while the attractive contribution, controlled by the parameter , is varied. The swelling behavior of the microgel obtained with this model is reported in Fig. 1(c,d) for two representative solvent densities. The swollen-to-collapsed transition is well reproduced in both cases.

So far, we have assessed the “atomistic” type of solvent. We further examine the possibility to use a coarse-grained solvent by means of DPD simulations, which correctly reproduce hydrodynamic interactions at long times[26]. In order to establish a meaningful comparison with the implicit solvent case and avoid unphysical crossing of the chains, we retain the bead-spring model for monomer-monomer interactions and we apply the DPD treatment only to monomer-solvent and solvent-solvent interactions.

The results of DPD simulations, for the parameters specified in Methods, are reported in Fig. 1(e). In this case, the VPT transition is modulated by the monomer-solvent repulsion quantified by the parameter in Eq. (4): for small values of the microgel is swollen, while it contracts when increases. We notice that is systematically larger at comparable swelling for MD-solvents than for DPD results, which, on the other hand, quantitatively reproduce the values obtained in the implicit solvent description. This is due to the softness of the DPD interactions which, contrarily to the MD treatment, do not introduce significant solvent-monomer excluded volume effects, thereby not affecting the microgel size.

In order to establish a correspondence between different models, we rescale the explicit solvent data onto the implicit one, where is the solvophobic parameter (see Methods). Figure 2(a) shows the normalized , where is the value of the gyration radius at maximum swelling, as a function of the effective swelling parameter . The latter corresponds to the solvophobic parameter of the implicit solvent simulations. We report the comparison for the two cases where the agreement is found to be fully satisfactory for all , namely the DPD and MD models. Of the latter, we consider only the case with the highest solvent density, , since deviations with respect to the implicit solvent case are observed with lower densities: the swelling range of the microgel would be shortened, as can be observed in Figure 1(c). Thus, it appears that, while is definitely superior to the simple LJ potential to model the VPT of the microgel, the density of the solvent particles is a key parameter in tuning the details of the transition: a lower density will have a smaller effect on the microgel, resulting in a more limited contraction with respect to the implicit solvent model. From now on we will discard the LJ potential and we will refer to MD simulations as those performed with the interaction. A similar effect can be obtained in DPD simulations by changing the cutoff radius and the interaction parameters of the conservative force, which represents the length scale in DPD and the size of the solvent beads (see Methods).

To verify the robustness of our protocol, we now repeat the above analysis on a microgel configuration assembled with a smaller confinement radius, . Fig. 2(b) reports the swelling behaviour of the more compact microgel for the DPD and MD models at the optimal solvent density identified above. Together with the data, we also report snapshots of the two microgels (insets) in their maximally swollen state, showcasing the very different topology of the networks. The good agreement between the rescaled swelling curves for both studied microgel configurations allows us to conclude that the developed models are robust and both can faithfully reproduce the swelling behavior observed with the implicit model[15]. Fig. 2(c.I-c.III) further highlights the arrangement of solvent particles inside the microgel for MD simulations at different values of across the VPT. The microgel remains very permeable to the solvent even close to the transition temperature, finally expelling it only in the fully compact state. In the next sections we focus on MD and DPD to study the effects of the solvent on the microgel structural features and on the kinetics of the volume phase transition.

Structural features of a loose microgel in an explicit solvent

Figure 3: Density profiles for a loose microgel configuration across the VPT. Monomer radial density profile for a microgel as a function of the distance from its center of mass. Full lines refer to the implicit-solvent model, while symbols are used for MD (circles) and DPD (triangles) simulations. Each sub-panel refers to a different swelling state as in Fig. 2(a).

We now discuss the structural features of the microgel at relatively large confinement, corresponding to the swelling curve in Fig. 2(a). First, we show results for the density profile of the microgel in Fig. 3, for several values of the swelling parameter across the VPT for both MD and DPD simulations. We find that, in general, both solvent models yield density profiles that are very similar to the implicit solvent case. This is particularly true for the swollen states, where the typical core-corona structure of the microgels is clearly distinguishable. Under these conditions, DPD simulations are even more accurate than MD ones in reproducing the results of the implicit model. When increases and the microgel becomes more compact, the difference between the three models becomes more evident. Specifically, as the microgel collapses MD simulations produces lower density profiles in the core region with respect to the implicit-solvent case at the same , while the DPD model generates more compact structures.

We notice that low density profiles exhibit a non-flat behavior in the inner core region of the microgel. These inhomogeneities, that are stronger for smaller microgels, can be removed out by averaging over independent topologies[15]. Here we do not perform such an average because we aim to compare the behavior of a given microgel configuration with and without solvent. Beyond the VPT the oscillations are suppressed by the higher density, and hence the profiles are much flatter within the core.

While density profiles provide real-space information on the microgel structure, they are not easily accessible in experiments, except for very recent super-resolution microscopy investigations[27, 28]. Instead, they can be indirectly obtained from fitting the form factors to the fuzzy sphere model[29]. The form factors can be measured by small angle neutron or x-ray scattering experiments. Thus, in contrast to density profiles, numerical can be used to make a direct comparison with experiments, without having to rely on fits to specific models. Indeed, while the fuzzy-sphere model correctly describes the core-corona structure, it does not take into account the presence of dangling chains in the outer corona shell[30, 15]. We thus directly evaluate the form factors of the microgel across the VPT and present them in Fig. 4 as a function of wavevector for the same values of swelling parameters used in Fig. 3.

Figure 4: Microgel form factors for a loose microgel across the VPT. as a function of . Full lines refer to the implicit-solvent model, while symbols are used for MD (circles) and DPD (triangles). Each sub-panel refers to a different swelling state according to Fig. 2(a).

We find that the use of an explicit solvent does not considerably alter the form factors with respect to the implicit solvent case for all values of the swelling parameters. As increases and the solvent quality decreases, shows an increasing number of oscillations which become more and more pronounced. Furthermore, the position of the first peak, which is related to the microgel overall size, shifts to larger and larger wavevectors, indicating the shrinking of the microgel. However, a subtle difference is present between the two types of employed models: while DPD results are perfectly superimposed to the implicit solvent case for all , the MD results are found to be always shifted to a slightly smaller -value with respect to them. This is a reflection of the overall microgel size, which is a bit larger for MD explicit-solvent simulations with respect to DPD and implicit solvent, due to stronger excluded volume effects, as evident from Fig. 1. We further notice that at relatively large wavevectors () the MD form factor systematically overestimates the DPD and implicit-solvent ones for intermediate and large values of . However, all curves superimpose again at , where a small peak is found, independently of the swelling parameter value. The latter corresponds to the monomer-monomer nearest-neighbour peak and is a feature associated to the excluded-volume interactions included in the bead-spring model for polymers and to the finite size of the simulated microgel. Indeed, for larger and larger microgel size, this peak would become more and more separated from the first one, allowing for a larger number of oscillations. In experiments, such a peak is not generally noticeable because of the soft intrinsic nature of the monomers. Thus, it is a limitation of the present modelling, which on very small length scales becomes inaccurate.

We now turn to analyze the solvent density profile inside the microgel. The normalized profile , where is the bulk solvent density, is shown in Fig. 5 as a function of the distance from the center of mass of the microgel. Clearly, the distribution reflects, as a mirrored image, the one of the microgel monomers. Indeed, when the core of the microgels becomes denser and denser, more and more solvent gets expelled. It is interesting to note that, beyond the VPT and except for the very collapsed states, a significant fraction of solvent is retained within the polymer network, even well inside the core region. At the VPT, which takes place at , the density of the solvent inside the core is larger than 50% of the bulk value. Finally, we notice that there seems to be a consistent trend of the MD solvent to be more excluded from the network region with respect to the DPD results, again a feature associated to the larger excluded volume of the MD model. However, the two models yield qualitatively very similar results and reinforce the common view that microgels, despite their inhomogeneous structure and dense core region, retain of water in their swollen configuration and still contain a large amount of water well beyond VPT, in qualitative agreement with the experiments results of Ref. [31].

Figure 5: Solvent density profiles for a loose microgel configuration across the VPT. We show the solvent density profile normalized with respect to the bulk solvent density , as a function of the distance from the center of mass of the microgel. Circles and triangles refer to MD and DPD solvent, respectively. Each sub-panel refers to a different swelling state according to Fig. 2(a).

Results for a more confined microgel

Figure 6: Microgel density profiles, solvent density profiles and form factors for a compact microgel across the VPT. (a-c) microgel density profiles as a function of the distance from the center of mass of the microgel; (d-f) solvent density profiles normalized with respect to the solvent bulk density as a function of ; (g-i) microgel form factors as a function of the wavenumber. Data are reported for a swollen state (), a state close to the VPT () and a compact state (). Full lines refer to the implicit solvent (), while symbols are used for DPD (triangles) and MD (circles). The insets in panels g and h show an enlargement of the high wavevector region where solvent-monomer excluded volume interactions induce an excess of signal for the MD data.

We now repeat the above structural analysis for a more compact microgel topology obtained with a smaller confining radius (), whose swelling curve was reported in Fig.2(b).

The density profiles of the microgel are reported in Fig. 6(a-c) for a few selected values of the swelling parameter and again for both MD and DPD explicit solvents. We find that the DPD model reproduces very well the implicit-solvent data, particularly for the more swollen conditions. When increases, the DPD monomer density in the core is slightly larger than for the implicit case. However, the corona profiles of the two microgel representations are identical. On the other hand, the MD solvent results underestimate the microgel density profile in the core and also display a different corona profile for all . If compared to the findings for the looser microgel configuration (Figure 3), the DPD solvent model behaves similarly for both types of networks and well reproduces the implicit model data in all cases. By contrast, the MD results present systematic differences with respect to the other two sets of data making the agreement not completely satisfactory. This is a consequence of the “atomistic” treatment of the solvent, which interacts via excluded volume with the polymer. Especially for compact microgels, when excluded volume becomes more and more relevant, these assumptions in the model may become unrealistic. Thus, while for looser networks both MD and DPD explicit solvents provide a good description of the microgel, for more compact microgels the DPD model has definitely the upper hand. This is also shown in the behavior of the solvent density profiles reported in Fig. 6(d-f). Again we find that the MD solvent is much more excluded from the interior of the microgels at all . On the other hand, we see that, notwithstanding the relative higher compactness of this microgel, a significant amount of solvent remains inside the core in the swollen states, being roughly 60% of its bulk value close to the VPT, in agreement with what found for the less confined microgel configuration and with experimental estimates[31].

The form factors, shown in Fig. 6(g-i), further confirm that DPD results are in good agreement with the implicit model ones. However, the MD outcomes display a clear shift in the peak position which is much more evident than for the looser configuration (see Fig. 4). In addition, we observe an excess of signal, highlighted in the insets of Fig. 6(g,h), at in swollen conditions, which is absent in the DPD and implicit solvent simulations. This difference occurs at a length that is roughly twice that of the monomer-monomer peak, thus being associated to monomers that are apart, i.e. with a solvent particle in between them. Such a feature is smeared out at increasing , when the microgel collapses and monomer-monomer interactions become dominant. We notice that the excess signal is not observed for the looser microgel as, at the same value, excluded volume interactions are far less important. Overall, this further shows that the MD model, while still acceptable for not too dense and open microgels, becomes more inaccurate for rather compact ones.

Collapse kinetics

Figure 7: Collapse kinetics. Radius of gyration as a function of time for a loose microgel () for (a.I) and (a.II) for implicit (, full line), MD (dashed and dotted lines) and DPD solvents (dashed lines); (b) cluster size distribution for (indicated as III in a.I) for implicit and DPD solvents. In order to improve statistics data are averaged over six different microgels configurations; (c.I-III) simulation snapshots for state points I-III (circles in a.I). Clusters are highlighted by different colors according to their size (as indicated in the color bar). Light grey monomers are either found in small clusters () or belong to the main network ().

After having established the explicit solvent models and having analyzed the properties of microgel and solvent particles in equilibrium for different values of the swelling parameters, we now turn our attention to the kinetics of collapse of the microgel in the presence of the solvent. Employing the same approach adopted in Refs. [32, 33, 34, 35] for linear polymers, we start from a swollen microgel in a loose configuration and perform a sudden quench to a different state. In particular, we examine two final states whose value of correspond to an almost fully collapsed state () and to a state close to the VPT (). We then assess whether the collapse transition is affected by the presence of the solvent by comparing the kinetics of the implicit-solvent model with that obtained using MD and DPD ones. Figure 7(a.I-II) shows the time evolution of the radius of gyration of the microgel for the three different types of simulations at two different . In all cases the curves reach at long time the same value of but, in these simulation conditions, the time taken to equilibrate is different, being faster in implicit solvent simulations compared to those of DPD and MD (the slowest). All curves display a sharp one-step collapse with no trapping phenomena in metastable states. This is qualitatively in agreement with experiments in which microgels with a similar core-corona structure to ours are subjected to an abrupt temperature jump from low (swollen state) to high temperature (globular state) [36].

In order to highlight the role of the solvent, we perform a cluster analysis to identify how the microgel structure evolves during the collapse. To this aim, we detect clusters of non-bonded monomers (see Methods), and calculate their size distribution for state points having the same but simulated with different models. Remarkably, we find the same cluster distribution for both implicit and DPD solvent, indicating that the solvent plays no significant role on the folding dynamics of the microgel, as shown in Fig. 7(b). To visualize the restructuring of the microgel following the instantaneous decrease in the solvent quality, snapshots of the microgel are reported in Fig. 7(c.I-III) for three different times. The microgel, while shrinking, first reorganizes by grouping monomers into small clusters (panel c.I). Each cluster is connected to the others via single or multiple links so that the structure, at an intermediate shrinking stage, displays a large number of holes and becomes increasingly inhomogeneous. As the shrinking proceeds, the clusters start to merge, becoming larger and larger in size and joining the main network (panels c.II-III). Finally, at long times, all non-bonded monomers are connected and only a single cluster is left. We stress that this pattern is also found for the implicit model simulations and for the more confined microgel (not shown). These results strongly indicate that the solvent plays a minor influence on the structure of the microgel during the collapse transition. Indeed, at each swelling stage, the microgel has a similar structure regardless of the solvent employed, suggesting that deswelling occurs via the same sequence of transient states. It would be interesting to compare these findings with more accurate solvent treatments such as Multi-Particle-Collision-Dynamics simulations[13, 14].


The tunable swelling of the microgel particles has been, since their discovery, one of the most relevant features of these colloids. Indeed, the opportunity to tune the particle volume fraction without changing their number density, but only the temperature, is a formidable advantage for experimental investigations. However, this poses a computational challenge in choosing a suitable model that best describes their swelling-deswelling transition. The recent assembly of realistic microgel networks in Ref.[15] correctly reproduces experimental density profiles and form factors through an implicit solvent treatment. However, the inclusion of the solvent grants additional information, such as the uptake of solvent within the polymer network or surface tension effects. For these reasons, in this work we have compared the implicit solvent results to explicit solvent ones by employing two common approaches to simulations that allows for an atomistic and a coarse-grained approach, namely MD and DPD. We found that we can reproduce the implicit solvent swelling behavior by tuning the monomer-solvent interaction potentials after having adjusted the solvent density. This stems from the fact that, when the solvent is treated explicitly, the external pressure exerted by the solvent needs to be adjusted. In DPD simulations, the same effect can be obtained by regulating the cut-off radius.

We considered two microgels differing in the degree of compactness, which can be obtained by different synthesis protocol[37] and/or by varying the number of crosslinkers. We found that, particularly when the network is denser, excluded volume interactions play a relevant role in the description of the microgel. Indeed, in the full MD simulations an additional peak in the structure appears at small length scales. At the same time, the internal density profile of the microgel is also affected, resulting in a less dense core and a modified corona behavior, which is more significant in the collapsed state. Despite reducing the size of the solvent may solve excluded volume issues, doing so would dramatically increase the number of particles required to observe the same swelling behavior, as the box size is fixed by the dimensions of the microgel particle.

By contrast, DPD results better describe the implicit model ones for both microgel density profiles and form factors, at all swelling conditions. Furthermore, the DPD model reproduces the behaviour of the radius of gyration of the implicit model at different swelling conditions in an almost quantitative fashion. We have also investigated to what extent the solvent penetrates into the microgel, and we found that in the MD simulations much less solvent is present in the interior of the network, whereas DPD results seem more realistic in comparison to experimental estimates. Indeed, we find that, in the swollen state, the network is completely hydrated, retaining more than 90% of the solvent (with respect to the bulk density) in the core of the microgel. Even above the VPT the microgel contains a large fraction of solvent, which is finally excluded only at very large , amounting to temperatures C according to the mapping established in Ref. [15] for PNIPAM microgels.

We also examined the collapse kinetics and assessed how the presence of the solvent affects it. We observed that, in the conditions we performed simulations, a slowing down of the collapse dynamics occurring for the more structured solvent (MD simulations) and to a smaller extent for the coarse-grained solvent (DPD simulations) with respect to the implicit simulations. However, we also found that the system, when compared at the same swelling degree (quantified by the radius of gyration of the microgel), always presents a similar structure, regardless of the model. In particular, at first the network becomes rather inhomogenous, with regions where monomers have clustered together and empty regions. Later on, the clusters merge together and become larger and larger, until the collapse is complete and the microgel is essentially a fully folded network. Such transient behavior, featured by the appearing of crumples, has also been observed previously in simulations[33, 34, 12]. The similarity between these results with those found for an implicit solvent treatment suggests that hydrodynamic interactions do not play a major role in the swelling-deswelling transition, which is instead mainly controlled by the quality of polymer-solvent interactions.

In summary, in this work we have established that DPD simulations with a coarse-grained solvent constitute the most suitable method to include explicitly a generic solvent in the simulation of a microgel colloid. Even though a partially satisfactory description can be also obtained with the use of an MD solvent, the atomistic description allows for the presence of significant excluded volume interactions that brings unphysical features in the model. On the other hand, DPD simulations do show a full agreement with the implicit model and provides a realistic description of the solvent arrangement within the network. Thus, our model of realistically assembled microgels in DPD explicit solvent opens up the possibility to tackle those phenomena where the physical presence of the solvent is crucial. In particular, our model may serve as a starting point to numerically investigate the so-called “Mickering” emulsions[38] and in the fascinating case of microgels at fluid-fluid interfaces[22, 39, 23, 24, 21]. Finally, a realistic description of how the solvent is trapped within the polymer network may lead to future advances in the field of drug delivery and controlled-release, and can provide further insights into its mechanism[40, 41].


Microgel assembly. The starting configuration of a microgel particle is prepared as in Ref. [15]. First, we produce a fully-bonded disordered network by self-assembling a binary mixture of bi- and tetravalent patchy particles with an applied spherical confinement. Once the network has assembled, we replace the patchy interactions with the classical bead-spring model for polymers[17], in which bonded monomers interact via the sum of a Weeks-Chandler-Andersen [42] (WCA), , and a Finite-Extensible-Nonlinear-Elastic[43, 44] (FENE), , potentials:


with an adimensional spring constant and the maximum extension value of the bond. Non-bonded monomers only experience a repulsive WCA potential. Regarding units, lengths are given in units of , which corresponds to the diameter of a monomer of unit mass , energy in units of and time in units of .

In this work, we build networks of monomers confined within a sphere of two different radii , namely and . The change in allows to vary the topology of the network, which becomes more compact for small and looser (and with more dangling ends) for larger  [15, 16]. The number of crosslinks is fixed to of the total number of monomers.

Implicit solvent. The implicit solvent is modeled through the addition of an attractive potential that acts between all monomers, either bonded or non-bonded, of the microgel [44, 45]:


where and  [44]. The potential is modulated by the parameter , which controls the solvophobicity of the monomers and plays the role of an effective temperature. For , monomers do not experience any attraction and good solvent conditions are reproduced while, by increasing up to 1.5, the monomers become fully attractive, mimicking bad solvent conditions. Previous analysis has shown that the volume phase transition takes place at  [15]. MD simulations are performed at constant reduced temperature (where is the Boltzmann constant) using the Nosè-Hoover thermostat and a timestep .

Adding an explicit solvent in MD simulations. We take a configuration of the microgel assembled as described above and perform MD simulations in the presence of a varying number of additional spheres that mimic the solvent particles which, for efficiency reasons, have also a diameter . The number of solvent particles varies between and in a simulation box of size , yielding solvent number densities , for which the LJ solvent is in the fluid regime. Lower densities would bring the LJ solvent to phase separate, while higher would lead to a crystallization of the solvent particles. All MD simulations with explicit solvent are performed with the LAMMPS simulation package [46] at making use of the Nosè-Hoover thermostat and a timestep . The center of mass of the microgel is fixed in the center of the simulation box. To model solvent-solvent interactions we use a Lennard-Jones (LJ) potential, . Here is the same as the one used in the WCA of monomer-monomer interactions.

The choice of the monomer-solvent (ms) interactions is crucial in order to implement the solvophobic effect, giving rise to the volume phase transition of the microgel. In this respect, we test different approaches. Our first choice is to employ again the LJ potential in which its depth is varied, so that the attractive contribution can be tuned. A weaker attraction would thus give rise to a more repulsive monomer-solvent interaction that should cause the shrinking of the microgel. However, as explained by analyzing the swelling curves in the Results section, a decrease of this parameter causes the repulsive barrier to be less efficient, giving rise to unphysical consequences on the swelling behavior. To fix this problem, we consider a -dependent Lennard-Jones potential [47], , defined as:


where is the same as for the monomer-monomer (Eq. (2)) and LJ solvent-solvent interactions, while plays the role of an inverse temperature (analogue to the inverse of in the implicit solvent model). Indeed, for large values of there is an attractive contribution between a monomer and a solvent particle, mimicking good solvent conditions, while for , the WCA potential is recovered and monomer-solvent interactions are purely repulsive. The potential is truncated and shifted at . The advantage of using such a potential with respect to the simple LJ interactions is that it allows to alter the monomer-solvent interactions, and thus the “quality” of the solvent, without affecting the excluded-volume part; this remains encoded in the term and it does not depend on .

Adding an explicit solvent in DPD simulations. DPD is a mesoscale simulation technique[26, 48] that treats solvent particles as coarse-grained beads and is able to describe hydrodynamic interactions through a momentum-conserving thermostat. In DPD simulations, particles and interact by three pairwise additive forces: a conservative force , a dissipative force and a random force where


Here with the position of particle , , , with the velocity of particle , and are weight functions, is a Gaussian random number with zero mean and unit variance and is the friction coefficient (here ); to ensure that the correct Boltzmann distribution is achieved at equilibrium, and . The interaction region for the dissipative force is defined in the same way as for the conservative force, i.e. . We refer the reader to Ref. [26] for more details on the DPD simulation technique.

Although the DPD protocol is often applied to all the species that are present in the simulation, here we make use of DPD only for the solvent-solvent and the solvent-monomer interactions, while keeping the microgel model unaltered. This hybrid technique allows to simulate a “fast” solvent by retaining important features of the microgel particle, such as the topology and the excluded volume interactions among monomers, already investigated in the implicit solvent case [15]. In DPD simulations, we fix the solvent number density at , using particles in a simulation box of side , and we tune the interaction parameters and the radius of the solvent beads until the swelling curve of the implicit solvent model is reproduced, in this case at . The same curve may be found by using different combination of these parameters in the limit in which the size of the solvent bead is comparable with that of the microgel monomer. All DPD simulations are also performed with LAMMPS at using the velocity-Verlet algorithm to integrate the equations of motion; the DPD thermostat is applied to the solvent particles only. Moreover, the center of mass of the microgel is fixed in the center of the simulation box as in MD simulations. The monomer-solvent interaction parameters , that for simplicity we call , plays the role of an effective temperature and, depending on its value, controls the volume phase transition of the polymer network. The repulsion coefficient for solvent-solvent interactions is fixed at .

Rescaling of swelling curves. The swelling degree of the microgel is expressed via the ratio between the absolute value of the gyration radius and its maximum value, obtained in good solvent conditions. The gyration radius is computed as , where indicates the position of the microgel center of mass. Since each model depends on a different “swelling parameter”, that we call generically , we scale all curves onto the implicit model one, using as the reference swelling parameter. For those explicit solvent models where a small value of the swelling parameter corresponds to a collapsed state of the microgel, i.e. and , the scale has to be inverted. In order to properly rescale the axes onto each other for two curves and , we consider two points on the first ( and ) and on the second curve ( and ), respectively. The rescaled -coordinate is calculated using the following relationship: , where and with .

Form factors. The microgel form factor is calculated as where the angular brackets indicate an average over different equilibrium configurations of the same microgel and over different orientations of the wavevector .

Cluster analysis in the kinetics of deswelling. To investigate the structural changes during the transient kinetics of deswelling, we define clusters within the microgel that are formed by non-bonded monomers only: two such monomers belong to the same cluster when their distance is smaller than , which roughly corresponds to the first peak of the radial distribution function. The cluster size distribution of clusters of size is calculated by averaging over six independent configurations of the microgel.


We acknowledge support from the European Research Council (ERC Consolidator Grant 681597, MIMIC).

Author contributions statement

FC, NG, LR and EZ performed simulations, analyzed results and wrote the paper.

Additional information

The authors declare no competing interests.


  • [1] Yunker, P. J. et al. Physics in ordered and disordered colloidal matter composed of poly (n-isopropylacrylamide) microgel particles. \JournalTitleRep. Progress in Physics 77, 056601 (2014).
  • [2] Lyon, L. A. & Fernandez-Nieves, A. The polymer/colloid duality of microgel suspensions. \JournalTitleAnnual Review of Physical Chemistry 63, 25–43 (2012).
  • [3] Fernandez-Nieves, A., Wyss, H., Mattsson, J. & Weitz, D. A. Microgel suspensions: fundamentals and applications (John Wiley & Sons, 2011).
  • [4] Pelton, R. Temperature-sensitive aqueous microgels. \JournalTitleAdvances in Coll. and Interf. Science 85, 1–33 (2000).
  • [5] Seiffert, S. Microgel capsules tailored by droplet-based microfluidics. \JournalTitleChemPhysChem 14, 295–304 (2013).
  • [6] Tavagnacco, L., Zaccarelli, E. & Chiessi, E. On the molecular origin of the cooperative coil-to-globule transition of poly(n-isopropylacrylamide) in water. \JournalTitlePhys. Chem. Chem. Phys. (2018). DOI 10.1039/C8CP00537K.
  • [7] Saunders, B. R. & Vincent, B. Microgel particles as model colloids: theory, properties and applications. \JournalTitleAdvances in Coll. and Interf. Science 80, 1–25 (1999).
  • [8] Jha, P. K., Zwanikken, J. W., Detcheverry, F. A., De Pablo, J. J. & De La Cruz, M. O. Study of volume phase transitions in polymeric nanogels by theoretically informed coarse-grained simulations. \JournalTitleSoft Matter 7, 5965–5975 (2011).
  • [9] Ghavami, A., Kobayashi, H. & Winkler, R. G. Internal dynamics of microgels: A mesoscale hydrodynamic simulation study. \JournalTitleJournal of Chemical Phys. 145, 244902 (2016).
  • [10] Kobayashi, H., Halver, R., Sutmann, G. & Winkler, R. G. Polymer conformations in ionic microgels in the presence of salt: Theoretical and mesoscale simulation results. \JournalTitlePolymers 9, 15 (2017).
  • [11] Ahualli, S., Martín-Molina, A., Maroto-Centeno, J. A. & Quesada-Pérez, M. Interaction between ideal neutral nanogels: A monte carlo simulation study. \JournalTitleMacromolecules (5), 2229–2238 (2017).
  • [12] Nikolov, S., Fernandez-Nieves, A. & Alexeev, A. Mesoscale modeling of microgel mechanics and kinetics through the swelling transition. \JournalTitleApplied Mathematics and Mechanics 39, 47–62 (2018).
  • [13] Ghavami, A. & Winkler, R. G. Solvent induced inversion of core–shell microgels. \JournalTitleACS Macro Letters 6, 721–725 (2017).
  • [14] Kobayashi, H. & Winkler, R. G. Structure of microgels with debye–hückel interactions. \JournalTitlePolymers 6, 1602–1617 (2014).
  • [15] Gnan, N., Rovigatti, L., Bergman, M. & Zaccarelli, E. In silico synthesis of microgel particles. \JournalTitleMacromolecules 50, 8777–8786 (2017).
  • [16] Rovigatti, L., Gnan, N. & Zaccarelli, E. Internal structure and swelling behaviour of in silico microgel particles. \JournalTitleJ. Phys.: Cond. Matter 30, 8pp (2018).
  • [17] Grest, G. S. & Kremer, K. Molecular dynamics simulation for polymers in the presence of a heat bath. \JournalTitlePhysical Review A 33, 3628 (1986).
  • [18] Pham, T. T., Schiller, U. D., Prakash, J. R. & Dünweg, B. Implicit and explicit solvent models for the simulation of a single polymer chain in solution: Lattice boltzmann versus brownian dynamics. \JournalTitleJournal of Chemical Physics 131, 164114 (2009).
  • [19] Spaeth, J. R., Kevrekidis, I. G. & Panagiotopoulos, A. Z. A comparison of implicit-and explicit-solvent simulations of self-assembly in block copolymer and solute systems. \JournalTitleJournal of Chemical Physics 134, 164902 (2011).
  • [20] Neyt, J.-C., Wender, A., Lachet, V., Ghoufi, A. & Malfreyt, P. Quantitative predictions of the interfacial tensions of liquid–liquid interfaces through atomistic and coarse grained models. \JournalTitleJournal of Chemical Theory and Computation 10, 1887–1899 (2014).
  • [21] Rumyantsev, A. M., Gumerov, R. A. & Potemkin, I. I. A polymer microgel at a liquid–liquid interface: theory vs. computer simulations. \JournalTitleSoft Matter 12, 6799–6811 (2016).
  • [22] Isa, L., Buttinoni, I., Fernandez-Rodriguez, M. & Vasudevan, S. Two-dimensional assemblies of soft repulsive colloids confined at fluid interfaces. \JournalTitleEPL 119, 26001 (2017).
  • [23] Scheidegger, L. et al. Compression and deposition of microgel monolayers from fluid interfaces: particle size effects on interface microstructure and nanolithography. \JournalTitlePhysical Chemistry Chemical Physics 19, 8671–8680 (2017).
  • [24] Brugger, B., Vermant, J. & Richtering, W. Interfacial layers of stimuli-responsive poly-(n-isopropylacrylamide-co-methacrylicacid)(pnipam-co-maa) microgels characterized by interfacial rheology and compression isotherms. \JournalTitlePhysical Chemistry Chemical Physics 12, 14573–14578 (2010).
  • [25] Schwenke, K., Isa, L., Cheung, D. L. & Del Gado, E. Conformations and effective interactions of polymer-coated nanoparticles at liquid interfaces. \JournalTitleLangmuir 30, 12578–12586 (2014).
  • [26] Groot, R. D. & Warren, P. B. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. \JournalTitleJournal of Chemical Physics 107, 4423–4435 (1997).
  • [27] Conley, G. M., Aebischer, P., Nöjd, S., Schurtenberger, P. & Scheffold, F. Jamming and overpacking fuzzy microgels: Deformation, interpenetration, and compression. \JournalTitleScience Advances 3, e1700969 (2017).
  • [28] Bergmann, S., Wrede, O., Huser, T. & Hellweg, T. Super-resolution optical microscopy resolves network morphology of smart colloidal microgels. \JournalTitlePhysical Chemistry Chemical Physics 20, 5074–5083 (2018).
  • [29] Stieger, M., Pedersen, J. S., Lindner, P. & Richtering, W. Are thermoresponsive microgels model systems for concentrated colloidal suspensions? a rheology and small-angle neutron scattering study. \JournalTitleLangmuir 20, 7283–7292 (2004).
  • [30] Boon, N. & Schurtenberger, P. Swelling of micro-hydrogels with a crosslinker gradient. \JournalTitlePhysical Chemistry Chemical Physics 19, 23740–23746 (2017).
  • [31] Bischofberger, I. & Trappe, V. New aspects in the phase behaviour of poly-n-isopropyl acrylamide: systematic temperature dependent shrinking of pnipam assemblies well beyond the lcst. \JournalTitleSci. Rep. 5, 15520 (2015).
  • [32] Reddy, G. & Yethiraj, A. Implicit and explicit solvent models for the simulation of dilute polymer solutions. \JournalTitleMacromolecules 39, 8536–8542 (2006).
  • [33] Chang, R. & Yethiraj, A. Solvent effects on the collapse dynamics of polymers. \JournalTitleJournal of Chemical Physics 114, 7688–7699 (2001).
  • [34] Pham, T. T., Bajaj, M. & Prakash, J. R. Brownian dynamics simulation of polymer collapse in a poor solvent: influence of implicit hydrodynamic interactions. \JournalTitleSoft Matter 4, 1196–1207 (2008).
  • [35] Nikolov, S., Fernandez-Nieves, A. & Alexeev, A. Mesoscale modeling of microgel mechanics and kinetics through the swelling transition. \JournalTitleApplied Mathematics and Mechanics 39, 47–62 (2018).
  • [36] Seiffert, S. Impact of polymer network inhomogeneities on the volume phase transition of thermoresponsive microgels. \JournalTitleMacromolecular Rapid Communications 33, 1135–1142 (2012).
  • [37] Habicht, A., Schmolke, W., Lange, F., Saalwächter, K. & Seiffert, S. The non-effect of polymer-network inhomogeneities in microgel volume phase transitions: Support for the mean-field perspective. \JournalTitleMacromolecular Chemistry and Physics 215, 1116–1133 (2014).
  • [38] Schmidt, S. et al. Influence of microgel architecture and oil polarity on stabilization of emulsions by stimuli-sensitive core–shell poly (n-isopropylacrylamide-co-methacrylic acid) microgels: Mickering versus pickering behavior? \JournalTitleLangmuir 27, 9801–9806 (2011).
  • [39] Geisel, K., Isa, L. & Richtering, W. Unraveling the 3d localization and deformation of responsive microgels at oil/water interfaces: a step forward in understanding soft emulsion stabilizers. \JournalTitleLangmuir 28, 15770–15776 (2012).
  • [40] Li, H., Yan, G., Wu, S., Wang, Z. & Lam, K. Numerical simulation of controlled nifedipine release from chitosan microgels. \JournalTitleJournal of Applied Polymer Science 93, 1928–1937 (2004).
  • [41] Bysell, H., Månsson, R., Hansson, P. & Malmsten, M. Microgels and microcapsules in peptide and protein drug delivery. \JournalTitleAdvanced Drug Delivery Reviews 63, 1172–1185 (2011).
  • [42] Weeks, J. D., Chandler, D. & Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. \JournalTitleJournal of Chemical Physics 54, 5237–5247 (1971).
  • [43] Bernabei, M., Moreno, A. J., Zaccarelli, E., Sciortino, F. & Colmenero, J. Chain dynamics in nonentangled polymer melts: A first-principle approach for the role of intramolecular barriers. \JournalTitleSoft Matter 7, 1364–1368 (2011).
  • [44] Soddemann, T., Dünweg, B. & Kremer, K. A generic computer model for amphiphilic systems. \JournalTitleEuropean Physical Journal E 6, 409–419 (2001).
  • [45] Verso, F. L., Pomposo, J. A., Colmenero, J. & Moreno, A. J. Simulation guided design of globular single-chain nanoparticles by tuning the solvent quality. \JournalTitleSoft Matter 11, 1369–1375 (2015).
  • [46] Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. \JournalTitleJournal of Computational Physics 117, 1–19 (1995). URL
  • [47] Rovigatti, L., Capone, B. & Likos, C. N. Soft self-assembled nanoparticles with temperature-dependent properties. \JournalTitleNanoscale 8, 3288–3295 (2016).
  • [48] Keaveny, E. E., Pivkin, I. V., Maxey, M. & Em Karniadakis, G. A comparative study between dissipative particle dynamics and molecular dynamics for simple-and complex-geometry flows. \JournalTitleJ. Chem. Phys. 123, 104107 (2005).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description