Particle Interactions Mediated by Dynamical Networks: Assessment of Macroscopic Descriptions

We provide a numerical study of the macroscopic model of Barré et al. (Multiscale Model Simul, 2017, to appear) derived from an agent-based model for a system of particles interacting through a dynamical network of links. Assuming that the network remodeling process is very fast, the macroscopic model takes the form of a single aggregation–diffusion equation for the density of particles. The theoretical study of the macroscopic model gives precise criteria for the phase transitions of the steady states, and in the one-dimensional case, we show numerically that the stationary solutions of the microscopic model undergo the same phase transitions and bifurcation types as the macroscopic model. In the two-dimensional case, we show that the numerical simulations of the macroscopic model are in excellent agreement with the predicted theoretical values. This study provides a partial validation of the formal derivation of the macroscopic model from a microscopic formulation and shows that the former is a consistent approximation of an underlying particle dynamics, making it a powerful tool for the modeling of dynamical networks at a large scale.


Introduction
Complex networks are of significant interest in many fields of life and social sciences. These systems are composed of a large number of agents interacting through local interactions, and self-organizing to reach large-scale functional structures. Examples of systems involving highly dynamical networks include neural networks, biological fiber networks such as connective tissues, vascular or neural networks, ant trails, polymers, economic interactions. (Boissard et al. 2013;Mogilner and Edelstein-Keshet 2007;DiDonna and Levine 2006;Broedersz et al. 2010). These networks often offer great plasticity by their ability to break and reform connections, giving to the system the ability to change shape and adapt to different situations (Boissard et al. 2013;Chaudury et al. 2007). For example, the biochemical reactions in a cell involve proteins-DNA, RNA, gene promotors linking/unlinking to create/break large structures-complex of molecules (Kupiec 1997). Because of their paramount importance in biological functions or social organizations, understanding the properties of such complex systems is of great interest. However, they are challenging to model due to the large amount of components and interactions (chemical, biological, social, etc). Due to their simplicity and flexibility, individual-based models are a natural framework to study complex systems. They describe the behavior of each agent and its interaction with the surrounding agents over time, offering a description of the system at the microscopic scale (see, e.g., Barré et al. 2017;Boissard et al. 2013;Degond et al. 2016). However, these models are computationally expensive and are not suited for the study of large systems. To study the systems at a macroscopic scale, mean-field or continuous models are often preferred. These last models describe the evolution in time of averaged quantities such as agent density and mean orientation. As a drawback, these last models lose the information at the individual level. In order to overcome this weakness of the continuous models, a possible route is to derive a macroscopic model from an agent-based formulation and to compare the obtained systems, as was done in, e.g., Barré et al. (2017), Boissard et al. (2013), Degond et al. (2016) for particle interactions mediated by dynamical networks.
A first step in this direction has been made in Barré et al. (2017), following the earlier work . In this work, the derivation of a macroscopic model for particles interacting through a dynamical network of links is performed. The microscopic model describes the evolution in time of point particles which interact with their close neighbors via local cross-links modeled by springs that are randomly created and destructed. In the mean-field limit, assuming large number of particles and links as well as propagation of chaos, the corresponding kinetic system consists of two equations: for the individual particle distribution function and for the link densities. The link density distribution provides a statistical description of the network connectivity which turns out to be quite flexible and easily generalizable to other types of complex networks.
In the large-scale limit and in the regime where link creation/destruction frequency is very large, it was shown in Barré et al. (2017), following Degond et al. (2016), that the link density distribution becomes a local function of the particle distribution density. The latter evolves on the slow time scale through an aggregation-diffusion equation. Such equations are encountered in many physical systems featuring collective behavior of animals, chemotaxis models, etc. (Topaz et al. 2006;Blanchet et al. 2006;Carrillo et al. 2010;Golestanian 2012;Kolokolnikov et al. 2013) and references therein. The difference between this macroscopic model and the aggregation-diffusion equations studied in the literature Carrillo et al. (2003), Topaz et al. (2006), Bertozzi et al. (2009) lies in the fact that the interaction potential has compact support. As a result, this model has a rich behavior such as metastability in the case of the whole space (Burger et al. 2014;Evers and Kolokolnikov 2016) and exhibits phase transitions in the periodic setting as functions of the diffusion coefficient, the interaction range of the potential, and the links equilibrium length (Barré et al. 2017). By performing the weakly nonlinear stability analysis of the spatially homogeneous steady states, it is possible to characterize the type of bifurcations appearing at the instability onset (Barré et al. 2017). We refer to Barbaro and Degond (2014), Chayes and Panferov (2010), Degond et al. (2015), Barbaro et al. (2016) for related collective dynamics problems showing phase transitions.
If numerous macroscopic models for dynamical networks have been proposed in the literature, most of them are based on phenomenological considerations and very few have been linked to an agent-based dynamics. On the contrary, the macroscopic model proposed in Barré et al. (2017) and its precursor  have been derived via a formal mean-field limit from an underlying particle dynamics (see also Degond et al. 2014). However, because the derivation performed in Barré et al. (2017) is still formal, its numerical validation as the limit of the microscopic model as well as the persistence of the phase transitions at the micro-and macroscopic level as predicted by the weakly nonlinear analysis in Barré et al. (2017) needs to be assessed. This is the goal of the present work.
More precisely, we show that the macroscopic model indeed provides a consistent approximation of the underlying agent-based model for dynamical networks, by confronting numerical simulations of both the micro-and macromodels. Moreover, we numerically check that the microscopic system undergoes in one-dimensional a phase transition depicted by the values obtained for the limiting macroscopic aggregation-diffusion equation. Furthermore, we numerically validate the weakly nonlinear analysis in Barré et al. (2017) for the type of bifurcation in the twodimensional setting, where simulations for the microscopic model are prohibitively expensive.
In summary, the main contributions of this work are as follows: (i) It provides a numerical validation of the macroscopic model in 1D as its derivation from the microscopic one in Barré et al. (2017) was only formal. It justifies its further use in 2D where the microscopic model is too computationally intensive. (ii) It also provides the experimental validation of the formal bifurcation analysis for the macroscopic model performed in Barré et al. (2017). In particular, it confirms that the two types of bifurcations revealed in Barré et al. (2017)-subcritical and supercritical-do actually occur. (iii) Finally, it shows that this bifurcation structure is indeed relevant for the microscopic model, for which no theoretical analysis exists to date.
The paper is organized as follows. In Sect. 2, we present the microscopic model and sketch the derivation of the kinetic and macroscopic models from the agent-based formulation. In Sect. 3, we focus on the one-dimensional case: We first summarize the theoretical results on the stability of homogeneous steady states of the macroscopic model from Barré et al. (2017) and show that both the macroscopic and microscopic simulations are in good agreement with the theoretical predictions made by nonlinear analysis of the macroscopic model. We then compare the profiles of the steady states between the microscopic and macroscopic simulations and show that the two formulations are in very good agreement, also in terms of phase transitions. Finally, in Sect. 4 we provide a numerical study of the two-dimensional case for the macroscopic model. The two-dimensional numerical simulations on the macroscopic model are able to numerically capture the subcritical and supercritical transitions as predicted theoretically. Because of the computational cost of the microscopic model, the macroscopic model is not only very competitive and efficient in order to detect phase transitions, but also it is almost the only feasible choice showing the main advantage of the limiting kinetic procedure.

