Using Patchy Particles to Prevent Local Rearrangements in Models of Non-equilibrium Colloidal Gels.

Simple models based on isotropic interparticle attractions often fail to capture experimentally observed structures of colloidal gels formed through spinodal decomposition and subsequent arrest: the resulting gels are typically denser and less branched than their experimental counterparts. Here, we simulate gels formed from soft particles with directional attractions ("patchy particles"), designed to inhibit lateral particle rearrangement after aggregation. We directly compare simulated structures with experimental colloidal gels made using soft attractive microgel particles, by employing a "skeletonization" method that reconstructs the three-dimensional backbone from experiment or simulation. We show that including directional attractions with sufficient valency leads to strongly branched structures compared to isotropic models. Furthermore, combining isotropic and directional attractions provides additional control over aggregation kinetics and gel structure. Our results show that the inhibition of lateral particle rearrangements strongly affects the gel topology and is an important effect to consider in computational models of colloidal gels.

Introduction

The aggregation of suspended particulate matter due to attractive interactions can cause the formation of a volume-spanning network,[1] which has been observed in various materials such as foodstuffs,[2] proteins,[3] and cement.[4] The process of gelation relies on spinodal decomposition into particle-poor and particle-rich phases, where the latter arrests to form a percolating network whose precise nature depends on parameters such as the particle volume fraction ϕ, the attraction strength, and the quench speed.[5−10] This characteristic spinodal decomposition and subsequent arrest behavior can be recreated with computational models involving strong, isotropic attractions between spherical particles.[11] However, using such models often leads to discrepancies compared to the experimentally observed gel structures. For example, when using isotropic models, the formation of percolated network structures has been found to require ϕ > 0.07,[12] while, experimentally, percolation can be found for volume fractions as low as ϕ = 0.025.[13−16] Several computational studies have considered the absence of hydrodynamic interactions (HIs) as an explanation for the observed discrepancies between “dry” models and experiments.[15,17−21] While HIs have been shown to slow down the compactification of clusters, leading to more elongated local clusters in gels and percolation at ϕ = 0.04,[15] a recent analysis of simulations with and without HIs by de Graaf et al.[22] indicates that the corresponding gel structures are in fact very similar as long as the observation time is rescaled to account for the significantly different dynamics induced by HIs. Thus, it appears that structural discrepancies between experiments and modeling are not fully accounted for by HIs, and that other effects need to be incorporated into the models.

One possible explanation for the more branched network structures seen in experiments is the inhibition of lateral particle rearrangements after initial aggregation. Particle surface roughness can lead to interlocking and an effective rotational friction between two particles,[23] hampering lateral particle movement. Such rearrangements are not penalized when using interaction potentials that only depend on the separation between particles, leading to a stronger coarsening and densification of the particle network after aggregation.[15,20,24] Furthermore, for soft particles composed of cross-linked polymer networks in water, so-called microgels,[25,26] polymer ends on the particle surfaces can exhibit solvent-avoiding behavior at elevated temperatures.[27] The nature of the interparticle attraction between microgel particles is complex and contains contributions from van der Waals attraction[28] and the interaction between solvent-avoiding polymers on particle surfaces,[29] where the latter can furthermore cause polymer entanglement upon colloidal aggregation. This “locking in” of the local configuration upon aggregation is not effectively modeled using isotropic potentials. In this article, we propose a model of such soft microgel particles based on an overlay between an isotropic repulsion and a directional attraction, so-called “patchy particles”, as a tool to hinder lateral particle reorganization after gelation. This class of models has been extensively studied in the context of equilibrium systems such as polymer gels, associating fluids, proteins, and patchy colloids,[30−38] in addition to one very recent study of non-equilibrium colloidal gels formed via depletion attraction.[24] We show here that such patchy particle models are capable of capturing the experimentally observed structures in non-equilibrium colloidal gels formed through arrested spinodal decomposition of soft microgel particles interacting via short-range attractive interactions and soft repulsion.[16] Using a novel analysis technique based on a “skeletonization” of the gel backbone,[39] we compare structural properties between gels formed in simulations with those obtained from experiments on soft poly(N-isopropylacrylamide) (pNIPAm) particles that aggregate following an initial temperature-induced collapse. We show that using a directional model with sufficient number of attractive patches promotes the formation of highly branched gel structures, which reflect the experimentally obtained structures even in the absence of hydrodynamic interactions. Our results suggest that inhibited lateral particle rearrangements due to polymer entanglement and surface roughness play an important role in determining the structural signatures of colloidal gels made from soft particles, and that models based on directional attractions are suitable for describing such systems.

