Collective dynamics and self-organization of active particles on random networks

TL;DR

This study models how active particles on a network, representing cells in ECM, undergo a phase transition from isolated clusters to a large, dynamically reorganizing cluster driven by force percolation and durotaxis.

Contribution

Introduces a minimal active-particle-on-network model revealing a critical transition in collective behavior driven by force network percolation and cell motility mechanisms.

Findings

Transition from isolated to large clusters at critical density

Radical particles dominate dynamic reorganization

Force percolation underpins collective cell behavior

Abstract

Collective cell migration in 3D extracellular matrix (ECM) is crucial to many physiological and pathological processes. Migrating cells can generate active pulling forces via actin filament contraction, which are transmitted to the ECM fibers and lead to a dynamically evolving force network in the system. Here, we elucidate the role of such force network in regulating collective cell behaviors using a minimal active-particle-on-network (APN) model, in which active particles can pull the fibers and hop between neighboring nodes of the network following local durotaxis. Our model reveals a dynamic transition as the particle number density approaches a critical value, from an "absorbing" state containing isolated stationary small particle clusters, to a "active" state containing a single large cluster undergone constant dynamic reorganization. This reorganization is dominated by a subset of highly dynamic "radical" particles in the cluster, whose number also exhibits a transition at the same critical density. The transition is underlaid by the percolation of "influence spheres" due to the particle pulling forces. Our results suggest a robust mechanism based on ECM-mediated mechanical coupling for collective cell behaviors in 3D ECM.

Full text

Collective cell migration is crucial to many physiological and pathological processes such as tissue regeneration, immune response and cancer progression ref1 ; ref2 ; ref3 ; ref4 . Cell migration in 3D extracellular matrix (ECM) is a complex dynamic process involving a series of intra-cellular and extra-cellular activities ref12 ; ref13 , and can be regulated by a variety of cell-ECM interactions via chemotaxis ref14 , durotaxis ref16 , haptotaxis ref17 , and contact guidance ref18 . A migrating cell also generates active pulling forces, which are transmitted to the ECM fibers via focal adhesion complexes ref21 and consistently remodel the local ECM (e.g., by re-orienting the collagen fibers, forming fiber bundles and increasing the local stiffness of ECM) ref24 ; ref29 ; chen2019 . In a multi-cell system, the pulling forces generated by individual cells can give rise to a dynamically evolving force network (carried by the ECM fibers) ref8 ; Frey07 . The force network can further influence the migration of the cells, which in turn alters the ECM and the force network ref5 ; ref6 ; ref7 ; ref8 ; ref9 ; ref10 ; ref11 . This feedback loop between the force network and cell migration could lead to a rich spectrum of collective migratory behaviors.

The cell-ECM system is an example of complex many-body systems in which the individuals (e.g., migrating cells) communicate and interact with one another through their environment (e.g., ECM), while simultaneously re-shaping the environment, altering the means (e.g., the force network) to pass information among themselves. Other examples of such complex systems include flocks of birds, schools of fish, and active swimmers in crowded environment Popkin16 ; Bechinger16 . In this letter, we investigate the collective cellular dynamics and self-organizing multi-cellular patterns in 3D ECM resulted from the dynamically evolving force network, using a minimal active-particle-on-network (APN) model. Although focusing on the cell-ECM system, the physical insights obtained here are also valuable to understanding active-particle systems dominated by environment-mediated particle-particle interactions.

Active-particle-on-network model: The 3D ECM is modeled as a discrete network with a “graph” (i.e., node-bond) representation in a cubic simulation domain with linear size $L$, which is composed of $M_{n}$ nodes and $M_{b}$ bonds. The average coordination number $Z$, i.e., the average number of bonds connected to each node, is given by $Z=2M_{b}/M_{n}$. We have used both the periodic boundary (PB) conditions and fixed boundary (FB) conditions (i.e., the nodes within a certain distance $\delta L$ from the boundaries of the simulation domain are fixed) in our simulations, and find that even for a moderate system size (e.g., $M_{n}\sim 5000$) the boundary conditions do not affect the results. In the subsequent discussions, we will mainly present the results obtained using the fixed boundary conditions, under which $M_{n}$ denotes the number of free (non-fixed) nodes.

Next, $N_{p}$ active particles (e.g., congruent spheres) are introduced in the network such that each particle occupies a randomly selected un-occupied node (i.e., each node can be occupied by only one particle). The number density $\rho$ of the particles is defined as $\rho=N_{p}/M_{n}$, i.e., the fraction of nodes occupied by the active particles. Each particle can generate a contraction, which pulls all of the bonds connected to the node it occupies towards the particle center (i.e., the node) by a fixed amount $\delta l$, leading to different pulling forces in the bonds and thus, a force network in the system. We consider the particles can “migrate” from its original node to an un-occupied neighboring node following local durotaxis, i.e., along the bond with the highest stiffness, which is also the bond that carries the largest pulling force among all neighboring bonds (see Fig. 1 for illustration). The diameter of the particles is not essential in our model and thus, is not explicitly considered here.

