Aggregation-diffusion Equations for Collective Behavior in the Sciences

A surprisingly large number of complex phenomena can be modeled as systems of point particles that interact under short- and long-range forces. While researchers typically use this modeling framework to study many-body systems (such as gravitational collapse, electron transport in semiconductors, and kinetic equations in statistical physics), it has recently found applications in other areas of physics as well as biology, social science, machine learning, and even optimization. These applications are diverse and span both microscopic and macroscopic scales; examples include ion channel transport, chemotaxis, bacterial orientation, cellular adhesion, angiogenesis, animal herding, and human crowd motion. One can approximate the interaction effects between the “individuals” in each case with long-range attractive forces (e.g., ligand binding, electrical interactions, or social preferences) and short-range repulsions (due to volume constraints or crowding).

Several agent-based (microscopic) modeling approaches—such as cellular automata and driven Brownian particles—can describe these phenomena and encode complex behavior. However, the best way to understand and rigorously analyze such behavior is to study the mean-field limit (or its variants) — i.e., the limit of the microscopic model as the number of agents becomes large. For instance, consider two types of cells, both of which are able to express different surface proteins and ligands (such as cadherins or nectins). Assume that the nuclei of the first cell type are located at positions \(\{x_i\}\), \(i=1,\dots,N\), and the nuclei of the second type are located at \(\{y_j\}\), \(j=1,\dots,M\); let also \(N=M\) for simplicity. We consider cells that interact either by attraction at medium distances (e.g., through filopodia, which are protrusions of the cellular membrane) or by strong repulsion at short distances due to volume and size constraints (see Figure 1). If the force exerted by cells of type \(b\) onto type \(a\) is radial and conservative (with potential \(W^N_{ab}\), where \(a,b=1,2\)), we may pose the elementary agent-based model 

\[\dot x_i = - \frac{m_1}N \sum_{j\neq i}\nabla W^N_{11}(x_i-x_j) - \frac{m_2}N\sum_{j} \nabla W^N_{12}(x_i-y_j),\]

\[\dot y_i = - \frac{m_2}N \sum_{j\neq i} \nabla W^N_{22}(y_i-y_j) - \frac{m_1}N\sum_{j} \nabla W^N_{21}(y_i-x_j),\]

where \(m_1\) and \(m_2\) are the typical cellular masses of each cell type.

Under suitable assumptions, the model’s associated empirical measures (weighted sums of Dirac deltas that are supported at each particle’s position) can approximate the macroscopic normalized cell densities in the many-particle limit:

\[\rho_1(t,x) \simeq \frac{m_1}N \sum_{i=1}^N \delta_{x_i(t)}{(x)} \quad \textrm{and} \quad \rho_2(t,x) \simeq  \frac{m_2}N \sum_{i=1}^N \delta_{y_i(t)}{(x)} \quad \textrm{as} \quad N\to \infty.\]

Given the biological context, it is reasonable to suppose that repulsion is purely localized and the attraction forces only act within a cutoff radius \(R\). A natural scaling choice for the potentials is therefore \(W^N_{ab}\simeq \varepsilon \delta_0 + W_{ab}\) as \(N \to \infty\) [7], where \(\varepsilon\) relates to the typical volume that each cell occupies (here we assume that \(\varepsilon\) remains equal across species and interactions). In this scaling, the limit \(N \to \infty\) yields the mean-field model