Microscopic Model
The two-dimensional microscopic model features N particles located at points X i ∈ , i ∈ [1, N ] linking/unlinking-dynamically in time-to their neighbors which are located in a ball of radius R from their center. The link creation and suppression are supposed to follow Poisson processes in time, of frequencies ν N f and ν N d , respectively (see Fig. 1).
Each link is supposed to act as a spring by generating a pairwise potential where κ is the intensity of the spring force and the equilibrium length of the spring. We define the total energy of the system W related to the maintenance of the links: where i(k)and j (k) denote the indexes of particles connected by the link k. Particle motion between two linking/unlinking events is then supposed to occur in the steepest descent direction to this energy, in the so-called overdamped regime:

Kinetic Model
To perform the mean-field limit, following Barré et al. (2017) and Degond et al. (2016), we define the one-particle distribution of the N particles, f N (x, t), and the link distribution of the K links, g K (x 1 , x 2 , t). Postulating the existence of the following limits: the kinetic system reads: where we have postulated that the distribution of pairs of particles reduces to f (x 1 , t) f (x 2 , t), and We refer the reader to Barré et al. (2017) for details on the mean-field limit. In the equation for the limit distribution of particles (4), the first term on the righthand side is a linear diffusion term which is an effect of the random motion of the particles on the microscopic level. The second term is the attractive-repulsive part due to a spring-like force between the particles that are linked. Its counterpart appears in the equation for the limit distribution of links (5). This equation has also a diffusion part and the production term (the last two terms on the right-hand side) which is due linking processes taking place between the particles that are not yet connected, and unlinking processes breaking the existing links.