Results and Discussion

As described in further detail in Methods, the soft repulsion between particles is modeled as an isotropic Hertzian potential with an interaction diameter σ, which we use as our basic unit of length; this model has been used extensively to model soft microgel particles.[40,41] The repulsive particles are decorated with Np attractive patches (4 ≤ Np ≤ 12) of strength 40kBT, with kBT as the thermal energy, which is strong enough to ensure irreversible aggregation. Each patch is described through a rounded square-well-like potential, centered around regularly spaced points on the particle surface (see Figure ). In addition, we studied the corresponding isotropic attraction model where the same attractive potential instead acts along the entire particle surface, leading to each bond having the same strength as for the patchy model but with no directionality. The model was solved using Brownian dynamics (BD) simulations of N = 104 particles at a volume fraction of ϕ = 0.05, matching the approximate experimental density; for further details, see Methods. Note that the directionality imposed on the bonds in our model is a way to prevent lateral particle rearrangements after aggregation: the underlying experimental interaction potential is still effectively isotropic. Thus, to accurately model the experimental system, it is crucial that Np is large enough to ensure that virtually all particle collisions lead to a bond (i.e., that aggregation remains diffusion limited), while still being small enough to prevent lateral rearrangements through the formation of multiple bonds per particle pair. Alternative ways of constraining the particle bonds after aggregation proposed previously is the use of three-body potentials[30,42,43] or more complex simulation algorithms that create surface-localized, irreversible bonds after collision.[44,45] While aiming for the same goal, these approaches are computationally more complex than models based on attractive patches.

Figure 1

(A) Patchy particle with Np = 12 patches arranged in an icosahedral geometry. (B) Repulsive Hertzian interaction potential as a function of the center-to-center separation r (red), and the corresponding total (repulsive + attractive) potential, with the latter representing the case of two patches perfectly facing each other.

The simulation results were compared to experiments on gels formed from fluorescently labeled pNIPAm particles at ϕ = 0.05 with approximately 10 min aging time. The gel backbone was reconstructed from confocal laser scanning microscopy (CLSM) images using a skeletonization algorithm.[39,46] The experimental backbone structures were then compared with those obtained from BD simulations after the same effective aging time of 10 min. Further details of the experimental and analysis procedures are given in Methods.

Figure shows the gel backbone structures retrieved from both experimental gels and BD simulations of the isotropic and patchy particle models. Here, NBR is the number of branch points in the backbone, averaged over four simulations or 200 experimental image stacks. Comparing the latter two images (Figure B,C), it is clear that the directional interaction yields a significantly more branched and less open structure than the isotropic interaction, due to the inhibition of lateral diffusion after initial aggregation, and can better approach the experimental results. To further investigate the importance of directional bonds on the final gel structure, we performed simulations using four different valencies: Np = 4, 6, 8, and 12, arranged in tetrahedral, octahedral, cubic, and icosahedral geometries on the particle surface. In Figure A, we follow the time evolution of the gelation process by calculating the average number of nearest neighbors ⟨Z⟩ within a distance σ of each particle. It is clear that the introduction of patches strongly decreases the long-time value of ⟨Z⟩ by inhibiting particle reorganization after initial aggregation, an effect not due to the finite number of binding sites alone: even with Np = 12, ⟨Z⟩ saturates at a value of ≈4.5, which is significantly lower than the corresponding value of ⟨Z⟩ ≈ 7.5 observed for the isotropic model. The initial aggregation kinetics are somewhat slowed down for the patchy particles compared to the isotropic ones, indicating that the directionality of the interactions alters the diffusion-limited dynamics of the isotropic system, since the patchy particles need to align their binding sites before aggregating; we will see below how this unwanted effect can be overcome through the overlay between isotropic and directional interactions. The structure factors S(q) of the simulated gel structures are shown in Figure B. In accordance with our previous observations, S(q) for the isotropic system shows significantly higher values than the corresponding patchy particle gels at intermediate q, reflecting the thicker strands in the former. At high q, the specific patch geometries induce local order that results in peaks reflecting the various local structures. Finally, the low-q values of S(q) again clearly show the smaller mesh size of the gels formed from directionally interacting particles as compared to isotropic ones.

