1 Introduction

Hydrogels have become a highly promising platform for developing novel microscale and nanoscale technologies, especially in biomedical engineering, where they are used as advanced drug delivery vehicles, microscopic manipulators/sensors, and scaffolds to regrow bones, muscles, and organs^[1]. Recent experiments show that by varying the relaxation rates of the underlying hydrogel matrix, mesenchymal stem cells can be predisposed towards osteogenic differentiation (vs. adipogenic differentiation), which is important for bone tissue formation^[2]. Methacrylated gelatin (MG) scaffolds have been used to regrow cardiac muscle cells^[3]. It has been shown that by enriching the MG network with carbon nanotubes, the beating frequency of the cells could be varied. Hydrogels have also been used extensively to create novel drug delivery devices. A new hydrogel-based micro-device capable of aiding diabetics by delivering a sustained release of insulin in response to enzymatic oxidation of glucose has been recently demonstrated^[4]. Microscopic hydrogel particles have been used to encapsulate living cells to create three-dimensional (3D) micro-environments for cell culturing^[5] and for controlled release of DNA targets^[6]. Ultrasoft microgels have been shown to act as artificial platelets facilitating blood clotting^[7], whereas dense colloidal suspensions of such microgels enable assays to evaluate cell invasiveness^[8]. In addition, stimuli-responsive multi-layered gel structures are promising for microactuator applications^[9-11] and for designing self-propelling robots^[12-14].

Given the overwhelming potential of hydrogels, our aim is to develop a mesoscopic computational model, which can accurately capture the mechanical properties and swelling kinetics of microgel particles through the volume-phase transition, thereby facilitating the use of hydrogels in practical applications. The development of mesoscale models has attracted significant attention thanks to their ability to consider intermediate spatial and time scales that pose difficulties for more conventional atomistic and macroscopic (continuous) descriptions. Atomistic models provide a great wealth of detailed information regarding gel hydration and swelling^[15]. However, the applications in the relevant nm-μm length scales are typically not computationally accessible with atomistic simulations, hindering the potential utility in solving engineering problems. Continuous approaches enable access to much larger length and time scales compared with atomistic models. In this case, however, connecting individual chain parameters to bulk material properties is not an easy task for hydrogel materials^[16] due to the complex interplay among elastic, hydrodynamic, and chemical interactions. As a result, mesoscale models of polymeric materials are gaining popularity due to their ability to tackle problems inaccessible by other methods^[17-19].

In this work, we examine a mesoscale hydrogel model based on dissipative particle dynamics (DPD). We focus on neutral hydrogel networks, and probe how individual chain parameters affect bulk gel properties in solvents of different quality. In particular, we examine how the bulk modulus changes through the volume-phase transition of microgel particles, and compare the simulation results with predictions from Flory-Rehner's theory. Finally, we employ our mesoscale model to simulate the kinetics of spherical microgel particles during swelling and deswelling transitions and pinpoint important differences in microgel microstructure along these processes.

2 Methodology

We use DPD, a particle-based method where beads represent mesoscale molecular clusters and fluid volumes^[20-22], to model a mesoscale polymeric network immersed in a viscous solvent. The beads follow Newtonian dynamics and interact via soft potentials allowing a larger time step for the integration of the equation of motion. This in turn enables modeling systems over extended time scales, which are unachievable with atomistic methods. Furthermore, DPD uses pairwise interactions that preserve local hydrodynamics critical for modeling systems involving fluid flows. DPD has long been successfully tested against scaling theories for modeling simple systems^[23-30] and experimental data^[31-39], showing that it can capture relevant physical effects.

The governing dynamics between beads in DPD is set by three main forces

acting between a given bead i and its neighboring bead j located within a cutoff radius r_c. There is a conservative repulsive force

which leads to the bead excluded volume, a dissipative force

representing the effects of viscosity, and a random force

associated with thermal fluctuations. The latter two forces are related by

due to the fluctuation-dissipation theorem. In all forces, is a weighing function with

Additionally, , ξ_ij is a standard normal variable with zero mean, k_B is the Boltzmann constant, and T is the temperature.

To model the polymeric chains in the microgels, we use a bead-spring model (see Fig. 1(a)), where DPD beads are connected using a harmonic bond potential, a bending angle potential, and a segmented-repulsive potential (SRP). The energy for the bond potential is given by

where k_bond is the bond stiffness, and r_eq is the equilibrium separation length between beads. The energy for the angle potential is given by