\[\begin{equation}
\left\{
\begin{aligned}
\partial_t \rho_1 &= \nabla \cdot\Big(\rho_1 \nabla \big(\varepsilon (\rho_1+\rho_2) + W_{11}\ast\rho_1 + W_{12}\ast\rho_2 \big)\Big),\\
\partial_t \rho_2 &= \nabla \cdot\Big(\rho_2 \nabla \big(\varepsilon (\rho_1+\rho_2) + W_{22}\ast\rho_2 + W_{21}\ast\rho_1 \big)\Big).
\end{aligned}
\right.
\end{equation} \tag 1\]

Each equation in this system is a particular case of the aggregation-diffusion equation

\[\frac{\partial\rho}{\partial t}+\nabla \cdot \left(\rho u\right)=0, \quad u=-\nabla U' \left ( \rho \right )-\nabla V -\nabla W*\rho, \tag2\]

which is a generic continuum equation for the kinematic evolution of a population density \(\rho(t,x)\). The velocity field \(u\) drives the density’s evolution through a balance of three forces: (i) diffusion due to localized repulsion effects or noise, which is modeled by the convex function \(U(\rho)\); (ii) drift due to external forces, which is described by the potential \(V(x)\); and (iii) a symmetric interaction force (aggregation, repulsion, or both) that is given by the potential \(W(x)\). Linear diffusion is commonly written as \(U(s)=s\log s\), while nonlinear diffusion is written as \(U(s)={s^m}/{(m-1)}\) where \(m>0\) [14]. The diffusion coefficient in these scenarios is \(m \rho^{m-1}\). For \(m > 1\), diffusion is stronger when the concentration is large; for \(m < 1\), it is stronger when the concentration is small. So-called confining potentials \(V = |x|^p\) for \(p > 0\) often appear in the Fokker-Planck and McKean-Vlasov equations.

Gradient Flows

We can understand the aggregation-diffusion equation in \((2)\) as the gradient flow of the free-energy functional

\[\mathcal F[\rho] = \int_{\mathbb{R}^d} U(\rho (x)) \mathrm{d}x + \int_{\mathbb{R}^d} V(x) \rho(x) \mathrm{d}x + \frac{1}{2}\int_{\mathbb{R}^d}\int_{\mathbb{R}^d} \rho(x) W(x-y) \rho(y) \mathrm{d}y  \mathrm{d}x \tag3\]

in a suitable metric space, i.e., the 2-Wasserstein space of probability densities over \(\mathbb{R}^d\) [15]. This space arises when we equip the set of probability densities with a finite second moment \(\mathcal{P}_2(\mathbb{R}^d)\) that has 2-Wasserstein distance \(d_2\). This distance originates from the theory of optimal transport and has found fruitful applications in the study of mean-field limits because it induces a topology where the evolution of empirical measures in time is continuous (unlike for Sobolev spaces). As a result, \(d_2\) can help to quantify particle systems’ convergence to their mean-field limits. This notion of distance for probabilities is also useful when defining transport-based interpolations between densities as geodesic curves in the 2-Wasserstein space; Figure 2 illustrates the geodesic interpolation between the normalized characteristics of two sets. 

Due to the lack of linear structure in the 2-Wasserstein space, it is difficult to precisely construct solutions to the aggregation-diffusion equation in \((2)\) as gradient flows of the free-energy functional in \((3)\). We can generate a discrete-in-time sequence of measures from \(\rho_0\in{\mathcal{P}}_2(\mathbb{R}^d)\) via the iteration

\[\rho_{k+1}= \underset{{\rho\in {\mathcal P}_2(\mathbb{R}^d)}}{\textrm{argmin}} \left\{ \frac{1}{2 \Delta t}d_2^2(\rho,\rho_{k})+\mathcal F[\rho] \right\}\]

for \(k\in \mathbb{N}\) and any fixed \(\Delta t>0\). A suitable time interpolation leads to a curve of measures \(\rho_{\Delta t}\) that converge to the unique solution of \((2)\) as \(\Delta t \to 0\) [1, 10].

If the functional \(\mathcal F[\rho]\) has a unique global minimizer in \(\mathcal{P}_2 (\mathbb{R}^d)\), then the solution to \((2)\) converges over time to the minimizer whenever the functional satisfies a suitable notion of convexity along the geodesics that are induced by \(d_2\). However, this condition may not be accurate; the purely diffusive setting \((V=W=0)\) yields asymptotic spreading of \(\rho(t)\) in self-similar profiles for large times [14]. The asymptotic balance between aggregation and diffusion is generally quite subtle, and researchers have devoted significant efforts to its study. Under certain conditions, they have established the existence of radially decreasing, compactly supported global minimizers of \(\mathcal{F}[\rho]\) [6]. Furthermore, the minimization of \(\mathcal{F}\) in some extreme cases causes the formation of Dirac deltas that accumulate only a fraction of the total mass. Understanding these features is a considerable challenge for the dynamical problem [5, 9]. Almost-sharp conditions on \(W\) exist such that \(\rho(t)\) behaves like the solution to the purely diffusive problem. Yet in some settings, the solution can asymptotically develop the Dirac delta that is suggested by the free energy’s minimization.

Mathematical Biology Applications 

Mathematical biology is a prominent application area for aggregation-diffusion systems. Cell population models are particularly relevant and often use the aggregation-diffusion equation (in either scalar or system form) to describe the motion of cells as they concentrate or separate from a target. The aforementioned diffusion terms are consistent with population pressure effects, whereby groups of cells naturally gravitate away from areas of high concentration. In addition, the aggregation term has classically served as a model for chemotaxis, as cells often direct their movement along the increasing gradient of a chemical agent that is produced by either external sources or the cells themselves. The Keller-Segel model for chemotaxis [12] is an archetypal example of the aggregation-diffusion equation in mathematical biology.

A new paradigm in cell adhesion models suggests that short-range interaction potentials can capture ligand binding through filopodia (as in Figure 1). One example of such a model [7] considers two species of cells, each of which is attracted to both itself and the other species in the form of \((1)\). Despite its apparent simplicity, this model accounts for several realistic behaviors. If all of the attraction strengths are equal, the cells simply mix evenly. But if the relative strengths vary, intricate phenomena will arise; for instance, one species of cells might completely envelop the other and form a “halo.” Figure 3 depicts the astonishing match between the model solutions and the experimental results [11]. The model demonstrates that the differential adhesion hypothesis [13] is the driving mechanism for cell sorting [8], which we modeled via the asymmetry on the interaction potentials between the two cell types. Aggregation-diffusion systems \((1)\) with this type of asymmetry do not typically have a gradient flow structure and can develop a wealth of complex phenomena, such as traveling wave solutions, free boundaries, and mixing or segregation between species. Each of these phenomena presents a challenging open problem in the analysis of partial differential equations.

José A. Carrillo delivered an invited lecture about this research at the 10th International Congress on Industrial and Applied Mathematics, which took place in Tokyo, Japan, last year.

Acknowledgments: The authors were supported by the Advanced Grant Nonlocal-CPD (Nonlocal PDEs for Complex Particle Dynamics: Phase Transitions, Patterns and Synchronization) of the European Research Council Executive Agency under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 883363). In addition, David Gómez-Castro was partially supported by grants RYC2022-037317-I and PID2023-151120NA-I00 from the Spanish government.

References
[1] Ambrosio, L., Gigli, N., & Savaré, G. (2008). Gradient flows in metric spaces and in the space of probability measures. In Lectures in mathematics ETH Zürich. Basel, Switzerland: Birkhäuser.
[2] Bailo, R., Carrillo, J.A., & Hu, J. (2020). Fully discrete positivity-preserving and energy-dissipating schemes for aggregation-diffusion equations with a gradient-flow structure. Commun. Math. Sci., 18(5), 1259-1303.
[3] Carrillo, J.A., Colombi, A., & Scianna, M. (2018). Adhesion and volume constraints via nonlocal interactions determine cell organisation and migration profiles. J. Theoret. Biol., 445, 75-91.
[4] Carrillo, J.A., Craig, K., Wang, L., & Wei, C. (2021). Primal dual methods for Wasserstein gradient flows. Found. Comput. Math., 22, 389-443.
[5] Carrillo, J.A., Craig, K., & Yao, Y. (2019). Aggregation-diffusion equations: Dynamics, asymptotics, and singular limits. In N. Bellomo, P. Degond, & E. Tadmor (Eds.), Active particles, Volume 2: Advances in theory, models, and applications (pp. 65-108). Modeling and simulation in science, engineering and technology. Cham, Switzerland: Birkhäuser.
[6] Carrillo, J.A., Hittmeir, S., Volzone, B., & Yao, Y. (2019). Nonlinear aggregation-diffusion equations: Radial symmetry and long time asymptotics. Invent. Math., 218(2), 889-977.
[7] Carrillo, J.A., Murakawa, H., Sato, M., Togashi, H., & Trush, O. (2019). A population dynamics model of cell-cell adhesion incorporating population pressure and density saturation. J. Theoret. Biol., 474, 14-24.
[8] Falcó, C., Baker, R.E., & Carrillo, J.A. (2024). A local continuum model of cell-cell adhesion. SIAM J. Appl. Math., 84(3), S17-S42.
[9] Gómez-Castro, D. (2024). Beginner’s guide to aggregation-diffusion equations. SeMA J.
[10] Jordan, R., Kinderlehrer, D., & Otto, F. (1998). The variational formulation of the Fokker-Planck equation. SIAM J. Math. Anal., 29(1), 1-17.
[11] Katsunuma, S., Honda, H., Shinoda, T., Ishimoto, Y., Miyata, T., Kiyonari, H., … Togashi, H. (2016). Synergistic action of nectins and cadherins generates the mosaic cellular pattern of the olfactory epithelium. J. Cell Biol., 212(5), 561-575.
[12] Keller, E.F., & Segel, L.A. (1970). Initiation of slide mold aggregation viewed as an instability. J. Theoret. Biol., 26(3), 399-415.
[13] Steinberg, M.S. (1962). On the mechanism of tissue reconstruction by dissociated cells. I. Population kinetics, differential adhesiveness, and the absence of directed migration. Proc. Natl. Acad. Sci., 48(9), 1577-1582.
[14] Vázquez, J.L. (2006). The porous medium equation: Mathematical theory. Oxford, U.K.: Oxford University Press.
[15] Villani, C. (2003). Topics in optimal transportation. In Graduate studies in mathematics (Vol. 58). Providence, RI: American Mathematical Society.