Scaling and Macroscopic Model
In this paper, the space and time scales are chosen such that μ = 1 and the variables are scaled such that: The spring force κ is supposed to be small, i.eκ = ε −1 κ, the noise D is supposed to be of order 1, and the typical spring length and particle detection distance R are supposed to scale as the space variable, i.e.,˜ = ε 1/2 ,R = ε 1/2 R. Finally, the main scaling assumption is to consider that the processes of linking and unlinking are very fast, i.e., ν f = ε 2 ν f ,ν d = ε 2 ν d . In the example of cell dynamics, mentioned in Introduction, see Kupiec (1997), the frequency of linking/unlinking depends on the size of the molecule, and it is a very fast process (order of seconds) while the macroevolution of the cell such as growth of the cell is much slower (order of minutes). For the sake of simplicity, we will consider in this paper thatν f ν d = 1, andκ = 2. For such a scaling, it is shown in Barré et al. (2017) that in the limit ε → 0, if we suppose ( f ε , g ε ) → ε→0 ( f, g), then: for some compactly supported potential V such that: In this paper, we take κ = 2 in (1); hence, V has the form: It is worthy to mention that in the context of cell dynamics several authors have already advocated for non-local terms of the form in (6) to model cell adhesion, see, for instance, Domschke et al. (2014), Painter et al. (2015), and the references therein. Therefore, we can reinterpret this fast linking/unlinking as a way of modeling cell adhesion at the microscopic level.
In the following, we aim to study theoretically and numerically both the macroscopic model given by Eq. (6) and the corresponding microscopic formulation given by Eq. (3) and rescaled with the scaling introduced in this section. We first focus on the one-dimensional case and we show that the numerical solutions behave as theoretically predicted, and that we obtain-numerically-a very good agreement between the micro-and macroformulations.

Theoretical Results
In this section, we apply the results of Barré et al. (2017) to the one-dimensional periodic domain [−L , L], to study the stability of stationary solutions of the macroscopic model given by Eqs. (6a) and (7) with R < L.

Identification of the Stability Region
We first linearize equation (6a) around the constant steady state ρ * = 1 2L , so that the total mass is equal to 1, we denote the perturbation by ρ, so we have f = ρ * + ρ, that satisfies where V is given by (7). We will further decompose f into its Fourier modes and the Fourier transform is given bŷ Applying the Fourier transform to (8), a straightforward computation gives where the Fourier modes of the potential V are given bŷ Here, we denoted Therefore, the stability of the constant steady state will be ensured if the coefficient in front ofρ k on the r.h.s. of (9) has a non-positive real part for k = 1. Indeed, as observed in Barré et al. (2017), this condition implies that all the other modes are then also stable. This condition is related to the H-stable/catastrophic behavior of interaction potentials that characterizes the existence of global minimizers of the total potential energy as recently shown in Cañizo et al. (2015), Simione et al. (2015).

Characterization of the Bifurcation Type
As shown in Barré et al. (2017), it is possible to distinguish two types of bifurcation as functions of the model parameters. Indeed, if we define: we have the following proposition (see Barré et al. 2017): Proposition 1 Assume that λ > 0 and λ k < 0, ∀k = ±1. Then, • if 2V 2 −V −1 > 0, the steady state exhibits a supercritical bifurcation; • if 2V 2 −V −1 < 0, the steady state exhibits a subcritical bifurcation.
Note that the above criterion only involves the potential, but does not involve the parameter D, and it only restricts the values of α or .

Choice of Numerical Parameters
In the linearized equation (8), there are four parameters that may vary: D, , R, and L. In this part of the paper, we focus on the case where the potential is of comparable range R to the size of the domain L, and fix the value of the following parameters: L = 3 and R = 0.75; therefore, z 1 = π 4 . Using (9) and the discussion from the end of Sect. 3.1.1, we can identify the region where the constant steady state is unstable. ComputingV 1 < 0 and V 1 < −D, respectively, leads to the following restriction for two remaining parameters of the system and D: which allows to approximate the instability region for this particular case as D < D( ) = 0.1781(0.4948 − ). We also introduce a notation c = Rα c , which in this case gives c = 0.4948. The parameter c denotes the value of above which the constant steady state is always stable independently of the value of the parameter D. Using (10) and Proposition 1, we check that the bifurcation changes its character ≈ 0.4530. Recall that our criterion did not involve the parameter D; therefore, the bifurcation is supercritical if only ∈ ( * , c ) ≈ (0.4530, 0.4948), and subcritical if ∈ (0, l * ) ≈ (0, 0.4530). The value of parameter D corresponding to the instability threshold for l = l * ≈ 0.4530 is denoted by D * , and it is equal to 0.0074. All of these parameters are presented in Fig. 2.
Remark 1 Choosing R comparable to L allows us to observe two types of bifurcation: continuous and discontinuous one. It was observed in Barré et al. (2017) that taking R L would cause that for most values of the bifurcation would be subcritical (discontinuous). This effect is captured in Fig. 8, for the two-dimensional case.