where k_bend is the bending stiffness, and θ is the angle between two polymer bonds sharing a common bead. The segmented repulsive potential prevents polymer chains from crossing each other. It is implemented by imposing a soft repulsive force between chain bonds. The associated force is given by

where r_cSRP is the cutoff distance for the SRP potential, r_min is the minimum distance between two bonds, and C is the strength of the potential^[40].

In our simulations, we set

Here, ρ is the number density of DPD beads in the system, and Δt is the simulation time step. The repulsions between solvent beads and between beads within polymer chains is

To model swelling and deswelling of the polymer gel, we alter the polymer-solvent repulsion a_p-s in the range between 20 and 35. For the sake of simplicity in what follows, we use a without a subscript to indicate the polymer-solvent repulsion, i.e., . The bond stiffness for the polymer chains is set to k_bond = 35 with r_eq = 0.6. The bending stiffness is set to k_bend = 10, except for the single chain simulations where k_bend varies in the range between 0 and 10. Note that all dimensional parameters are expressed in DPD units, unless indicated otherwise.

To create the gel network, we randomly distribute 3 600 crosslink points in a cubic box^[38]. We then find all crosslink points separated by a distance in the range between 0.8r_eqN and 1.1r_eqN, and randomly connect them with connectivity up to 4 with straight polymer chains. Here, N is the average number of monomers per chain. The number of monomers in each chain is set such that the distance between beads is close to r_eq. A microgel with spherical shape is formed from the initial cubic gel by cutting out the chains outside the gel radius, with disconnected chain segments removed via a depth-first search algorithm^[41] (see Fig. 1(b)). The final average connectivity is found to be about 3.4, mainly due to free chains at the outer gel surface.

For all the microgel networks in our simulations, we use an identical initial crosslinking configuration. We investigate the effect of the number of beads per chain by setting N = 6, 9, 12, and 24. To create networks with different N, we rescale the relative positions of crosslink points such that the average number of the monomers in the chains connecting them is close to the desired value. Thus, all networks share the same characteristics including the number of crosslinks, connectivity, and entanglements, and only differ by the number of monomers per chain. This corresponds to changing the crosslink-to-monomer ratio, which is often done in experiments^[42-43]. We conduct our simulations in a periodic computational domain. In the simulations with microgel networks, we choose a domain size in the range from 100×100×100 to 170×170×170 such that any interactions of the microgel with its periodical images are prevented. The domain size for single chain simulations is 20×20×20. At the beginning of the simulations, we equilibrate the microgel network in solvent until a steady state situation is reached.

3 Results and discussion 3.1 Single chain model

Neutral poly (N-isopropylacrylamide) (PNIPAM) gels are composed of polymer chains that exhibit flexible behavior^[44]. We use this property as a guiding criterion for constructing our mesoscale model for a polymer chain. Flory's theory allows us to characterize the flexibility of individual polymer chains in a dilute solution in terms of the relation between the end-to-end distance R_end and the number of monomers N. This relation is given by

where l_k is the Kuhn length, and v is Flory's exponent equal to 3/5 for flexible, self-avoiding chains in good solvent^[45].

When the bead-spring model is used without bending stiffness, we reproduce the flexible behavior of a polymer chain with just a few beads per chain. In Fig. 2(a), we show a log-log plot of R_end as a function of the number of beads for chains with different k_bend and the chain solvent repulsion a=20. The lines in the figure correspond to v =3/5. When k_bend = 0, the bead-spring model conforms to the expected scaling behavior, even for very short chains. Including bending stiffness, however, leads to deviations for chains with a small number of beads (see k_bend = 5 and k_bend = 10 in Fig. 2(a)). These chains behave as flexible chains only if the number of beads in the chain is sufficiently large. In Fig. 2(b), we show the dependence of the minimum number of beads in the chain N_min for the chain to be flexible as a function of the bending stiffness. This number increases as k_bend increases. Hence, chains with N > N_min behave as flexible chains, whereas chains with N < N_min are semi-flexible. Note that in our mesoscale microgel model, we are interested in using shorter chains to reduce the size of the system and accelerate the simulations.

