# Spontaneous repulsion in the reaction on coupled networks

## Abstract

We study the transient dynamics of an process on a pair of randomly coupled networks, where reactants are initially separated. We find that, for sufficiently small fractions of cross-couplings, the concentration of (or ) particles decays linearly in a first stage and crosses over to a second linear decrease at a mixing time . By numerical and analytical arguments, we show that for symmetric and homogeneous structures where is the mean degree of both networks. Being this behavior in marked contrast with a purely diffusive process—where the mixing time would go simply like —we identify the logarithmic slowing down in to be the result of a novel spontaneous mechanism of repulsion between the reactants and due to the interactions taking place at the networks’ interface. We show numerically how this spontaneous repulsion effect depends on the topology of the underlying networks.

###### pacs:

89.75.Hc, 05.40.-a, 82.20.-w^{1}

Nearly one hundred years ago, Marian von Smoluchowski introduced a mathematical model to describe coagulation phenomena in terms of diffusion-controlled reaction processes Smo917 ().
Despite its apparent simplicity, the kinetics of this model was found to yield a wealth of intriguing phenomena, whose analysis have widely enriched our understanding of pattern formation in chemical compounds, biological systems Murray (); Kaneko () and elsewhere.

From the Statistical Mechanics perspective, reaction-diffusion (RD) processes represent a fertile groundwork where to analyze the emergence of spontaneous mechanism by starting from microscopic rules kampen1992 ().
Most studies in this direction aimed to unveil the effect that dynamical correlations and geometrical (or topological) constraints of the underlying structures have on the spatiotemporal evolution of the reactants’ concentrations ShlomoBen ().

The process, in particular, is known to exhibit anomalous kinetics on low dimensional and fractal geometries, where density fluctuations yield the formation of self-segregation domains composed of particles of the same type OZ (); Wilczek (); Redner (); Front ().
These phenomena result in a drastic slowing down in the rate of the reactions, forcing the system in a long-lived non-equilibrium state with a sub-linear decay in the density of the surviving particles.
Since this type of process grasps the essential kinetics featured by the spreading of pathogen-antipathogen agents Castellano2009 (); Vespignani (), or underlying the pattern-formation of diverse chemical reaction Galfi (); Ovchinnikov (); Taitelbaum (), the appearance of sub-diffusive dynamics may become inefficient for practical applications, which typically aim in fast mixing regimes. With this goal, it was later proved that the adoption of Lévy processes indeed washes out segregation phenomena, leading to super-diffusive dynamics Levy1996 ().

After this “classical” period, the study of RD processes has experienced a relevant boost with the inception of network theory as a new field for characterizing the structures underlying real-world complex systems Newman (); Reuvenbook ().
In fact, besides the focus on networks’ topological properties, a mainstream has been (and still is) to understand the interplay between their structure and the dynamics of processes taking place on them Sbrocca ().

Consistently with the scenario observed in other models Calda (); collective (), complex networks significantly influence the collective properties of RD processes pastorsatorras2015 ().
Numerical gallos1 (); gallos2 () and theoretical catanzaro () results have showed that the small-world property of these substrates mitigates the local fluctuations in the particles’ density, facilitating their reactions.
Notably, on scale-free (SF) networks (i.e. random graphs with connectivity distribution and ) the kinetics of the process exhibit jamming effects at early stages and then super-diffusive behaviors weber (), with the latter becoming stronger as the network heterogeneity increases.
Homogeneous structures—like random regular (RR), Erdős-Rényi (ER) or scale-rich (SR, i.e. and ) networks—result instead in a linear decay of the density, in accordance with the mean-field predictions gallosRev (); catanzaroRev ().

Though the portrait of RD processes on isolated structures is nowadays clear, not much is known regarding their dynamical behaviors on multilayer networks gomez (); garas ().
The influence that their mesoscopic organization has on the collective behaviors of processes acting on them is attracting significant interest havlin (); boccaletti (); Kivela (); ginestra (); garasbook (), and has already produced interesting results interdep1 (); interdep2 (); bornholdt17 ().
Increasing evidence, in fact, is showing that the existence of multiple layers, joined with the possibility of modeling different types of cross-system interactions, results in novel collective behaviors whose analysis is still in its infancy charoSA (); arenas17 (); micha17 (). Following this mainstream, we study here the dynamics on a pair of randomly coupled networks with initially separated reactants, where we find a novel spontaneous dynamical phenomenon.