Macroscopic Model
We now make use of the numerical scheme developed in Carrillo et al. (2015) to analyze the macroscopic equation (6a) with the potential (7) in the unstable regime. The choice of the numerical scheme is due to its free energy decreasing property for equations enjoying a gradient flow structure such as (6a). Keeping this property of gradient flows is of paramount importance in order to compute the right stationary states in the long time asymptotics. In fact, under a suitable CFL condition the scheme is positivity preserving and well balanced, i.e., stationary states are preserved exactly by the scheme.
To check the correctness of the criterion from Proposition 1, we consider two cases corresponding to two different types of bifurcation, as depicted in Fig. 2: In order to trace the influence of the diffusion on the type of bifurcation, for fixed 1 , 2 , we will be looking for the values of diffusion coefficients D 1,λ , D 2,λ such that Recall that according to Barré et al. (2017), the parameter λ defined in (11a) measures the distance from the instability threshold. We will use this information to determine the values of parameters D 1,λ = D 1,λ (λ) and D 2,λ = D 2,λ (λ) computed from (11a). We consider 14 different values for subcritical and supercritical case, as specified in Table 1. Moreover, in Barré et al. (2017) the authors proved that the perturbation ρ(t) of the constant steady state satisfies the following equation and Equation (13) means that for the supercritical bifurcation we can observe a saturation. This means that before stabilizing A(t) first grows exponentially until the r.h.s. of (13) is equal to zero, i.e., for Using this information to estimate the r.h.s. of (12), we obtain that This condition gives us the upper estimate for the amplitude of perturbation ρ when the steady state is achieved, that is, after the saturation. The derivation of Proposition 1 in Barré et al. (2017) assumes sufficiently small perturbation of the steady state. Therefore, the initial data for our numerical simulations should be least smaller than the value of |A| corresponding to the saturation level. It turns out that |A| computed in (14) is always less than √ λ, so the size of initial perturbation of the steady state should be also taken in this regime. If we choose the initial data for the numerical simulations of the supercritical case in this regime, we should see a continuous decay of the saturated amplitude of perturbation to 0, as λ decreases. We will perturb the initial data for the subcritical case similarly, showing that even though the smallness restriction is respected, the saturated amplitude of perturbation is a discontinuous function of λ.
In what follows, we perturb the constant initial condition by the first Fourier mode: In the numerical simulations, we consider the case δ = 0.01. In order to distinguish between the homogeneous steady states (corresponding to the stable regime) and the aggregated steady states (corresponding to the unstable regimes), we compute the following quantifier Q on the density profiles of the numerical solutions: where where T max corresponds to the formation of the steady state. Note that (i) if the steady state is homogeneous in space then Q = 0 and (ii) if f is a symmetric function with respect to x, then Q = c 1 .
To estimate T max , we use the following criterion. From the theory (Carrillo et al. 2003), we know that steady states are positive everywhere and the quantity ξ = D log + V * is equal to some constant C. We then compute the distance of ξ from its mean value: The steady state is achieved if ξ * is sufficiently close to 0, and in our numerical scheme we continue the computations until t = T max for which ξ * (T max ) < 10 −7 .
The computed values are presented in Tables 5 and 6 in "Appendix 2". In Fig. 3, we show the values of the order parameter Q as a function of the noise intensity D for both types of bifurcation. As shown in Fig. 3, the quantifier Q indeed undergoes a discontinuous transition around D = 0.0347 for = 0.3 (subcritical case, Fig. 3a) and a continuous transition around D = 0.004 for = 0.4725 (supercritical case, Fig. 3b). These results show that the numerical solutions are in very good agreement with the theoretical predictions.
In order to check the accuracy of our prediction of the value of T max , we show in Fig. 4 the graph of ξ * (t) for several values of D in the supercritical and the subcritical cases (see Table 1). As shown in Fig. 4, we observe a very sharp change of ξ * for the subcritical bifurcation and much smoother one for the supercritical case. The amplitude change of ξ * is also a good indication of the type of bifurcation. As for the order parameter, we see that for the subcritical bifurcation it is on similar level (Fig.  4a) for all values of D, while for the supercritical bifurcation it decays to 0 (Fig. 4b).
We will use this observation to analyze the results of the two-dimensional simulations later on. Finally, we can also check how the theoretical prediction of the size of perturbation from (15) is confirmed by our numerical results. For this purpose, we compute the maximum of the perturbation once the steady state is achieved: for all the points of supercritical bifurcation. The results are presented in Fig. 5 and in Table 2.
We now aim at performing the same stability analysis on the microscopic model from Sect. 2.1-the starting point of the derivation of the macroscopic model.