In addition to changing the values of N_min, bending stiffness also affects the magnitude of the Kuhn length, which represents the effective length of orientationally independent segments within a polymer chain. For k_bend = 0, the Kuhn length is relatively short and is equal to about 0.6, which is the equilibrium length r_eq for stretching springs connecting the beads in the chain. The addition of bending stiffness leads to an increased Kuhn length which increases nearly linearly with k_bend within the tested range (see Fig. 3). As we show below, larger values of Kuhn length are critical for modeling highly porous swollen gels.

3.2 Microgel swelling

Hydrogels are characterized by a high porosity in good solvents^{[44, 46]}, which plays a fundamental role in gel transport properties, mechanics, and swelling kinetics. It is, therefore, important that our mesoscale microgel model can exhibit the relevant values of porosity in good solvent conditions.

In Fig. 4(a), we plot the gel porosity P as a function of k_bend for chains with different N. To determine the microgel porosity, we first measure the microgel volume V by constructing a surface mesh enclosing the microgel network (see Fig. 4(b))^[47]. The surface mesh is defined using a spherical probe with a radius of about half of the microgel mesh size. We then calculate microgel porosity based on the number of solvent beads located within the surface mesh. We find that P increases with chain length. However, when k_bend = 0, the porosity does not exceed 0.87 even for relatively long chains with N = 24. A further increase in the chain length can lead to higher porosity. However, this also significantly extends the computation time, as the domain size scales as N³.

Long computational times can be avoided by noting that P can be increased by using shorter chains with a greater k_bend and, thus, a greater Kuhn length. Indeed, we find that when k_bend = 10, the porosity is about 0.95 for the chains with N = 24, which is close to typical experimental values for hydrogels. Furthermore, these chains behave as flexible chains, which is also typical in experiments. We, therefore, use the chains with k_bend = 10 in our mesoscale model. Using chains with higher bending stiffness will further increase microgel porosity. However, this will require longer chains in order to ensure that fully flexible behavior is maintained, which in turn will lead to excessive computational time.

We conduct simulations, in which we place a spherical microgel network in the solutions with different solvency that is varied by changing the network solvent repulsion coefficient a. Larger values of a correspond to bad solvents, whereas lower values of a correspond to good solvents. At each a, we calculate the network porosity (see Fig. 5(a)) for gels with different N. In bad solvents, all microgel networks have a nearly zero porosity, indicating that the microgel expels practically all of the solvent and adopts a nearly collapsed state. Note, that P>0, which agrees with what is found in experiments^[48]. In good solvents, the microgel swells, resulting in values of P approaching 0.95 for large N.

The microgel volume change, as the solvent changes from good to bad, strongly depends on N. For N = 6, the volume increases approximately 7 times, when the microgel changes from the collapsed state to the swollen state. For N = 24, the volume increases about 30 times between collapsed and swollen states. We find that, as the chain length increases, the volume-phase transition becomes sharper, as indicated by the greater values of the porosity gradient, |ΔP/Δa|, for the microgels with longer chains in Fig. 5(b). This behavior is related to the microgel softening due to the decreasing crosslink density ρ_c =η/V_a in microgels with larger N. Here, η is the total number of crosslinks, and V_a is the gel volume for a given value of a. This result is in qualitative agreement with expectations from Flory-Rehner's theory of polymer gels^[49].

Larger volume changes can be also obtained by using polymer chains with greater bending stiffness. In good solvent, stiffer chains lead to a higher microgel porosity, whereas in bad solvent conditions, the size of the collapsed gels is nearly insensitive to k_bend and is mostly defined by the excluded volume of the chains.

3.3 Bulk modulus

We measure the bulk modulus of our mesoscale gels by quasi-statically varying the total osmotic pressure imposed on the microgel network. To change the osmotic pressure, we create a semi-permeable spherical shell around the microgel (see Fig. 6). The spherical shell interacts with the gel particles via a harmonic potential, while the solvent particles are allowed to pass freely through the shell. Varying the radius of the spherical shell, thus, effectively mimics a change in the osmotic pressure. For each solvency a, we vary the radius of the semi-permeable shell, and evaluate the corresponding network volume and pressure. To find the microgel volume V, we construct a surface mesh around the gel. The network pressure is calculated as

where W is the number of microgel beads in the network, and σ_ii^ave (i=x, y, z) are the average normal stresses per bead in the x-, y-, and z-directions. The stress is averaged over all beads in the microgel network.

Figure 7(b) shows the resultant bulk moduli K for microgels with different N through the volume phase transition. The pressure is calculated for different values of solvency a in the linear region of microgel compression, where Π_network is a linear function of V (see the inset in Fig. 7(a)). The variation of Π_network with V in this linear region is used to evaluate the microgel bulk modulus as follows:

