Volume explored by a branching random walk on general graphs

Branching processes are used to model diverse social and physical scenarios, from extinction of family names to nuclear fission. However, for a better description of natural phenomena, such as viral epidemics in cellular tissues, animal populations and social networks, a spatial embedding—the branching random walk (BRW)—is required. Despite its wide range of applications, the properties of the volume explored by the BRW so far remained elusive, with exact results limited to one dimension. Here we present analytical results, supported by numerical simulations, on the scaling of the volume explored by a BRW in the critical regime, the onset of epidemics, in general environments. Our results characterise the spreading dynamics on regular lattices and general graphs, such as fractals, random trees and scale-free networks, revealing the direct relation between the graphs’ dimensionality and the rate of propagation of the viral process. Furthermore, we use the BRW to determine the spectral properties of real social and metabolic networks, where we observe that a lack of information of the network structure can lead to differences in the observed behaviour of the spreading process. Our results provide observables of broad interest for the characterisation of real world lattices, tissues, and networks.

www.nature.com/scientificreports www.nature.com/scientificreports/ To extract the number of distinct sites visited, we introduce an immobile tracer particle species with occupation numbers m x . They are spawned as offspring by the active walkers with rate γ at the sites they are visiting, thereby leaving a trail of tracers behind, similar to the breadcrumbs left by Hänsel and Gretel 12 , Fig. 1a. We impose the constraint that at most a single tracer can reside at any given site, which means that the spawning of a tracer is suppressed in the presence of another tracer. It is that suppression that generates significant complications from the point of view of the stochastic process. Yet, only with this restriction in place is the number of tracers a measure of the number of distinct sites visited by the walkers, as pictured by the cloud of visited sites in Fig. 1b.
There is no interaction between active and tracer particles, other than at the spawning of immobile tracers by active walkers. In principle, the spawning (attempt) rate γ has to diverge in order to mark every single site visited by the walkers. However, it turns out that this limit is irrelevant as far as the asymptotic features of this process at large system sizes and long times are concerned 5 .
By definition, the sets {n} and {m} of occupation numbers n x and m x , respectively, for each site x of a given graph, are Markovian and a master equation can be written for the joint probability  n m t ({ }, { }; ) to find the graph in a certain configuration of occupation numbers at time t , and the terms on the right-hand side, ), indicate the contributions from branching s, extinction of active walkers e and tracer particles ε′, hopping H and deposition γ, respectively (see Sec. S1 for details). We constructed a statistical field theory from the master Eq. (1) using the ladder operators introduced by Doi 13 and Peliti 14 (methods Sec. V A). To regularise the propagators of the immobile particles in the field theory, we allow for the extinction of immobile particles with rate ε′ in Eq. (1), not dissimilar to the birds that foiled Hänsel and Gretel's plans (Fig. 1a). The propagators for active and tracer particles do not renormalise, and the limit ε′ → 0 is taken before any observable is evaluated. Through field-theoretic renormalisation in dimensions d = 4 − ε we can then determine the exact scaling behaviour of the number of distinct sites visited by the walkers. The branching process described by Eq. (1) has three regimes, as becomes evident in the field-theoretic formulation, where a net extinction rate r = e − s appears. This net extinction rate is not renormalised in the field-theory and therefore no mass shift appears. The BRW is subcritical for r > 0, critical for r = 0 (onset of epidemics) and supercritical for r < 0. Hereafter, we focus on the critical case, where fluctuations dominate the dynamics, and the behaviour becomes unpredictable and highly volatile. Furthermore, for both analytical and numerical computations we consider the initial condition of a single walker at t = 0. Extensions to different initial conditions are straight-forward.

Results for Regular Lattices
Following the field theoretic approach (details in Secs. V A and Sec. S2) of the bulk critical behaviour in the continuum limit, where hopping is replaced by diffusion by introducing a diffusion constant D, we find that in the thermodynamic limit at long times t, the expected number of distinct sites visited or the volume explored, 〈a〉(t, L), scales like t (d−2)/2 in dimensions d < 4. In dimensions d < 2 this volume remains finite in large t. The scaling of the p-th moment of the number of distinct sites visited follows,  Tables S1 and S2). The process is free beyond d c = 4 dimensions, where the probability of any walker or any of its ancestors or descendants ever to return to a previously visited site drops below unity, and the scaling becomes independent of the dimension, The gap-exponent in dimensions greater than d c = 4 is thus 4, as confirmed by numerical observations in dimension d = 5 (see Fig. 2d and Table S2). As correlations become irrelevant, this is usually referred to as mean-field behaviour. The set of sites visited may thus be regarded as a four-dimensional object, projected into the d-dimensional lattice considered. Focusing on dimensions below d c = 4, the distribution of the number of distinct sites visited, a, follows a power law, with metric factor A and cutoff a c~( Dt) d/2 for Dt ≪ L 2 and a c ~ L d otherwise. These results show how increasing the dimensionality of the lattice promotes the appearance of larger events, evidencing the relevance of dimension on the spreading. In dimensions d ≥ d c = 4 the resulting scaling of the distribution is that of Eq. (4) at d = 4, where the probability distribution decays like a −3/2 . Numerically, we recorded, for each realisation, the total number of distinct sites visited by the process in order to construct the distribution, a ( )  , of sites visited. The numerical results coincide with our theoretical predictions, as shown in Fig. 3a.
The exponents found above for d = 1 are in agreement with the exact solution by Ramola et al. 8 , where  a ( ) decays as a −3 . In two dimensions the power-law tail decays as a −2 , which coincides with the decay of the 2d convex hull area distribution 4 .