Microscopic Model
Here, we aim at performing simulations of the microscopic model from Section 2.1, rescaled with the scaling from Section 2.3. After rescaling and if we consider an explicit Euler scheme in time (see "Appendix 1"), we can show that Eq. (3) between time steps t n and t n + t n reads (in non-dimensionalized variables): where W is defined by (2) and N (0, 1) is the normal distribution with mean 0 and standard deviation 1. Between two time steps, new links are created between close enough pairs of particles that are not already linked with probability P f = 1−e ν f t n /((N −1)ε 2 ) and the existing links disappear with probability P d = 1 − e −ν d t n /((N −1)ε 2 ) . Therefore, the rescaled version of the microscopic model features a very fast link creation/destruction rate, as the linking and unlinking frequencies are supposed to be of order 1/ε 2 , for small ε. Note also that to capture the right time scale, the time step t must be decreased with ε, which makes the microscopic model computationally costly for small values of ε. For computation time reasons, we also consider the limiting case ε = 0 of the microscopic model; we can show that it reads: where  Table 3.
As for the macroscopic model, the order of the particle system at equilibrium is measured by the quantifier Q defined by Eq. (16) where f i = f (x i ) and N i are, respectively, the density and the number of agents whose centers belong to the interval [−L + (i − 1)dx, −L + i dx]. Following the analysis of the macroscopic model, we explore the same two cases: 1 = 0.4725, D 1 = 0.0040, and 2 = 0.3, D 2 = 0.0347 to check whether they correspond to the super and subcritical bifurcations, respectively. In Fig. 6, we show the values of Q plotted as functions of the noise intensity D computed from the simulations of the scaled microscopic model (17)  These results show that the scaled microscopic model has the same properties as the macroscopic one, and that the values of the parameters ( , D) which correspond to a bifurcation in the steady states tend, as ε → 0, to the ones predicted by the analysis of the macroscopic model. Indeed for the limiting case "ε = 0" of the microscopic model, we obtain a very good agreement between the micro-and macroformulations, showing that the microscopic model behaves as predicted by the analysis of the macroscopic model. It is noteworthy that the small differences observed in the values of the transitional D (subcritical case, Fig. 6a) can be due to the fact that we use a finite number of N = 1000 particles for the microscopic simulations, whereas the macroscopic model is in the limit N → ∞. However, these differences are very small when we consider the limit case ε = 0 for the microscopic model. Indeed, determining visually the transitional D in the microscopic simulations with neglecting the slight increase appearing after (see Fig. 6b), the relative error between the microscopic and macroscopic transitional D, |D mic −D mac | D mac is 7% for = 0.3, and 5% for = 0.4725. In order to give a more quantitative analysis on the influence of the number of particles, we show in Fig. 7 the values of Q plotted as functions of the noise intensity D, for the macroscopic model (dashed red curves), and for the microscopic model with "ε = 0" [Eq. (18)] and different number of particles N : N = 500 (green curves), N = 1000 (blue curves), and N = 2000 (red curves). Figure 7a shows the case = 0.3 (subcritical bifurcation) and (b) the case = 0.4725 (supercritical bifurcation). As depicted in Fig. 7, as the number of particle increases, the value of the critical noise intensities D c get closer to the ones predicted by the macroscopic model for both the subcritical and supercritical transitions. Moreover, the values of Q corresponding to space homogeneous equilibria (for D > D c in both cases) get closer to zero as N increases, and its variations after the transitional D observed for N = 500 or N = 1000 get negligible for N = 2000. Altogether, these results show that the microscopic model is a good approximation of the macroscopic dynamics when considering a large number of particles and a small value of ε. It is noteworthy that the simulations of the microscopic model become very time-consuming when considering N = 2000 particles, and we refer the reader to "Computational Aspects of the Micro-and-Macroscopic Models" section of Appendix 1 for a detailed analysis of the computational time.
We now aim at comparing the profiles of the solutions between the microscopic and macroscopic models, to numerically validate the derivation of the macroscopic model from the microscopic dynamics.