Our results show that, for sufficiently small fractions of the cross-couplings between the layers, the concentration of both reactants decays as for a long transient, and then crosses over to a second linear regime where (with ) at a mixing time .
We interpret the initial transient () as the unmixed regime, where the reactions between and particles take place mainly at the boundaries between the networks (i.e. at the interconnected nodes).
At larger times (), the reactants penetrate more and more the two layers and start to react everywhere in the system.
After this fast mixing stage (, the remaining particles are uniformly distributed in the system and their dynamics is driven essentially by diffusion.
We find that the mixing time depends on the ratio , where is the average degree of the networks, according to the formula

(1) |

that we derive analytically for RR graphs and verified by extensive simulations on diverse synthetic networks.
Since Eq. (1) is in marked contrast with a purely diffusive process, where would simply scale as , we interpret the logarithmic factor as the reflection of a spontaneous repulsion mechanism between reactants due to the reactions taking place at the boundary between the networks.
We find that this mechanism becomes stronger with increasing heterogeneity of the underlying structures, in which case Eq. (1) holds only approximately.

The paper is organized as follows.
We derive analytically Eq. (1) for RR graphs, that we verify against extensive simulations on uncorrelated configuration model (UCM) networks.
After investigating the effects that the underlying topology has on and on the dynamical regimes observed, we give our conclusions.

### I. Analytic approach.

We consider two configuration model networks Bollobas (), composed by the same number of nodes and same structural properties, a condition that we will refer hereafter as the “symmetry” of the interconnected network.
Let us further assume the two layers are coupled by means of undirected interlinks, placed at random between a fraction of couples of nodes belonging to different layers (Fig. 1).
Two populations of and reactants are then randomly distributed with initially separated concentrations, so that all the particles of the same type are placed on the same layer.
For simplicity we assume that the initial concentrations of reactants are equal.
To track the evolution of each populations, let us denote by and the concentration of particles in network 1 and 2, respectively; similarly, let and be the concentration of particles in network 1 and 2, respectively.
Having assumed equal initial conditions, we have and .
Moreover, by symmetry, and for all .
It is worth to notice here that, whilst this symmetric condition certainly holds in the case of coupled layers with the same or mildly different topological features (say RR and ER layers, or two ER networks with different average degrees), it will in general require some adjustments for layers with different structures.

To further simplify the analysis, let us assume that the networks underlying the RD process have homogeneous topologies, so that the average degree of nodes is the only characteristic parameter of the structure. In this case, disregarding any dynamical effect due to the topological fluctuations, we assume that the overall behavior is captured by the average densities and . We can then describe the rate of change of the concentrations and which, by symmetry, is given by

(2) | ||||

(3) |

where the first two terms in both lines are due to diffusion, and the last term is due to reaction.
The effective diffusion rate is the probability at each node to move to the other network.

Subtracting Eq. (2) from Eq. (3), one obtains whose solution is .
Inserting this expression into, say Eq. (2), leads to

(4) |

which is a Riccati differential equation. Being the densities now decoupled, we can focus hereafter only on the solutions of one of them, dropping dumb labels.
Solving numerically Eq. (4) will give us the theoretical predictions for the evolution of the reactants’ concentrations.

Let us now notice that Eq. (4) is characterized by two distinct regimes: a reaction-dominated regime, where with solution which can be understood as the asymptotic behavior of the system, and a diffusion-limited regime where the exponentially decaying terms are dominant over the reaction one. By comparing these two regimes, we obtain an equation for the quasi-stationary behavior, namely

(5) |

Eq. (5) can be adopted in order to define the mixing time , that is the characteristic time it takes for the and reactants to mix and react like in a single network. To this aim, we make the assumption that the second summand on the right hand side of Eq. (5) is dominant, i.e. comment (). In this approximation, we get

(6) |

which, in combination with the solution of the reaction-limited regime, gives to the leading order Eq. (1).

### II. Simulations and results

To simulate the process on two randomly coupled networks, we have first generated two equal synthetic networks composed of nodes where the reactants will be initially distributed.
Once the networks are prepared, we randomly place the cross-couplings between them with probability by quenching the labels of interconnected nodes.
To avoid any dynamical effect due to degree-degree correlations, we have constructed two random networks according to the UCM.
Each network is hence assigned with a specific degree sequence having the desired connectivity distribution for the structure, with lower and structural cut-off given by , where is the minimum degree of each node UCM ().

Two population of reacting ( and ) species are then randomly distributed on the interconnected network with initially separated concentrations, meaning that all particles are placed on nodes of one layer, and all particles on nodes of the other.
Following our symmetric choice for the system, we fix the initial densities of reactants to be the same, so that .

Particles diffuse in the system by performing independent random walks, where hops are allowed only to nearest neighbor nodes.
Being interested in studying an process, we further assume that reactants of the same species do not interact with each other once they occupy the same site simultaneously, i.e. we adopt a bosonic version of the dynamics bosonic1 (); bosonic3 ().
A reaction occurs whenever an and a particle occupy the same site, in which case both reactants generate an inert species and are then removed from the system.
We monitor the time evolution of the concentrations of and particles, where the total time advances at each step as , being the number of particles currently present in the system.
Results are then averaged over a set of realizations, whose nominal cardinality is hereafter assumed to be , unless otherwise stated.

In Fig. 2 we present the inverse particle density as a function of time for decreasing values of in coupled RR (Fig. 2a), ER (Fig. 2b) and SF networks (Fig. 2c).
In all the cases, we find that for small enough values of the fractions of interconnected nodes, three distinct dynamical regimes exist.
One for short times (), where the particles diffuse inside their own layer and reactions occur mainly at the interface, one for intermediate times (), where particles start to cross and react in the opposite layer, and a third one for long times (), where eventually the survived reactants are well mixed and the coupled systems behave like a single network.
As shown in the inset of Fig. 2b, this kinetic pattern strongly depends on the choice of an initial separation of the reactants, and eventually disappears when the particles concentrations are initially mixed.
To demonstrate these results, we have compared in Fig. 2b the numerical solution obtained from the theory and given by Eq. (4) with the data collected from the simulation, in which case we observe an excellent agreement.

For homogeneous structures, we have also tested the effects that the average degree has on the mixing of the system.
In Fig. 2d, it is depicted again the inverse particle density as a function of time for ER networks, only this time with a fixed fraction and increasing values of the mean degree.
For these substrates, we find that the same repulsion mechanism is obtained by both increasing or decreasing , as one might expect since a particle hopping to a node which is cross-coupled to the other layer will diffuse to it with probability .
Finally we tested the sensitivity of this dynamical scenario with respect to different choices of the initial concentrations .
As shown in Fig. 2e,f, increasing values of slightly affect the dynamical behavior of the system for both ER and SF interconnected networks, mainly influencing the transients in which the system indulges before entering the first diffusive regime.
In particular, higher values of the densities make the particles in the two systems to experience earlier the repulsive effects.

At this point, it worth to notice that the phenomenology described so far partially agrees with the results presented earlier in Ref. garas () by A. Garas, where the same system was investigated from the more general perspective of different strategies of cross-systems interactions.
However, by contrast with the conclusions drown in Ref. garas (), here we have shown that in the limit of low enough fractions of interconnected nodes, an initial separation of reactants generally leads to a novel spontaneous mechanism of repulsion among the reactants which, at the best of our knowledge, has been so far overlooked.

Next, we investigate the effects that different values of the effective transmission probability has on the mixing properties of the system, so to verify the accuracy of the approximations adopted to derive the relation (1) for . To evaluate this quantity, we have plot the logarithmic derivative of the -axis in Fig. 2, and searched for the maxima of the corresponding curves (see Fig. 3). The time at which a maximum is reached defines the mixing time , sharply marking the dynamical crossover between the two regimes. As shown in Fig. 4a, we find that the time for the process on ER (and equivalently, though not shown, RR) networks depends on as predicted by Eq. (1). To further support this behavior, in the inset of Fig. 4a we have validated the logarithmic dependency of due to the repulsion mechanism by comparing the pattern observed with the constant behavior that one would have found in the case of a purely diffusive process.

By performing the same analysis on heterogeneous (SF and SR) networks, we find that the mixing time surprisingly follows the behavior predicted by the mean-field theory for homogeneous structures (Fig. 4b), though with slightly less accuracy.
We trace the origin of this result to the fact that, in our model, low-degree nodes are indeed the most important ones for reactions to occur between the two populations, as they most likely carry the cross-connections.
In fact, as percolation results would suggest SFperc1 (), the probability of picking at random a hub is very low in SF or SR networks, making even harder the chance of having hub-to-hub interconnections. In this respect, we expect that degree-based couplings among networks will likely tune the effects of the repulsion mechanism, leading to a faster mixing of the reactants for e.g. hub-to-hub couplings.

To complete the analysis of the model for this Rapid Communication, we have tested the response to the system to finite-size effects, by adopting networks with equal topologies, but different number of nodes.
In particular, for ER networks the three dynamical regimes encountered in Fig.2b remain unaltered in their main stages (Fig. 5a), and only exhibit an extinction point at earlier times for networks of decreasing sizes.
Finally in Fig. 5b we have repeated the test for SR networks, where we have found that the pattern observed in Fig. 2c is again mainly unchanged, except that the convergence of the mixing time to its thermodynamic value is monotonic in , enlightening the possible occurrence of finite-size effects in the case of power-law networks of small-size whose study will be performed elsewhere.

### III. Summary and conclusions

In this work we have studied, both numerically and analytically, the dynamics of an process on a pair of randomly coupled networks, where reactants are initially separated. For small enough values of the fraction of interconnected nodes between the layers, we have found that the inverse particle density scales linearly at short times and then crosses over to a second linear regime at time . As the crossover determines the time at which the two population start to extensively mix, we have analyzed the dependence of the mixing time on the effective diffusion rate , unveiling a novel repulsive mechanism whose spontaneous emergence delays the mixing of the reactants. We gave numerical evidence that, on randomly coupled synthetic networks, this effect does not show a sensitive dependence on the heterogeneity of the underlying topology, but it is in fact dominated by nodes with low connectivity. Whether or not the same behavior will appear on networks with targeted (e.g. hub-to-hub, or non-hub-to-hub) interconnections Aguirre () or on more realistic structure having e.g. a spatial embedding or degree-degree correlations, remains an intriguing question calling for further investigations. Moreover, since the diffusion-controlled annihilation process adopted in this work can be considered as an archetypal model for reaction kinetics Redner2 (), we believe that these results will inspire the investigation of the effects that the initial distribution of reactants, together with the mesoscopic architecture of the interconnected network, will have on the pattern formation in more elaborated and realistic spreading models Castellano2009 (); agliari (), enlarging in this way our understanding of the interplay between structure and dynamics.

## Acknowledgments

Results presented in this work have been produced using the European Grid Infrastructure (EGI) through the National Grid Infrastructures (HellasGrid) as part of the SEE Virtual Organisation. This research was supported by European Commission FP-FET project Multiplex # . S.H. acknowledges the ISF, ONR Global, the Israel-Italian collaborative project NECST, Japan Science Foundation, BSF-NSF, and DTRA (Grant # HDTRA-1-10-1-0014), together with the Israeli Ministry of Science, Technology and Space (MOST, grant #3-12072) for financial support.

### Footnotes

- thanks: F. L. and B. G. contributed equally to this paper.

### References

- M. V. Smoluchowski, Z. Phys. Chem. 92, 129 (1917).
- J. D. Murray, Mathematical Biology, Inter. App. Math. 17 (Springer, 2004).
- C. Furusawa & K. Kaneko, Science 338, 215 (2012).
- N. G. Van Kampen, Stochastic processes in chemistry and physics (Amsterdam: North Holland, 1992).
- D. Ben-Avraham & S. Havlin, Diffusion and reactions in fractals and disordered systems (Cambridge U.P., 2000).
- A. A. Ovchinnikov & Ya. B. Zeldovich, Chem. Phys. 28, 215 (1978).
- D. Toussaint & F. Wilczek, J. Chem. Phys. 78, 2642 (1983).
- F. Leyvraz & S. Redner, Phys. Rev. Lett. 66, 2168 (1991).
- M. Araujo, H. Larralde, S. Havlin & H. E. Stanley, Phys. Rev. Lett. 71, 3592 (1993).
- C. Castellano, S. Fortunato, V. Loreto, Rev. Mod. Phys. 81(2), 591 (2009).
- A. Vespignani, Nature Phys. 8(1), 32 (2012).
- L. Gálfi, Z. Rácz, Phys. Rev. A, 38(6), 3151 (1988).
- A. A. Ovchinnikov, S. F. Timashev, A. A. Belyĭ, Kinetics of Diffusion Controlled Chemical Processes (Nova Science Pub. Inc., 1989).
- H. Taitelbaum, S. Havlin, J. E. Kiefer, B. Trus, G. H. Weiss, Jour. Stat. Phys. 65(5-6), 873 (1991).
- G. Zumofen, J. Klafter & M. F. Shlesinger, Phys. Rev. Lett. 77, 2830 (1996).
- M. Newman, Networks: an introduction (Oxford U. P., 2010).
- R. Cohen & S. Havlin, Complex networks: structure, robustness and function (Cambridge U. P., 2010).
- S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, & D.-U. Hwang, Phys. Rep. 424, 175 (2006).
- A. Barrat, M. Barthelemy & A. Vespignani, Dynamical Processes on Complex Networks (Cambridge U. P., 2008).
- G. Caldarelli, Scale-free networks: complex webs in nature and technology (Oxford U. P., 2007).
- S. N. Dorogovtsev, A. V. Goltsev & J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
- R. Pastor-Satorras, C. Castellano, P. V. Mieghem & A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
- L. K. Gallos & P. Argyrakis, Phys. Rev. Lett. 92, 138301 (2004).
- L. K. Gallos & P. Argyrakis, Phys. Rev. E 72, 017101 (2005).
- M. Catanzaro, M. Boguñá, & R. Pastor-Satorras, Phys. Rev. E 71, 056104 (2004).
- S. Weber & M. Porto, Phys. Rev. E 74, 046108 (2006).
- L. K. Gallos & P. Argyrakis, J. Phys.: Cond. Matt. 19, 065123 (2007).
- M. Catanzaro, M. Boguñá, & R. Pastor-Satorras, Reaction-diffusion processes in scale-free networks in Handbook of Large-Scale Random Networks, 203 (2008).
- S. Gómez, A. Díaz-Guilera, J. Gómez-Gardeñes, C. J. Pérez-Vicente, Y. Moreno & A. Arenas, Phys. Rev. Lett. 110, 028701 (2013).
- A. Garas, Phys. Rev. E 92, 020801(R) (2015).
- S. Havlin et al., Eur. Phys. J., Special Topics 214, 273 (2012).
- S. Boccaletti, G. Bianconi, R. Criado, C. Del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendiña-Nadal, Z. Wang, & M. Zanin, Phys. Rep. 544, 1 (2014).
- M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno & M. A. Porter, J. Compl. Net. 2, 203 (2014).
- G. Bianconi, Eur. Phys. Lett. 111, 56001 (2015).
- A. Garas, Ed. Interconnected networks (Springer, 2016).
- S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley & S. Havlin, Nature 464, 1025 (2010).
- J. Gao, S. V. Buldyrev, H. E. Stanley % S. Havlin, Nat. Phys. 8, 40 (2012).
- D. F. Klosic, A. Grimbs, S. Bornholdt & M.-T. Hütt, Nat. Comm. 8:534 (2017).
- C. del Genio, J. Gómez-Gardeñes, I. Bonamassa & S. Boccaletti, Sci. Adv. 2, e1601679 (2016).
- V. Nicosia, P. S. Skardal, A. Arenas & V.Latora, Phys. Rev. Lett. 118, 138302 (2017).
- M. M. Danziger, I. Bonamassa, S. Boccaletti & S.Havlin, arXiv:1705.00241 (2017).
- B. Bollobás, Random graphs, in Modern graph theory (Springer, New York, NY, 1998).
- M. Catanzaro, M. Boguñá & R. Pastor-Satorras, Phys. Rev. E 71, 027103 (2005).
- After the crossover, the reactants become (so to say) oblivious of the “mesoscopic” structure of the interconnected network. This means that the fraction of interconnected nodes must be much larger then the density of the survived particles, which justifies the condition adopted in deriving Eq. (6).
- V. Colizza, R. Pastor-Satorras & A. Vespignani, Nature Physics 3, 276 (2007).
- A. Baronchelli, M. Catanzaro & R. Pastor-Satorras, Phys. Rev. E 78(1), 016111 (2008).
- R. Cohen, K. Erez, D. ben-Avraham & S. Havlin, Phys. Rev. Lett. 85, 4626 (2000).
- J. Aguirre, D. Papo, J. M. Buldú, Nature Physics, 9(4), 230 (2013).
- P. L. Krapivsky, S. Redner, E. Ben-Naim, A kinetic view of statistical physics (Cambridge Uni. Press, 2010).
- E. Agliari, R. Burioni, D. Cassi, F. Maria Neri, IMA Jour. Man. Math., 21(1), 67 (2009).