Modelling realistic microgels in an explicit solvent

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.


Introduction
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 10 − 1000 nm by standard synthesis methods 4 , and even reach much larger scale (up to 100 µm) with microfluidic techniques 5 . These particles undergo a Volume Phase Transition (VPT) in water at a temperature of ≈ 305K, 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 [12][13][14] .
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] . 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 V α , 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 Z = 25σ ) 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.

Results
Similarly to solvent-solvent interactions, a straightforward choice for the monomer-solvent ones is the LJ potential 25 where, by varying the energy minimum ε ms , 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 R g is shown as a function of ε ms : 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 R g is observed for ε ms → 0: 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 V λ 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 a ms in Eq. (4): for small values of a ms the microgel is swollen, while it contracts when a ms increases. We notice that R g 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.

3/13
In order to establish a correspondence between different models, we rescale the explicit solvent data onto the implicit one, V α where α is the solvophobic parameter (see Methods). Figure 2(a) shows the normalized R g /R max g , where R max g is the value of the gyration radius at maximum swelling, as a function of the effective swelling parameter χ eff . 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 χ eff , namely the DPD and MD V λ models. Of the latter, we consider only the case with the highest solvent density, ρ = 0.87, 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 V λ 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 V λ 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, Z = 15σ . 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 χ eff 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
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 χ eff 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 χ eff , 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 P(q) can be measured by small angle neutron or x-ray scattering experiments. Thus, in contrast to density profiles, numerical P(q) 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 15,30 . We thus directly evaluate the form factors of the microgel across the VPT and present them in Fig. 4 as a function of wavevector q for the same values of swelling parameters used in Fig. 3. 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 χ eff increases and the solvent quality decreases, P(q) 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 χ eff , the MD results are found to be always shifted to a slightly smaller q-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 (qσ 1) the MD form factor systematically overestimates the DPD and implicit-solvent ones for intermediate and large values of χ eff . However, all curves superimpose again at qσ ∼ 7, 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 ρ s inside the microgel. The normalized profile ρ s /ρ s,bulk , where ρ s,bulk 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

5/13
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 χ eff ∼ 0.6, 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 90% 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  normalized with respect to the bulk solvent density ρ s,bulk , as a function of the distance r 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
We now repeat the above structural analysis for a more compact microgel topology obtained with a smaller confining radius (Z = 15σ ), 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 χ eff 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 χ eff . 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 χ eff . 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 qσ ∼ 3.0 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 ∼ 2σ apart, i.e. with a solvent particle in between them. Such a feature is smeared out at increasing χ eff , 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 χ eff 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
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 χ eff correspond to an almost fully collapsed state (χ eff ∼ 1.1) and to a state close to the VPT (χ eff ∼ 0.7). 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 χ eff . In all cases the curves reach at long time the same value of R g 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 R g 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 .

Discussion
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 χ eff 1.4, amounting to temperatures 60 • 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 12,33,34 . 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 [21][22][23][24]39 . 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 .

Methods
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 with k F = 15 an adimensional spring constant and R 0 = 1.5 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 m, energy in units of ε and time in units of mσ 2 /ε. In this work, we build networks of N ∼ 5000 monomers confined within a sphere of two different radii Z, namely Z = 15 σ and Z = 25 σ . The change in Z allows to vary the topology of the network, which becomes more compact for small Z and looser (and with more dangling ends) for larger Z 15,16 . The number of crosslinks is fixed to 3.2% of the total number of monomers.
Implicit solvent. The implicit solvent is modeled through the addition of an attractive potential V α that acts between all monomers, either bonded or non-bonded, of the microgel 44,45 : where δ = π(2.25 − 2 1/3 ) −1 and β = 2π − 2.25δ 44 . The potential is modulated by the parameter α, which controls the solvophobicity of the monomers and plays the role of an effective temperature. For α = 0, 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 α ∼ 0.6 15 . MD simulations are performed at constant reduced temperature T * = k B T /ε = 1 (where k B is the Boltzmann constant) using the Nosè-Hoover thermostat and a timestep t * = 0.002.
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 2.5 × 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 ε ms 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 , V λ (r), 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 λ = 0, the WCA potential is recovered and monomer-solvent interactions are purely repulsive. The potential is truncated and shifted at 2.5σ . 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 V WCA 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 i and j interact by three pairwise additive forces: a conservative force F C i j , a dissipative force F D i j and a random force F R i j where Here r i j = r i − r j with r i the position of particle i, r i j = |r i j |,r i j = r i j /r i j , v i j = v i − v j with v i the velocity of particle i, w D (r i j ) and w R (r i j ) are weight functions, θ is a Gaussian random number with zero mean and unit variance and γ is the friction coefficient (here γ = 4.0); to ensure that the correct Boltzmann distribution is achieved at equilibrium, w D (r i j ) = [w R (r i j )] 2 and σ 2 R = 2γk B T . The interaction region for the dissipative force is defined in the same way as for the conservative force, i.e. w D (r i j ) = 1 − r i j /R c . 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 ρ = 0.73, using N = 2.5 × 10 5 particles in a simulation box of side 70σ , 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 r c = 1.75σ . 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 T * = 1.0 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 a ms i j , that for simplicity we call a ms , 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 a ss i j σ = a ss σ = 25ε.

10/13
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 R g = [N −1 ∑ N i (r i − r CM ) 2 ] 1/2 , where r CM indicates the position of the microgel center of mass. Since each model depends on a different "swelling parameter", that we call generically χ eff , 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. V LJ and V λ , the scale has to be inverted. In order to properly rescale the x axes onto each other for two curves A and B, we consider two points on the first (x A 1 and x A 2 ) and on the second curve (x B 1 and x B 2 ), respectively. The rescaled x-coordinate is calculated using the following relationship: x new = x − x A ∆x B /∆x A + x B , where x i = 0.5(x i 1 + x i 2 ) and ∆x i = x i 1 − x i 2 with i = A, B. Form factors. The microgel form factor P(q) is calculated as P(q) = 1 N ∑ i j exp (−iq · r i j ) where the angular brackets indicate an average over different equilibrium configurations of the same microgel and over different orientations of the wavevector q. 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 1.2σ , which roughly corresponds to the first peak of the radial distribution function. The cluster size distribution n(s) of clusters of size s is calculated by averaging over six independent configurations of the microgel.