Comparison of the Density Profiles in the Microscopic and Macroscopic Models
Here, we aim at comparing the profiles of the particle densities of the microscopic model with the ones of the macroscopic model as functions of time. As shown in the previous section, for ε small enough, we recover the bifurcation and bifurcation types observed from the macroscopic model with the microscopic formulation, with very good quantitative agreement when considering the limiting microscopic model (18) with "ε = 0." The simulations of the microscopic model are very time-consuming for small values of ε, because we are obliged to consider very small time steps (see "Computational Aspects of the Micro-and-Macroscopic Models" section of Appendix 1). Here, due to computational time constraints, we therefore compare the results of the macroscopic model (6) with ε = 0 for which the time step can be taken much larger and independent of ε.
In order to have the same initial condition for both the microscopic and macroscopic models, we initially choose the particle positions for both models such that: We send the reader to "Appendix 1" for the numerical method used to set the initial conditions of the microscopic model. Because of the stochastic nature of the model, the microscopic model does not preserve the symmetry of the solution, contrary to the macroscopic model (where noise results in a deterministic diffusion term). To enable the comparison between the macroscopic and microscopic models, we therefore recenter the periodic domain of the microscopic model such that the center of mass of the particles is located at x = 0 (center of the domain). To this aim, given the set of particles X j , j = 1 . . . N , we reposition all the particles at pointsX j , j = 1 . . . N such that: where X m is the center of mass computed on a periodic domain: Finally, in order to decrease the noise in the data of the microscopic simulations due to the random processes, the density of particles is computed on a set of several simulations of the microscopic model. In Fig. 8, we show the density distributions of the macroscopic model (continuous lines) and of the microscopic one with "ε = 0" (circle markers) at different times, for = 0.4725 and = 0.3, respectively. For each value of , we consider two values for the noise intensity D: For = 0.4725, we study the cases D = 0.003 and D = 0.0003, and for = 0.3, we choose D = 0.0338 and D = 0.0034. Note that all these values are in the unstable regime.
As shown in Fig. 8, we obtain a very good agreement between the solutions of the macroscopic model and of the microscopic one with "ε = 0." Close to the transitional D (Fig. 8a, b), the particle density converges in time toward a Gaussian-like distribution for both the microscopic and macroscopic models. Note that the microscopic simulations seem to converge in time toward the steady state faster than the macroscopic model (compare the orange curves on the top panels). This change in speed can be due to the fact that the microscopic model features finite number of particles while the macroscopic model is obtained in the limit of infinite number of particles. Therefore, in the macroscopic setting, each particle interacts with many more particles than in the microscopic model, which could result in a delay in the aggregation process.
When far from the transitional D in the unstable regime (Fig. 8c, d), one can observe the production of several bumps in the steady state of the particle density. The production of several particle clusters in these regimes shows that the noise triggers particle aggregation. For small noise intensity, local particle aggregates are formed which fail to detect neighboring aggregates. As a result, one can observe several clusters in the steady state, for small enough noise intensities. These bumps are observed for both the microscopic and macroscopic models, showing again a good agreement between the two dynamics. In the next section, we present a numerical study of the macroscopic model in the two-dimensional case. As mentioned previously, the microscopic model is in very good agreement with the macroscopic dynamics for small values of ε as in the onedimensional case. Its simulations are, however, very time-consuming, due to the need of very small time steps (see "Computational Aspects of the Micro-and-Macroscopic Models" section of Appendix 1). As a result, the microscopic model is not suited for the study of very large systems such as the ones considered in the two-dimensional case. We therefore provide a numerical two-dimensional study using the macroscopic model only.
[−L , L] × [−L , L], since the rectangular case can be, in agreement with the analysis in Barré et al. (2017), reduced to the one-dimensional case studied above.
The starting point for the phase transition analysis is the linearized equation in which the spatially homogeneous distribution ρ * is now equal to 1 (2L) 2 . Applying the Fourier transform to this equation, we obtain and we denote λ ±1,0 = λ 0,±1 = λ. The Fourier transform of the potential V is given bŷ where we denoted and J i are Bessel function of order i and H i are the Struve functions defined by Again, fixing the ratio R L ≤ 1, the relation between D, l, and R for the phase transition can be read from the condition λ = 0, which yields which due to (19) gives The relevant criterion for the type of bifurcation in the two-dimensional case then reads: (20), and such that λ k 1 ,k 2 remains negative for all k 1 , k 2 such that |k 1 | + |k 2 | > 1, let

then,
• if c < d, the constant steady state exhibits a supercritical bifurcation, and • if c > d, the constant steady state exhibits a subcritical bifurcation.
Note that in the two-dimensional case, the bifurcation criterion involves also parameter D. On the other hand, the instability threshold D is given as a function of α and can be calculated using (20).