This procedure is repeated for different a, covering the entire range of the microgel swelling transition.

Figure 7(b) shows the resultant bulk modulus K for the microgels with different N through the volume phase transition. We find that the bulk modulus monotonically increases as the solvency decreases. This increase in K is related to a decrease in the microgel porosity and an increase in the crosslink density ρ_c as a increases. In a good solvent with a < 25 (swollen state), the bulk modulus decreases with N. In this case, increasing N decreases the crosslink density, resulting in gel softening. This softening of the microgels with increasing N also leads to a sharper volume-phase transition (see Fig. 5).

Analytically, the microgel bulk modulus can be evaluated with the Flory-Rehner theory, which postulates that, for nonionic gels, the total osmotic pressure Π_network is composed of a mixing osmotic pressure Π_mix due to polymer-solvent mixing and an elastic osmotic pressure Π_e due to network elasticity^[49]. In equilibrium with pure solvent, Π_mix balances Π_e. Using that

and Flory's descriptions of Π_mix and Π_e, we can obtain an expression for the microgel bulk modulus as follows:

where ϕ = 1 -P is the polymer volume fraction, N_c/V_o is the number of polymer chains in the preparation state per volume, which corresponds to the nearly collapsed state in our simulations, v_s is the molar volume of the solvent, ϕ_o is the volume fraction in the preparation state, and χ is the solvency parameter. For our microgels, we evaluate the solvent molar volume as

where r_c is the cutoff radius of the DPD potential, N_A is Avogadro's number, and ρ is the density in our simulations^[50].

We fit our simulation results for K versus a with Flory-Rehner's theory, using the relation between Π_network, and thus ϕ, with a. We consider the solvency parameter as a function of the polymer volume fraction

and set k_BT/v_s, N_cv_s/V_o, χ₀, χ₁, and χ₂ as the fitting parameters. In our simulations, the number of polymer chains per volume in the preparation state is N_c/V_o ~ 0.25 chains per unit volume. The obtained values of the fitting parameters are summarized in Table 1. We find that the values are nearly independent of the chain length N.

The resulting fits to our bulk modulus data correctly describe our simulation results (see Fig. 7(b)). We note that if the fits are performed under the assumption that χ is linearly dependent on ϕ, good agreement between our simulations and the theory cannot be obtained. Thus, the quadratic dependence of χ on ϕ is essential. This is consistent with previous experiments, and suggests that many-body interactions become important near the swelling transition^[49].

Using the results for χ from the fits, we can establish how this parameter depends on the network-solvent repulsion coefficient a. This is shown for different N in Fig. 8. We find that the χ-a dependence is not sensitive to the value of N. Furthermore, in good solvents, the solvency parameter is about 0.5, whereas in bad solvents, it increases to about 1. This behavior agrees reasonably well with the experimental data^[51-53].

3.4 Swelling kinetics

When the gel undergoes a volume transition, the kinetics of this process is defined by a balance among osmotic pressure, gel elasticity, and viscous drag due to the solvent penetrating in or out of the gel network. The swelling kinetics of a spherical gel can be analyzed analytically^[54]. The change in the gel radius is given by

where Δr_o is the total increase in the radius defined by

In the above equation, r_final and r_ini are the final and initial equilibrated gel radii, respectively, τ_s is the swelling time constant, and t is time. The swelling time constant is given by

In the above equation, D is the collective diffusion of the gel characterizing its elastic relaxation in a viscous solvent, and is expressed by

where K is the bulk modulus, μ is the shear modulus, and f is the friction coefficient between the polymer network and the fluid.

To examine whether our mesoscale gel model can properly describe swelling kinetics, we simulate the swelling of a spherical microgel that is initially in the collapsed state (a = 35). At t = 0, the collapsed gel is placed in a good solvent (a = 20). In our simulations, we track how the microgel radius changes as the microgel undergoes swelling in good solvent. The results of the simulations for different chain lengths N are shown in Fig. 9(a). We fit the simulation data to Tanaka's theory to determine the relaxation constant τ_s and the corresponding value of the collective diffusion coefficient D, which are given in Table 2. The data indicates that τ_s increases roughly proportional to N², whereas the collective diffusion D decreases with the increase in N^[55]. As a result, networks with longer chains swell slower than the networks with smaller chains. Figure 9(a) shows close agreement between the theoretical results and the simulation results for all N, indicating that our microgel model can properly capture not only microgel mechanics but also microgel interactions with the viscous solvent.