Figure 2

Gel backbones reconstructed using the skeletonization algorithm applied to (A) a stack of averaged CLSM images and (B, C) simulations using either (B) an isotropic interaction potential or (C) a directional potential with Np = 12. The number of branch points in the backbone NBR is shown in the insets, where the system volume is (47.1σ)3 in both experiments and simulations. The bottom row shows the corresponding particle configurations, with the CLSM data showing a number of 2D images from the 3D stack. The scale bar is 12σ long.

Figure 3

Time evolution of (A) the average number of nearest neighbors per particle ⟨Z⟩. (B) Structure factors S(q) of the final configuration of the gel.

In Figure A, we analyze the tortuosity of the gel skeleton, ξ, defined aswhere A and B are start and end points of paths along the gel backbone intersecting, respectively, the upper or the lower box faces in any Cartesian direction, λ(A, B) is the length of this path, while λEuc(A, B) is the Euclidian distance between A and B; the angular brackets denote averaging over all such paths for a given configuration. This measure is often used to characterize flow rates through porous catalysts[47] or blood vessel networks,[48] where a large density of branch points, high degree of branching and nonerratic strands decrease the tortuosity of the network.[49−51] In Figure A, we observe a clear decrease in the tortuosity as the gel coarsens due to a gradual straightening of gel strands. At long times, ξ relaxes to a value that is lower for higher Np as a consequence of the higher node density and degree of branching. The experimentally obtained value after 10 min of aging (dashed line in Figure A) is furthermore very close to the asymptotic value from the simulations with Np = 12, indicating a macrostructural similarity in strand geometry between experiment and simulation. We further quantify the gel structure by viewing it as a set of nodes connected by links, with ρL(Λ) as the density of links of length Λ, as shown in Figure B. The patchy potential clearly leads to a significant increase in the density of short links, corresponding to a more strongly branched network, an effect that becomes stronger with increasing Np. The experimental ρL(Λ) curve and the total number of branch points NBR are furthermore significantly higher than those of its isotropic model counterpart, reflecting the sparse network structure seen in Figure B. The slight shift in peak position and width in simulations compared to the experimental curves is likely caused by optical limitations inherent to CLSM.[52]

Figure 4

(A) Tortuosity ξ, defined in eq , as a function of time for the patchy particle simulations with Np = 12 and Np = 8, and isotropic particle simulations. The experimentally obtained value, as obtained from approximately 200 experimental image stacks at 10 min aging time, is indicated by the dashed line. (B) Link density ρL(Λ) as a function of the link length Λ for the final gel structures formed using the isotropic and the directional potential with Np = 12 and Np = 8, as obtained from the final 150 frames (approx. 30 s) averaged over four separate simulations. The experimental value (dashed line) was obtained as in (A). NBR denotes the average total number of branch points in the system. Data for low-valency patchy particle simulations are omitted due to high statistical noise.