Numerical Results
We first compute the approximate instability regime for the following three cases: 1. For R/L = 1, z 1,0 = π , the constant steady state is unstable for R < α c := 0.6620 and D < 0.2334R 2 α c − R .
This computation confirms the theoretical prediction from Barré et al. (2017) that the smaller R L is, the more likely it is that the bifurcation is of the subcritical type. The different types of bifurcation happen only for R L close to 1; otherwise, the bifurcation is always subcritical. To understand this behavior, we compute V * =V 1,0 (2V 2,0 −V −1,0 ) D +V 1,1 + 4 D +V 2,0 V 1,0V1,1 .
From Proposition 2, it follows that if V * > 0 the bifurcation is subcritical; otherwise, it is supercritical. We depict the function V * (α), where α = R for different values of R L ∈ [0.9, 1] in Fig. 9a. We see that decreasing the ratio R L causes that more and more of the graph of V * (α) lies above 0. This means that for most of the values of α ∈ [0, α c ] the bifurcation is subcritical.
We will study all of the five cases from Fig. 9b corresponding to different values of R L , but the same value of α = R = 0.3. The theoretical prediction is that the first two cases R L = 1 and R L = 0.975 correspond to a supercritical (continuous) bifurcation while the cases 3-5 correspond to the subcritical (discontinuous) bifurcation.
We perturb the constant initial data as in the one-dimensional case; namely, we take with δ = 0.01, and similarly to the one-dimensional case we compute the value of the order parameter Q where we used the empirical observation that the steady state is always symmetric with respect to (x, y) = (0, 0). For the stopping time criterion, we take the same as in one-dimensional case, namely ξ (T max ) < 10 −7 . In Fig. 10, we show the values of the order parameter Q as function of the noise intensity D for both types of bifurcation for cases 1 and 5, based on Tables 7 and 11 from "Appendix 2". As shown in Fig. 10, we indeed obtain a supercritical (continuous) transition in the values of Q as function of the noise D in case 1 (Fig. 10b), while the transition is discontinuous (subcritical) in case 5 (Fig. 10a). These results therefore show that the numerical results are in good agreement with the theoretical predictions and provide a validation of the numerical approximation and simulations of the macroscopic model.
The difference between the bifurcation types is also reflected in the amplitude of the steady state. For both types of bifurcation, i.e., for cases 1 and cases 5 we plot the  Table 4 final steady states in Fig. 11. The density profile for the supercritical bifurcation ( Fig.  11b) is much lower and rounded than the one for the subcritical bifurcation (Fig. 11a). Moreover, as in the one-dimensional case, we can check that the different bifurcation diagrams correspond to different shapes of ξ * . In Figs. 12 and 13, we present the graphs of ξ * (t) for all five cases depicted in Fig. 9b. For each of the cases, we present the graph of ξ * (t) for five different values of diffusion parameter D as specified in Table 4.
As shown in Fig. 12, the graph of ξ * (t) undergoes smooth changes for the different values of the noise D, highlighting a bifurcation of supercritical type. Figure 13 shows that ξ * (t) undergoes sharp changes for the different values of the noise D, highlighting the subcritical type of bifurcation, as predicted by the theoretical analysis of the macroscopic model in the two-dimensional case. Close to the transition zone (case 3, R L = 0.95, as shown in Fig. 13a), the changes in ξ * are smoother than for smaller values of R L (Fig. 13b, c), but the transition is still subcritical as can be confirmed by the values of order parameter Q given in Table 9.

Conclusion
In this paper, we have provided a numerical study of a macroscopic model derived from an agent-based formulation for particles interacting through a dynamical network of links. In the one-dimensional case, we were first able to recover numerically the subcritical and supercritical transitions undergone by the steady states of the macroscopic model, in the regime predicted by the theoretical nonlinear analysis of the continuous model. Moreover, the numerical simulations of the rescaled microscopic model revealed the same bifurcations and bifurcation types as obtained with the macroscopic model, with very good precision as ε goes to zero in the microscopic setting. Finally, when considering the limiting case "ε = 0" in the microscopic model, we obtained a very good agreement between the profiles of the solutions of the micro-and macromodels. It is noteworthy that both models also feature the same dynamics in time, with a slight delay in the macroscopic simulations compared to the microscopic dynamics. This delay may be due to the fact that the microscopic simulations are performed with a finite number of particles while the macroscopic model is in the limit of infinite number of individuals. However, as for very small values of ε the simulations of the microscopic dynamics are very time-consuming, we were not able to extend the numerical study to a higher number of particles. For the sake of completeness, we finally presented numerical simulations of the model in the two-dimensional case. For computational reasons, we were not able to perform two-dimensional simulations of the microscopic model, and we chose to focus on the macroscopic model. In the two-dimensional case, we were once again able to numerically recover supercritical and subcritical transitions in the steady states, as function of the noise intensity D, in the same regime as predicted by the theoretical analysis of the macroscopic model. These results validate the theoretical analysis, the numerical method, and the simulations developed for the macroscopic model.
By providing a numerical comparison between the micro-and macrodynamics, this study shows that the macroscopic model considered in this paper is indeed a relevant tool to model particles interacting through a dynamical network of links. As a main advantage compared to the microscopic formulation, the macroscopic model enables to explore large systems with low computational cost (such as two-dimensional studies) and is therefore believed to be a powerful tool to study network systems on the large scale. Direct perspectives of these works include the derivation of the macroscopic model in a regime of non-instantaneous linking-unlinking of particles. The hope is to understand deeper how the local forces generated by the links are expressed at the macroscopic level. The model could be improved by taking into account other phenomena such as external forces and particle creation/destruction. Finally, rigorously proving the derivation of the macroscopic model from the particle dynamics will be the subject of the future research.
ticles is large enough so that we can set ν N f = ν f N −1 and ν N d = ν d . Finally, as explained in the main text, we initially choose the particle positions for both models such that: For the microscopic model, the initial positions of the particles are set such that given a random position X ∈ [−L , L] and a random number η 2 ∈ [−1, 1], we let the probability of creating a new particle at position X + η 2 x − x 2 be x 2L + δ(λ) x cos π X L .

