We study the transient dynamics of an A+B→0 process on a pair of randomly coupled networks, where reactants are initially separated. We find that, for sufficiently small fractions q of cross couplings, the concentration of A (or B) particles decays linearly in a first stage and crosses over to a second linear decrease at a mixing time tx. By numerical and analytical arguments, we show that for symmetric and homogeneous structures tx∝(〈k〉/q)log(〈k〉/q) where 〈k〉 is the mean degree of both networks. Being this behavior is in marked contrast with a purely diffusive process, where the mixing time would go simply like 〈k〉/q, we identify the logarithmic slowing down in tx to be the result of a spontaneous mechanism of repulsion between the reactants A and B 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.
Nearly 100 years ago, Marian von Smoluchowski introduced a mathematical model to describe coagulation phenomena in terms of diffusion-controlled reaction processes . Despite its apparent simplicity, the kinetics of this model was found to yield a wealth of intriguing phenomena, whose analyses have widely enriched our understanding of pattern formation in chemical compounds, biological systems [2,3], and elsewhere.
From the statistical mechanics perspective, reaction-diffusion (RD) processes represent a fertile groundwork in which to analyze the emergence of spontaneous mechanisms by starting from microscopic rules . 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 .
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 [6–9]. These phenomena result in a drastic slowing down in the rate of the reactions, forcing the system in a long-lived nonequilibrium state with a sublinear 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 [10,11], or underlying the pattern formation of diverse chemical reaction [12–14], the appearance of subdiffusive 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 superdiffusive dynamics .
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 [16,17]. 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 .
Consistently with the scenario observed in other models [19–21], complex networks significantly influence the collective properties of RD processes . Numerical [23,24] and theoretical  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 superdiffusive behaviors , 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 [27,28].
Although the portrait of RD processes on isolated structures is nowadays clear, not much is known regarding their dynamical behaviors on multilayer networks [29,30]. The influence that their mesoscopic organization has on the collective behaviors of processes acting on them is attracting significant interest [31–35], and has already produced interesting results [36–38]. 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 [39–41]. Following this mainstream, we study here the dynamics on a pair of randomly coupled networks with initially separated reactants, where we find a 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
This Rapid Communication is organized as follows. We derive analytically Eq. (1) for RR graphs, which 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.
Analytic approach. We consider two configuration model networks , composed by the same number of nodes and same structural properties, a condition that we will refer to 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 population, let us denote by and the concentration of particles in networks 1 and 2, respectively; similarly, let and be the concentration of particles in networks 1 and 2, respectively. Having assumed equal initial conditions, we have and . Moreover, by symmetry, and for all . It is worth noting here that, while 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
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 quasistationary behavior, namely,
Equation (5) can be adopted in order to define the mixing time , which 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., . In this approximation, we get
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 cutoff given by , where is the minimum degree of each node .
Two populations 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 [45,46]. 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 300, 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. 2(a)], ER [Fig. 2(b)], and SF networks [Fig. 2(c)]. 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. 2(b), 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. 2(b) 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. 2(d), 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 Figs. 2(e) and 2(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 experience the repulsive effects earlier.
At this point, it is worth noting that the phenomenology described so far partially agrees with the results presented earlier in Ref.  by Garas, where the same system was investigated from the more general perspective of different strategies of cross-system interactions. However, by contrast with the conclusions drown in Ref. , here we have shown that in the limit of low enough fractions of interconnected nodes, an initial separation of reactants generally leads to a spontaneous mechanism of repulsion among the reactants which, to the best of our knowledge, has been so far overlooked. Next, we investigate the effects that different values of the effective transmission probability have on the mixing properties of the system, so as to verify the accuracy of the approximations adopted to derive the relation (1) for . To evaluate this quantity, we have plotted 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. 4(a), we find that the time for the process on ER (and equivalently, although not shown, RR) networks depends on as predicted by Eq. (1). To further support this behavior, in the inset of Fig. 4(a) 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. 4(b)], although 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 , the probability of picking a hub at random is very low in SF or SR networks, making the chance of having hub-to-hub interconnections even harder. 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. 2(b) remain unaltered in their main stages [Fig. 5(a)], and only exhibit an extinction point at earlier times for networks of decreasing sizes. Finally in Fig. 5(b) we have repeated the test for SR networks, where we have found that the pattern observed in Fig. 2(c) 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.
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 populations start to extensively mix, we have analyzed the dependence of the mixing time on the effective diffusion rate , unveiling a 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  or on more realistic structures 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 , 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 [10,50], 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 Organization. This research was supported by European Commission FP7-FET Project Multiplex No. 317532. S.H. acknowledges the ISF, ONR Global, the Israel-Italian collaborative project NECST, Japan Science Foundation, BSF-NSF, and DTRA (Grant No. HDTRA-1-10-1-0014), together with the Israeli Ministry of Science, Technology and Space (MOST, Grant No. 3-12072) for financial support.
F.L. and B.G. contributed equally to this work.