As previously observed in Figure A, a limitation of using a directional interaction model to describe intrinsically isotropic particles is the possibility of collisions not leading to bond formation, which stops the initial aggregation kinetics from being accurately described. In Figure , we therefore depict the results from an overlay of the patchy attraction (of strength UP) with an isotropic attraction (of strength UI), constraining the total energy per bond to UP + UI = 40kBT. The backbones, as seen in Figure A–F, clearly show that the branching density increases with decreasing UI. ⟨Z⟩ (Figure G) furthermore shows that the addition of a relatively small isotropic attraction (UI = 10kBT) changes the initial aggregation kinetics to closely resemble that of the purely isotropic system, while still reaching the same long-time value of ⟨Z⟩ as in the purely directional case. With further increase of UI, the late-stage dynamics will also change, as particles can now start rearranging by escaping from their attractive wells, defined by the strength of the directional bonds. The corresponding ρL(Λ) curves (Figure H) reflect the trends seen in Figure A–F and confirm that changing UI relative to UP allows for a detailed control over the gel structure. The addition of a small isotropic attraction (UI = 10kBT) causes a small increase in branching, likely caused by the transition from partially reaction-limited to purely diffusion-limited aggregation. From there, increasing UI decreases the degree of branching due to the lateral rearrangement mechanism discussed earlier. Thus, using a combination between isotropic and directional interactions provides a route to preserving the diffusion-limited aggregation present in the isotropic model, while being able to reproduce the more branched networks as observed in experiment. Furthermore, tuning the relative magnitudes of isotropic and directional contributions allows fine-tuning of the network branching.

Figure 5

Overlaying patchy and isotropic potentials allows for increased control of structure and aggregation kinetics. (A–F) Simulated backbone structures after 600 s aging time. (G) Time evolution of ⟨Z⟩ for the systems in (A–F). (H) Corresponding link density ρL(Λ); the dashed line shows the same experimental curve as in Figure B. In all cases, the number of patches is Np = 12, and the sum of the patchy and isotropic potential depths is kept constant at UP + UI = 40kBT.

Conclusions

In this article, we have described a computational model of colloidal gels formed through arrested spinodal decomposition based on particles with directional bonds (patchy particles). The simulation results were compared to experiments on gels formed from soft pNIPAm microgel particles through the use of a novel skeletonization algorithm that reconstructs a three-dimensional image of the gel backbone from a collection of microscopy images, providing a versatile method for analyzing gel structures, including systems that do not offer single-particle resolvability in microscopy. Our results show that the use of attractive patches strongly promotes branching as compared to a purely isotropic model. This strong effect points toward the importance of inhibiting lateral particle rearrangements after the initial aggregation due to surface roughness or entanglement between polymer strands on the touching particle surfaces. Overlaying isotropic and directional attractions furthermore allows for fine-tuning the structure, which can be adjusted to resemble experimental gel structures.

Unlike previous studies of the (equilibrium) phase behavior of patchy particles,[33,36] where the role of the directional bonds is to describe an underlying anisotropic interaction potential, the patchy potential here serves as a simple means of closing a kinetic pathway toward gel coarsening present in the experimental system, even though the underlying interaction potential is approximately isotropic. In common with previous approaches based on three-body potentials,[53] this induces an artificial short-range order dictated by the chosen patch geometry. One way to partially circumvent this drawback would be the use of randomly distributed patches; an approach that was very recently investigated in ref (24). The use of such irregular patch geometries however makes it difficult to prevent the formation of multiple bonds between patch pairs on the same two particles, and thus leading to a potentially ill-defined bond strength. The promising results presented here thus call for further investigations of the effects of different model approximations on the structure, coarsening dynamics, and rheology of model colloidal gels.

Model and Methods

Particle Model

The total pairwise interaction between microgel particles was described as the sum of a soft isotropic Hertzian repulsion UH, given byand a rounded square-well-like attraction UA describing short-range attractionswhich was either isotropic or localized. In the above equations, σ is the effective particle diameter, ϵH = 5 × 103kBT sets the strength of the soft repulsion, with kBT as the thermal energy, and r is the center-to-center distance between two particles (for UH) or two patches (for UA). In the attractive potential, ϵA describes the maximum attraction strength, while and specify the width and steepness of the attraction, respectively. The value of the repulsive strength ϵH was determined along the lines of previous work, after being adjusted to prohibit multiple-patch bonds,[16] while the attractive strength ϵA = 40kBT was chosen large enough to ensure effectively irreversible particle bonds, in accordance with what we observe in experiments. The attractive, massless patches were distributed on the particle surface at positions corresponding to vertices of regular polyhedra, with each patch centered around point lying 0.4875σ outside the particle center, to prevent bonds between multiple patches. For the isotropic interaction case, the potential described by eq instead acted over a spherical shell centered at the same distance from the particle center.