We also examine the microgel kinetics during deswelling. In this scenario, the swollen spherical gel (a = 20) is introduced to a bad solvent (a = 35). The time evolution of the microgel size is shown in Fig. 9(b). We find that the deswelling simulation data can also be fitted to Tanaka's theory, although the agreement is somewhat less accurate when compared with the swelling process. Comparing the fitting parameters in Table 2, we find that the relaxation time is considerably shorter for deswelling microgels, especially for longer chain lengths. This can be expected since the relaxation time is proportional to r_final², which is much smaller for collapsed microgels than for swollen microgels. We also find that the collective diffusion D is about 3 times smaller for deswelling microgels, indicating differences in the internal structure of the microgels during the swelling and deswelling transitions.

Figure 10 illustrates the microgels that undergo swelling (top row) and deswelling (bottom row) transitions. During swelling, the particle gradually increases in size while maintaining a nearly homogeneous structure (see Fig. 10(a)). In contrast, during deswelling, the microgel exhibits a remarkably different internal structure at intermediate times. Our simulations indicate that the deswelling from a fully swollen state proceeds through network coarsening that manifests as formation of chain bundles at the outskirts of the microgels (see Fig. 10(b)). Similar coarsening behavior has been reported for isochore phase separation in hydrogels^[56]. Eventually, however, the microgel volume decreases and the solvent is expelled, allowing the microgel to shrink to its collapsed size. Note that this process is different from the rapid deswelling of macroscopic gel, whereby a skin of collapsed gels is formed on the outer gel surface. This dense skin prevents solvent transfer from the hydrogel interior, thereby delaying the gel deswelling^[57].

4 Conclusions

We have examined the mechanics and swelling kinetics of neutral microgel particles using a mesoscale polymer network model. The model network is immersed in an explicit viscous solvent modeled using dissipative particle dynamics. We alter the solvency by changing the repulsion between the network and the solvent. To prevent chains from crossing each other, we employ a bond-bond repulsive potential. To model sufficiently large microgel networks, we identify single chain parameters leading to experimentally realistic microgel properties. We find that, to model swollen microgels with high porosity, bending stiffness needs to be incorporated into the polymer chains to increase their Kuhn length. This imposes a limitation on the minimum polymer chain length when flexible chain mechanics is required.

We show that microgels undergo a continuous volume transition when the solvency is changed. In bad solvents, the networks contract to a nearly collapsed state, whereas in good solvents, the microgels swell with a porosity that increases proportionally to the chain length. Microgels with higher porosity are characterized by a sharper volume transition. To probe the effect of the volume-phase transition on the microgel mechanics, we employ our simulation model to obtain the microgel bulk modulus. We find that the bulk modulus changes monotonically through the volume transition. In the collapsed state, the bulk modulus is independent of the chain length, and is defined by the chain-chain repulsion. In swollen networks, the modulus decreases with increasing porosity and decreasing crosslink density. We fit our bulk modulus data to the Flory-Rehner theory, and show that our mesoscale model is in good agreement with the theory. From the fitting, we obtain the Flory-Huggins solvency parameter. This procedure allows establishing a relation between the solvency parameter and the solvent-network repulsion parameter. We find that the solvency parameter exhibits a quadratic dependence on the polymer volume fraction, which is consistent with the experimental results, reflecting the importance of many-body interactions through the swelling transition.

Finally, we study the swelling kinetics of our microgel particles. Our simulations show close agreement to Tanaka's theory for spherical gels. This result indicates that the mesoscale model can properly capture the hydrodynamics of network-solvent interactions. Importantly, we find that, during swelling, microgels preserve a nearly homogeneous internal structure, whereas rapid deswelling leads to the formation of chain bundles in the outskirts of the particle and network coarsening. The difference in internal microgel structure is manifested by different kinetic constants, characterizing the speed of the swelling and deswelling transitions. In its current implementation, our mesoscale model is limited to modeling microgels exhibiting a continuous volume phase transition. Thus, our microgel model cannot simulate gels with a negative Poisson ratio and discontinuity in the bulk modulus, which are observed near the critical volume phase transition for certain microgels^[58]. More elaborate chain-chain and chain-solvent interactions need to be integrated into the model to capture these effects, which is an important direction for future model development.