The bonds of the network are modeled as elastic elements with only non-zero stretching modulus $E_{s}$ and free to rotate at the nodes. An active particle can generate pulling forces in the bonds connected to the node it occupies by contraction, i.e., $\delta l$. This contraction leads to a strain $\epsilon_{i}=\delta l/l_{i}$ in the bond $i$ with original length $l_{i}$, and thus, a pulling force $f_{s}=E_{s}A\delta l/l_{i}$, where $A$ is an effective cross-sectional area of the bonds. These pulling forces impose force boundary conditions for the ECM network, and the force-balance network configuration is obtained using an iterative force-based relaxation approach ref11 .

We note that many factors can affect the interactions between the dynamic force network and the collective dynamics of the active particles in our APN model. These may include the geometry/topology and mechanical properties of the network, as well as the number density, spatial distribution and the contractibility (i.e., $\delta l$) of the active particles. In this work, we mainly focus on disordered isostatic networks (i.e., $Z=6$) derived from maximally random jammed packings of congruent hard spheres ref32 . It is straightforward to generalize this study to random network models derived from confocal images of collagen gels ref33 . In addition, we use simple linear elastic network models (see Methods). This allows us to investigate the system in the elastic regime MacKintosh05 , in which the force network is mainly determined by the number density and spatial distribution of the active particles, and largely independent of particle contractibility. Our model can readily incorporate more realistic mechanical models for the ECM, taking into account non-linear responses of the fibers Safran12 ; nat_method15 and plasticity ref29 . Moreover, in an actual cell-ECM system, the cell migration might not be sensitive to individual stiffer fibers, but determined by certain meso-scale stiff structures emerged due to cell remodeling, such as bundles of high-stress fibers. Nonetheless, we believe that the general organizational principles of active particles on random networks obtained here are relevant to and can provide insights on the actual cell-ECM systems.

Dynamic phase transition in the APN system: We now describe the observed collective dynamics of the active particles on the random networks. In our simulations, we systematically vary the particle number density $\rho\in(0.05,0.95)$. For each $\rho$, the particles are initially randomly introduced in the network and the system is allowed to evolve according to the aforementioned APN dynamics. At low densities (i.e., $\rho<\rho_{c}\approx 0.114$), the particles rapidly aggregates into multiple isolated small clusters, which are randomly distributed within the ECM (see Fig. 2a and 2b). Here, we consider two particles belong to the same “cluster” if they occupy the two nodes connected by the same bond. As $\rho$ approaches $\rho_{c}$ from below, the maximal cluster size $n_{c}$ (i.e., the number of particles in the largest cluster of the system) increases dramatically (see Fig. 2c), indicating the majority of the particles are connected to form a single large cluster in the system.

In addition, we find that the isolated small clusters associated with $\rho<\rho_{c}$ are stationary, i.e., the particles in the clusters either do not move at all or hop between two adjacent nodes (typically at the boundary of a cluster). On the other hand, the dominant large clusters formed for $\rho>\rho_{c}$ undergo constant dynamic reorganization. To further quantify the dynamics of the clusters, we count the number of distinct nodes $m_{s}$ visited by a particle during a total of $s$ successive steps. The collected statistics for different particle densities are shown in Fig. 3a. We note that the $m_{s}$ statistics shown in Fig. 3a does not depend on $s$ and we have used $s=24$ here.

It can be seen from Fig. 3a that for small $\rho$, a particle can only visit one or two nodes, respectively indicating that the particle does not move or can hop between two adjacent nodes. As $\rho$ approaches $\rho_{c}$ from below, although the majority of particles are localized (indicated by the peak in the $m_{s}$ statistics associated with small node numbers), a subset of highly dynamic particles emerge which are able to visit many distinct nodes for a given number of steps (indicated by the emergence of the second peak associated with large node numbers in the $m_{s}$ statistics). We refer to these highly dynamic particles as “radicals”, i.e., those do not possess a periodic hopping pattern overall a finite number nodes. The trajectory of a small number of randomly selected radicals are shown in Fig. 3b for different $\rho$ values. Fig. 3c shows the number of radicals $N_{r}$ as a function of $\rho$. It can be seen that $N_{r}$ exhibits a clear transition as $\rho$ increases towards $\rho_{c}$. This is consistent with the transition observed in the maximal cluster size $n_{c}$ as $\rho$ increases (see Fig. 2c).