Simulation Details

All simulations were performed using the open-source LAMMPS molecular dynamics package,[54] using overdamped Langevin dynamics for the translational and angular velocities v and ω of particle iwhere F and T are the force and torque on particle i, respectively, due to interparticle interactions, Dt is the translational self-diffusion constant, and Dr = 3Dt/σ2 is the rotational diffusion, where we assume each particle to rotate as a solid sphere of diameter σ. ηt(t) and ηr(t) are the translational and rotational noises with unit variance, δ-correlated in space and time. In practice, we achieved the dynamics in eqs and 5 by solving the full (underdamped) Langevin equation, but with a damping time short enough to ensure that inertia is fully negligible, thus leading to effective overdamped dynamics. All simulations were run with N = 104 identical particles in a periodic box of linear dimensions L = 47.1σ, yielding a volume fraction of ϕ = 0.05 to mimic experiment. The initial configuration was one of randomly distributed particles, barring any unphysical particle overlap, thus resembling the experimental systems before aggregation due to particle collapse. All simulations were run four times with different initial configurations to obtain estimates of the variance in the different measured properties. Visualizations were obtained using the Open Visualization Tool OVITO.[55]

To map the units of our simulations to experimental time and length scales, we rescaled all simulation and experimental results by σ or the experimentally measured hydrodynamic diameter at 50 °C (see below) and the Brownian time τB, defined aswhere μ is the dynamic viscosity. The time step was set to 2.4 × 10–4 τB, corresponding to 1.2 × 10–5 s at 50 °C.

Experimental Procedures

The particles were quenched instantaneously from 20 to 50 °C in 10–1 M KCl, leading to the collapse of the individual particles followed by aggregation.[56] The gels were allowed to rest for 10 min, after which they were imaged using CLSM. The pNIPAm particles used in this work were prepared in accordance with a previously described protocol.[16] Hydrodynamic diameters were determined from dynamic light scattering using a modulated 3D cross-correlation instrument (LS instruments) with a 660 nm diode-pumped laser and determined to 620 and 354 nm at 20 and 50 °C, respectively. It should be noted that the hydrodynamic radii were obtained at 10–3 M KCl to prevent aggregation for temperatures above the volume phase transition temperature, leading to a slight overestimation of the radii compared to fully screened conditions. The number density of the stock solution was determined at 20 °C by counting the particles in at least 30 CLSM image stacks of fully crystalline samples. The number density together with the measured hydrodynamic radius was then used to calculate the volume fraction at 50 °C.[57] A CLSM sample slide containing a ϕ = 0.05 dispersion in 10–1 M KCl was heated to 50 °C by contact with a preheated metal block and imaged after 10 min using an inverted Leica TCS SP5 tandem scanner using a 100× oil immersion objective, in an enclosure that allows for temperature control with a 0.2 °C maximum variance using thermostatted air circulation.

Skeletonization Analysis

Raw CLSM images for skeletonization were treated with Gaussian blurring using cubic kernels of size 5 voxel lengths (0.6σ) to remove speckled noise at gel surfaces, and subsequently binarized with thresholding according to Otsu’s method.[58] BD simulations were transformed into similar binarized image stacks from coordinate diagrams by projecting particle positions and sizes onto image slices. A morphological closing step was performed using a structural element of size 1 voxel, followed by removal of small unconnected structures (<1% of largest cluster volume). These procedures involve Minkowski-based procedures adapted into MATLAB.[59] The skeletonization was then performed using previously developed algorithms[39,46] to create a pixel-thick backbone from binarized images. The skeleton was subsequently transformed into a set of nodes and links followed by a “cleaning procedure,” which repeatedly (i) removes branches connected to only one or two other nodes, excluding box edge nodes and (ii) fuses nodes within a separation smaller than a particle diameter, until a structure emerges that cannot be cleaned further. For ρL(Λ) calculations, links to the box edge nodes were not considered.