extension to General Graphs
In the field theoretic approach followed to find the scaling in Sec. II the spatial dimension of the lattice enters only in as far as its spectral dimension is concerned, which characterises the density of eigenvalues of the Laplace operator on the graph given. Our results extend naturally to all translational invariant lattices and graphs, by replacing the dimension d of the lattice in Eqs (2)-(4) by the spectral dimension d s of the graph, as detailed in Sec. S5. This holds true more generally as long as the lattice Laplacian itself does not undergo renormalisation, i.e. in the absence of an anomalous dimension 16 . In the study of networks the number of nodes N, is a more natural measure of the size of the network than the linear size L. Using L ~ N 1/ d s we can write the scaling of the BRW in time and number of nodes as Here, the gap-exponent for the scaling in number of nodes is always unity. This extension to graphs allowed us to predict the behavior of the BRW spreading in both artificial networks relevant for social and biological sciences, and complex systems in general 2,17-19 , as well as real networks. To illustrate this, we considered first the Sierpinski carpet (SCs) (Fig. 4a, methods Sec. V B 2), and random trees (RTs) (Fig. 4b and methods Sec. V B 3). Both of these graphs are widely applied in the context of porous media 20 and percolation 21 Fig. 4d,e. These theoretical predictions extend also to the distribution of visited sites (see Fig. 3b), by setting d = d s in (4). www.nature.com/scientificreports www.nature.com/scientificreports/ Furthermore, we studied the BRW behaviour on a class of scale free networks 24 . Since their introduction, scale free graphs have been observed to describe a plethora of natural phenomena, including the World-Wide-Web 25 , transportation 26 , and metabolic networks 27 , to name but a few. We considered a preferential attachment scheme 24 (Fig. 4c, see methods Sec. V B 4), to construct networks with power-law degree distribution (Fig. S1). The existence of a finite spectral gap in these networks, which indicates slow decay of return times 28,29 , suggested that the BRW process is bound to exhibit mean-field behaviour, i.e. d s ≥ 4. This was confirmed by numerical simulations, where the probability distribution of visited sites (Fig. 3b) has a power-law decay with exponent −1.52(2) ≈ −3/2, and the scaling in time and system size ( Fig. 4f and Table S1) follow mean-field behaviour as predicted by (3) for The spectral dimension gives information on the behaviour of dynamical processes on graphs. Here we use the BRW to characterise real-world networks through the power-law decay of the distribution of visited sites  a ( ), which according to Eq. (4) is − + a d (1 2/ ) s provided d s ≤ 4. For example, the BRW exhibits near mean field-behaviour on a subset of the Facebook network, which has been characterised as scale-free 30 . Hence, we derived a large effective (spectral) dimension, d s = 3.9(1), indicating a fast spreading of the viral process in this network (see Fig. 5).
We should emphasize that the spectral dimension is sensitive to changes in network topology and connectivity. To exemplify this we have considered two publicly available datasets for the yeast protein interaction network (see Fig. 5). We found that even though both network describe subsets of the same biochemical network, namely the complete yeast protein interactome, the spectral dimensions in both cases are significantly different, d s = 3.0(1) for the network with N = 1870 nodes 31 , and d s = 3.8(1) for the larger network of N = 2559 nodes 32 , leading to differences in properties of the spreading process among the two. The discrepancy points to differences in the connectivities of both networks and shows the importance of having access to the complete network in order provide a reliable analysis of its properties, which may have biological implications 33,34 .

outlook
The results presented above for the binary branching process, where walkers branch into exactly two new walkers, apply equally to more general branching processes, where the number of offspring in each birth event is given by a distribution (for details see Sec. S4). This can be seen, for example, in real-world scenarios where a single infected individual or device infects a whole neighbourhood around them, or in the case of signal propagation in  www.nature.com/scientificreports www.nature.com/scientificreports/ protein networks, where the activation of one node (or chemical reaction) can in turn activate a whole fraction of its neighboring nodes.
While the scaling behaviour does not depend on the initial position x 0 of a walker, provided it is located in the bulk and remains there as the thermodynamic limit is taken, the field theory has to be adjusted to account for more complicated boundary conditions 5 or the walker starting close to any such boundary. It may also be interesting to consider the case of initialising each site with an independent Poisson distributions of walkers 35 .
The approach followed in the present work provides a quantitative measure to explore and determine the spectral dimension of artificial and real networks. This is of particular interest when the spectral dimension is greater or equal to 2, where the traditional approach of exploring graphs, based on simple random walks 29,36 , fails. When simulating the BRW we made the observation that robust scaling is more easily obtained on small lattices if the  www.nature.com/scientificreports www.nature.com/scientificreports/ hopping rate H is clearly smaller than the rates of branching s and extinction e. For large values of the hopping rate particles leave the system during the initial transient, as seen in Fig. 2, thus boundary effects appear before any robust scaling can be observed. In graphs such as the PA network (Fig. 4), that does not have any boundaries, these artefacts are much less pronounced. In summary our results shed new light on the properties of spatial branching processes on general graphs, and their applicability in the study of real complex networks, and provide observables of broad interest for the characterisation of real world lattices, tissues, and networks.

Methods
Field theory of the BRW. In order to derive the main results for the scaling of distinct sites visited by the BRW (Sec. II) we work along established lines 37 , casting the master equation in a field theory of the annihilation fields φ(x, t) and ψ t x ( , ) for the active and the immobile particles, respectively, and of the corresponding (Doishifted) creation fields φ  t x ( , ) and ψ ∼ t x ( , ). The governing Liouvillian = + 0 1    consists of a harmonic part, and a non-linear part,  φ ψ φ ψ φ φ σψ φ φ λψ ψφ ξψ ψφ φ κψ ψφ φ χψ ψφ where we have taken the continuum limit. The space and time integrated Liouvillian produces the field-theoretic