{'fraction_non_alphanumeric': 0.025276533503976597, 'fraction_numerical': 0.013438157052747052, 'mean_word_length': 5.17005076142132, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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][6][7][8][9][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][14][15][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][18][19][20][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 modelling are not fully accounted for by HIs, and that other effects need to be incorporated into the models.
One possible explanation for 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 lat-eral 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 crosslinked 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][31][32][33][34][35][36][37][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-equlibrium colloidal gels formed through arrested spinodal decomposition of soft microgel particles interacting via short-range attractive interactions and soft repulsion. 16 By 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 a 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 the Methods section, the soft repulsion between particles is modelled 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 N p attractive patches (4 ≤ N p ≤ 12) of strength 40k B T , with k B T the thermal energy, 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 Fig. 1). 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 = 10 4 particles at a volume fraction of φ = 0.05, matching the approximate experimental density; for further details, see the Methods section. 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, in order to accurately model the experimental system it is crucial that N p 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 towards the same goal, these approaches are computationally more complex than models based on attractive patches.
The simulation results were compared to experiments on gels formed from fluorescently labeled pNIPAm particles at φ = 0.05 with approximately 10 minutes 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 minutes. Further details of the experimental and analysis procedures are given in the Methods section. The 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.
simulations of the isotropic and patchy particle models. Here, N BR is the number of branch points in the backbone, averaged over four simulations or 200 experimental image stacks. Comparing the latter two images (Figs. 2B and 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. In order to further investigate the importance of directional bonds on the final gel structure, we performed simulations using four different valencies: N p = 4, 6, 8 and 12, arranged in tetrahedral, octahedral, cubic and icosahedral geometries on the particle surface. In Fig. 3A, 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 N p = 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 Fig. 3B. 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 causes 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.
In Fig. 4A, we analyze the tortuosity of the gel skeleton, ξ , defined as  . N BR 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. more strongly branched network, an effect which becomes stronger with increasing N p . The experimental ρ L (Λ) curve and the total number of branch points N BR is furthermore significantly higher than its isotropic model counterpart, reflecting the sparse network structure seen in Fig. 2B.
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 As previously observed in Fig. 3A, 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 Fig. 5, we therefore depict the results from an overlay of the patchy attraction (of strength U P ) with an isotropic attraction (of strength U I ), constraining the total energy per bond to U P + U I = 40k B T .
The backbones, as seen in Fig. 5A-F, clearly show that the branching density increases with decreasing U I . Z , (Fig. 5G) furthermore shows that the addition of a relatively small isotropic attraction (U I = 10k B T ) 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. By increasing U I further, 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 (Fig. 5H) reflects the trends seen in Fig. 5A-F, and confirms that changing U I relative to U P allows for detailed control over the gel structure. The addition of a small isotropic attraction (U I = 10k B T ) causes a small increase in branching, likely caused by the transition from partially reaction-limited to purely diffusion-limited aggregation. From there, increasing U I 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.  Fig. 4B. In all cases, the number of patches is N p = 12, and the sum of the patchy and isotropic potential depths is kept constant at U P +U I = 40k B T .

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 3-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 towards the importance of inhibiting lateral particle rearrangements after 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 of 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 towards 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

Particle model
The total pairwise interaction between microgel particles was described as the sum of a soft isotropic Hertzian repulsion U H , given by and a rounded square-well like attraction U A describing short-range attractions, which was either isotropic or localized. In the above equations, σ is the effective particle diameter, ε H = 5 × 10 3 k B T sets the strength of the soft repulsion, with k B T the thermal energy, r is the center-to-center distance between two particles (for U H ) or two patches (for U A ). In the attractive potential, ε A describes the maximum attraction strength, while w = σ 20 and k = σ 200 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 = 40k B T 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, in order to prevent bonds between multiple patches. For the isotropic interaction case, the potential described by Eq. (3) 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 i and ω ω ω i of where the F i and T i is the force and torque on particle i due to interparticle interactions, D t is the translational self-diffusion constant and D r = 3D t /σ 2 is the rotational diffusion, where we assume each particle to rotate as a solid sphere of diameter σ . η η η t (t) and η η η r (t) represent translational and rotational noise with unit variance, δ -correlated in space and time. In practice, we achieved the dynamics in Eqs. (4)-(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 = 10, 000 identical particles in a periodic box of linear dimensions L = 47.1σ , yielding a volume fraction of φ = 0.05 in order 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 in order to obtain estimates of the variance in the different measured properties. Visualizations were obtained using the Open Visualization Tool OVITO. 55 In order 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 as where µ is the dynamic viscosity. The timestep was set to 2.4 × 10 −4 τ B , corresponding to 1.2 × 10 −5 seconds at 50 • C.

Experimental procedures
The particles were quenched instantaneously from 20 • C 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 minutes, 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 were determined to 620 nm and 354 nm at 20 • C and 50 • C, respectively. It should be noted that the hydrodynamic radii were obtained at 10 −3 M KCl in order to prevent aggregation for temperatures above the volume phase transition temperature, leading to a slight overestimation of 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 minutes 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.