A.2 Computational Aspects of the Micro-and-Macroscopic Models
The numerical scheme of the microscopic model is implemented in FORTRAN90. Parallelization is only used when running several simulations for the comparison with the macroscopic model (see Fig. 8 in the unstable regime than in the stable one (see discussion bellow), and finally, the CFL-type condition of the microscopic model (see previous paragraph) links the time step to the maximal number of links per fiber, also directly linked to the clustering of the particles. Altogether, these effects lead to a CPU time of order O( 1 3 ), dependent on the scaling parameter as well as on the value of the statistical quantifier Q which measures the "clustering" of the particles (compare Fig. 14 with Fig. 6). Note that in the limit case " = 0" (dashed lines in Fig. 14), the CPU time is much smaller than for = 0. Indeed, in this limit case there is no more linking dynamics as particles interact with all of their close neighbors. Therefore, in this case, the computational time does not depend on the linking and unlinking frequencies anymore, but still on the quantifier Q. Note that these results correspond to a fixed number of particles N = 1000, but that the CPU time also depends on this parameter.
In order to highlight the dependency of the CPU time on the clustering of particles and on the number of particles, we show in Fig. 15 Figure 15a shows the case N = 500 particles, and (b) and (c) correspond to N = 1000 and N = 2000, respectively. As shown in Fig. 15, the CPU time strongly depends on the amount of aggregation of the particles in the simulation. At the beginning of the simulation, particles are distributed almost homogeneously in space, leading to a small CPU time. As particles start aggregating, one observes an increase in CPU time, by a factor ranging from 3 to 6 as the number of particles increases from 500 to 2000 (Fig. 15a-c). These figures clearly show the simulation time at which particles start aggregating, which depends on the value of D as well as on the number of particles and on the regime we consider: Indeed, particles aggregate faster in the subcritical regime than in the supercritical one. The fact that formation of clusters in the unstable regime leads to a larger CPU time than in the stable one is directly linked to the presence of particle aggregates, as each particle interacts with more neighbors in the clustered case than in the homogeneous one. Altogether, these numerical results show that the microscopic model becomes very time-consuming when considering small values of and in the unstable regime or even when considering large number of particles. Our simulations reached 40 days for = 1 12 and D = 0.03 in the subcritical case ( = 0.3), whereas the CPU times for the macroscopic model ranged from few minutes to few hours. These results highlight the need for a continuous formalism when studying systems with large number of particles.
Let us finally mention that the numerical scheme for the macroscopic model is amenable to parallelization since the numerical cost bottleneck is due to the convolution terms and this is not difficult to parallelize as it is scalable using load balancing between processors. The communication costs to compute convolutions are highly reduced due to the compact support of the potential. The other only communication costs will come at the time of evaluating the new values of the density that needs fluxes from neighboring cells. However, since we only use two numerical values of the flux, the communication costs will be rather small. These aspects from the parallel programming perspective have been treated in equations with similar issues such as kinetic equations, see, for instance, [28]. The strategy there can be applied directly to our situation leads to speedups for two-dimensional computations of the macroscopic model.

B.1 Numerical Results for the One-Dimensional Case
It this section, we provide the data for the one-dimensional simulations of the supercritical bifurcation (Table 5) and of the subcritical bifurcation (Table 6) for the macroscopic model.

B.2 Numerical Results for the Two-Dimensional Case
It this section, we provide the data for the two-dimensional simulations of the supercritical bifurcation (Tables 7, 8) and of the subcritical bifurcation (Tables 9, 10, 11) for the macroscopic model.