The above analysis suggests that the system possesses a phase-transition-like behavior, as the particle number density $\rho$ increases, from an “absorbing” state in which the particles segregate into small isolated stationary clusters, to a “active” state, in which the majority of particles join in a single large dynamic cluster. This transition is also quantitatively manifested in the maximal cluster size $n_{c}(\rho)$ (Fig. 2c) and radial number $N_{r}(\rho)$ (Fig. 3c) as $\rho$ increases towards $\rho_{c}\approx 0.114$. In particular, our scaling analysis shows that as $\rho_{c}$ is approached from above, $n_{c}\sim(\rho-\rho_{c})^{\alpha}$, where the critical exponent $\alpha\approx 6.6\pm 0.1$. In addition, we find that $N_{r}\sim(\rho-\rho_{c})^{\beta}$, where the critical exponent $\beta\approx 1.63\pm 0.04$. The numerical values of $\alpha$, $\beta$ and $\rho_{c}$ are obtained by fitting the simulation data. We also note that the $\rho_{c}$ value is much lower than the site percolation threshold for the network ($\approx 0.310$) percolation . The absorbing-to-active transition has also been observed in a wide spectrum of “random-organizing” physical systems, such periodically driven colloids randorg01 ; randorg02 ; randorg03 , granular materials and amorphous solids randorg04 ; randorg05 ; randorg06 , vortices randorg07 and skyrmion systems randorg08 .

Influence sphere due to active pulling forces: We now investigate the mechanisms for the observed transition. Once a particle pulls the fibers, a stress gradient is built up surrounding this particle. When another particle “senses” the pulling force ref10 , it will tend to move up the stress gradient towards the contracting particle due to local durotaxis. This would lead to an effective mutual pulling between the particles.

It is reasonable to assume that the pulling forces generated by a specific particle can only influence other particles within a certain distance $R_{I}$. Due to the intrinsic network heterogeneity, $R_{I}$ may vary for different particles. Here, we take a “mean-field” approach and assign the same effective $R_{I}$ to all the particles in the system and introduce the concept of the influence sphere, which is a spherical region with radius $R_{I}$ centered at a contractile particle (see the inset of Fig. 4).

The influence-sphere radius $R_{I}$ is estimated from the cluster statistics of the APN systems at low $\rho$, i.e., those containing multiple small isolated stationary clusters in the final state. This is based on the assumption at low $\rho$, only particles which are within the influence region of one another would eventually aggregate. In particular, we first identify the particles within the same cluster in the final state of the system. Then the system is “re-winded” to the initial state, and the intra-cluster nearest-neighbor distances $d_{n}$ are computed for all clusters. We then use the mean nearest-neighbor distance $\bar{d}_{n}$ to estimate $R_{I}\approx 0.104L$ (where $L$ is the linear system size), which is roughly twice of the average fiber length (see SI for fiber length distributions).

Mean field theory: Percolation of influence sphere: We now investigate the percolation of the influence spheres as $\rho$ increases. For a given $\rho$, we randomly place particles on the ECM network. Instead of allowing the particles to move according to the APN dynamics, we place a virtual sphere with radius $R_{I}=0.104L$ at each particle, representing the influence spheres. We subsequently identify the clusters formed by the influence spheres , based on which the percolation of the system can be determined.

Fig. 4 shows the percolation probability analysis for the system ref34 (see SI for details), from which a percolation transition and the associated critical density (i.e., percolation threshold) $\rho^{*}\approx 0.109$ can be clearly identified. Interestingly, the percolation transition of the influence spheres coincides with the dynamic transition of the active particles at $\rho_{c}\approx 0.114$. This suggests that the dynamic transition of the active particles from the “absorbing” state to the “active” state is underlaid by and can be understood as the percolation transition of the influence spheres.

We also investigate the collective dynamics of active particles on other network models, including the networks derived from the diamond lattice, as well as random networks reconstructed based on confocal images (see SI for details). We find that although the critical transition density $\rho_{c}$ depends on the geometry and topology of the networks, the dynamic transition and the aggregation behaviors of the active particles occur in all of the studied model networks. This suggests the feedback loop between the evolving force network and particle dynamics provides a robust mechanism for regulating collective migratory behaviors. In future work, we will also explore the effects of fiber alignment and external mechanical cues.

Finally, we emphasize again that this minimal model does not take into account crucial mechanisms in actual cell migration such as ECM remodeling (e.g., orientation, bundling and degradation) and cell-cell adhesion. Interestingly, our studies indicate that, at least for the APN systems, the local durotaxis for the active particles is sufficient to induce and stabilize aggregations, even without adhesion. Nonetheless, we expect that the insights on the collective behaviors of active particles regulated by the dynamically evolving force network obtained here are helpful in understanding the collective dynamics emerged in actual multi-cellular-ECM systems, as well as in other active-particle systems dominated by environment-mediated particle-particle interactions.