From Logistic Growth to Exponential Growth in a Population Dynamical Model
FFrom Logistic Growth to Exponential Growth in aPopulation Dynamical Model
Inbar Seroussi and Nir Sochen Department of Mathematics, Weizmann Institute of Science, POB 26, Rehovot76100, Israel Department of Applied Mathematics, School of Mathematical Sciences, Tel-AvivUniversity, Tel-Aviv, 69978, IsraelAugust 2, 2019
Abstract
Dynamics among central sources (hubs) providing a resource and large number of com-ponents enjoying and contributing to this resource describes many real life situations. Mod-eling, controlling, and balancing this dynamics is a general problem that arises in manyscientific disciplines. We analyze a stochastic dynamical system exhibiting this dynamicswith a multiplicative noise. We show that this model can be solved exactly by passing tovariables that describe the mass ratio between the components and the hub. We derive adeterministic equation for the average mass ratio. This equation describes logistic growth.We derive the full phase diagram of the model and identify three regimes by calculatingthe sample and moment Lyapunov exponent of the system. The first regime describes fullbalance between the non-hub components and the hub, in the second regime the entire re-source is concentrated mainly in the hub, and in the third regime the resource is localizedon a few non-hub components and the hub. Surprisingly, in the limit of large number ofcomponents the transition values do not depend on the amount of resource given by thehub. This model has interesting application in the context of analysis of porous media usingMagnetic Resonance (MR) techniques.
Population dynamics on large scale networks has attracted a lot of attention due to its wideoccurrence in many disciplines, such as social sciences [1, 2], physics [3] and biology, communica-tion and control theory [4]. This dynamics is mainly affected by the topology of the network aswell as some internal stochastic noise. In many applications there are only a few nodes playing amajor role in the dynamical process, distributing and carrying most of the resources [5, 6]. Forexample, this can be the case in models describing population dynamics, economic growth [7],and distributed control systems [4]. An additional application of this problem is in the contextof diffusion measurements of porous systems, such as brain tissue, using Magnetic Resonance(MR) techniques. In this last case, the sensitivity of the MR signal to self-diffusion of watermolecules can be utilized to extract information about the network of cells (neurons) in thebrain. The concept of self-diffusion of molecules in a network of pores was already introduced in1 a r X i v : . [ phy s i c s . d a t a - a n ] J u l efs. [8, 9, 10]. The main challenge is how to determine the topology of the network based onthe MR measurements [8].Our model consists of a system of interacting sites on a graph G with N vertices and E edgesbetween them. We are interested in the stochastic dynamics of some characteristic property { m i ( t ) } i ∈G ,t ≥ . The property m i ( t ) is linked to a physical measurable quantity in the real worldand the graph is the underlying geometry/topology in which the property lives. The topology is acomplex network of sites. The model is described by the following family of stochastic differentialequations in the Stratonovich form on the graph G : dm i ( t ) dt = (cid:88) j J ij m j ( t ) − (cid:88) j J ij m i ( t ) + g i ( t ) m i ( t ) , (1)with the initial conditions m i (0) = m . The term g i ( t ) is a multiplicative white noise, suchthat (cid:104) g i ( t ) (cid:105) = f i , and (cid:104) g i ( t ) g j ( t (cid:48) ) (cid:105) = σ i δ ij δ ( t − t (cid:48) ). We choose the Stratonovich form, sinceits solution is a limiting case of a physical system involving white noise with short memory[11]. The topology of the network is encoded in the adjacency matrix J of the graph. Themodel consists of two parts: an interacting part, where the interaction strength depends on thelocation on the graph, and a non-interacting part, where each component follows a stochasticnoise with different variance σ i . The first part causes spreading, while the second pushes towardsconcentration (a.k.a localization or condensation). The model was already analyzed in the meanfield topology, i.e., when all the nodes are connected and interact at the same rate. In this case,the equilibrium distribution is a Pareto power-law [1]. It was also analyzed on trees [12, 13],and random graphs assuming separable probability distribution on the nodes [14]. The model onthe lattice is known in the mathematical literature as the time-dependent Parabolic AndersonModel (PAM) [15, 16]. The phase diagram of the model in this case depends on the dimensionof the lattice. On a general network, phase transitions depend on the spectral dimension of thenetwork [8].Here, we present and analyze a specific topology in which the model is shown to be solvable.Namely, we consider a directed graph with N + 1 nodes, one hub node interacting with N independent nodes. In the context of MR measurements of diffusion in a porous structure, theMR signal measured is assumed to be composed of two contributions: one coming from hindereddiffusion in the extracellular space and the other from restricted diffusion in the intracellularspace [17, 18, 19, 20]. The hub node represents the magnetization in the extracellular space(e.g., water), h ( t ), and the non-hub nodes represents N independent intracellular pores withmagnetization, m i ( t ). The motion of molecules between these regions changes the value of themagnetization as a function of time and is represented by the interaction term between nodes, J ij .The effect of the magnetic field gradient can be incorporated in the stochastic noise, for example,in f i , and/or its variance σ i . In the economic context, the system describes the dynamics ofthe money hold by the hub, which represents by the state/bank, and the money of each agent m i ( t ). In this case, the agents deposit money in the bank and the bank pays interest on it. Thestochastic noise represents the bank/state and the agent’s investments in the stock market andhousing [1, 7]. Analysis of the dynamics of the sums of the money held by the agents and thebank/state was curried out in Ref. [7].Here, we analyze the dynamics of the mass ratio between each agent and the bank/state(hub). We derive the equilibrium distribution in this case, and show that when the number ofnodes growth at least exponentially with time there exist a localization phase. We also generalizeour analysis to multiple number of hubs.Our main result is a full phase-diagram of the model. We show that this model can bedescribed by a stochastic equation for the mass ratio between each of the non-hub nodes and2 (𝑡) 𝑚 𝑖 (𝑡) 𝐽 out 𝑁𝐽 in 𝑁 Figure 1: An illustration of the system in the hub topology.the hub, and a deterministic non-linear equation for the average relative mass of all the non-hubnodes with respect to the hub. To identify the phases of the system, we calculate the sampleand moment Lyapunov exponents and identify a gap between them. The phase transitions arecharacterized by one parameter. This parameter takes into account the exchange rate betweenthe non-hub nodes and the hub and the variances of the multiplicative noises.
The basic hub topology is composed of an infinite number of nodes, { m i } , interacting at constantrate with a hub node, h , such that, J i = J out N , and J j = J N , respectively. Our normalization issuch that, the overall interaction between the nodes and the hub is finite in the limit of infinitenumber of nodes. The interaction among the non-hub nodes is characterized by the parameter δ ;when δ = 0, any interaction (transfer of mass) between the non-hub nodes is done only throughthe hub. The topology of the interaction between the non-hub nodes is defined by a Laplacianmatrix, L , satisfying (cid:80) i L ij = 0. Figure 1 presents an illustration of such a system for δ = 0.We assume that the stochastic noise acting on the non-hub nodes has the same variance andaverage for all nodes, σ i = σ out , f i = f . Eq. (1) in the Itˆo form reduces to the following systemof stochastic equations: dh dt = J in N (cid:88) j m j − J out h + f h + σ h + σ g h , (2) dm i dt = J out N h − J in N m i + f m i + σ m i − δ (cid:88) j L ij m j + σ out g i m i . (3)The Itˆo form will be useful latter on when one takes the limit N → ∞ . It is instructive to passto the following normalized variables by introducing the mass ratio between the non-hub nodesand the hub node, M i = m i h . This leads to the following system of N equations in the Itˆo form3See A for more details): dM i dt = J out N − (cid:18) J in N + ∆ f − J out − σ J in M ( t ) (cid:19) M i − δ (cid:88) j L ij M j + σξ i M i . (4)We introduce the average mass ratio M = N (cid:80) i M i , the effective variance σ = σ + σ , andthe average difference ∆ f = f − f . Here, ξ i ( t ) is a Gaussian process with mean zero andvariance one. The transformation of variables introduces a non-linear term, which accounts forthe interaction between any two nodes through the hub. In the limit N → ∞ , averaging over allthe nodes in Eq. (4) yields the following deterministic law for the average mass ratio: dMdt = σ β M (cid:18) α + 1 β − M (cid:19) , (5)where we use the dimensionless parameters α = 2 J out − ∆ fσ and β = J in σ . Since, Eq. (5) it is anon-linear equation, there are two steady-state solutions: M → J out − ∆ f + σ J in = α +1 β as N → ∞ ,and M → N → ∞ . Note that, in the case of δ = 0 the system geometry can be viewed asa directed tree structure with one level and infinitely many leaves. Therefore, the presence of athe non-linear term is caused by to the indirect interaction between the non-hub nodes followingthe tree topology of the system [12]. Convergence to each one of these fixed points depend on theinitial condition of the system and on the stability of these points. Stability analysis of these twopoints shows that the point M is stable when α + 1 >
0, while the point M = 0, is stablewhen α + 1 ≤
0. For example, in the context of MR measurements in porous media, this systemcan model a complex structure measured from a single voxel in the MR image. The value M eq describes the steady-state average magnetization ratio, between the extracellular space and theintracellular space. The first fixed point M is reached when the steady-state magnetizationratio is equal to the amount of molecules leaving the non-hub pores reduced by the magneticfield effects divided by the amount of molecules leaving the extracellular space. The second fixedpoint M represents the case in which on average most of the contribution to the magnetizationin a single voxel comes mainly from the extracellular space. Eq. (5) shows logistic growth and isa version of Lotka-Volterra equation, which describes many social phenomena in nature [21, 22].It can be solved exactly: for an initial condition M (0) = M , the solution is M ( t ) = M ( α + 1) β (cid:16) M + (cid:16) α +1 β − M (cid:17) e − σ α +1) t (cid:17) . (6) Given Eqs. (4) and (5) for the relative mass between the non-hub nodes and the hub, onecan derive an equivalent form describing the dynamics of the probability distributions of M i .Since the dynamics of the average mass ratio between the non-hub nodes and the hub node isgoverned by a deterministic non-linear equation, in the limit N → ∞ the system reduces to aset of stochastic independent equations for the relative mass in each node, M i . Therefore, wecan omit the index i , and look at the dynamics of the probability distribution of a typical node P ( M, t ). This dynamics of the probability distribution is described by the following Fokker-Plankequation: ∂P∂t = − σ ∂ (cid:16)(cid:16)(cid:16) α + 1 − β M − ˜ δ (cid:17) M + ˜ δ M ( t ) − ˜ δ (cid:17) P (cid:17) ∂M + σ ∂ (cid:0) M P (cid:1) ∂ M . (7)4o study the effect of interaction between nodes, we introduce a small mean field interactionbetween the non-hub nodes, δ , and define the dimensionless interaction rate, ˜ δ = δσ . By equatingthe left hand side to of Eq. (7) zero one can find the steady-state distribution. The solutionshows a Pareto power-law behavior: P eq ( M ) = A exp (cid:32) − ˜ δ M eq M (cid:33) M − µ ( M eq ) . (8)The power is a function of the steady-state average relative mass M eq , i.e., steady-state solutionof Eq. (5): µ ( M eq ) = βM eq + 1 − α + 2˜ δ . Substituting the value of the average mass ratio, wefind that µ ( α +1 β ) | α +1 ≥ = 2 + 2˜ δ , and µ (0) | α +1 < → − α + 2˜ δ . This shows that the systemhas two steady-states. The system collapses to one of them depending on the initial condition,i.e., the initial mass ratio between the hub and the non-hub nodes. Note that the value of µ isgreater than 2 when the system collapses to the state M eq = α +1 β , showing equality among thenon-hub nodes and the hub. For, M eq = 0, the mass is localized on a few nodes within the setof non-hub nodes. The power is, µ <
2, as long as 2˜ δ < α . Therefore, adding interactionbetween non-hub nodes reduces localization, as expected. The analysis above reveals the localization regime within the non-hub nodes when the influenceof the hub is renormalized. In this section, we analyze the regime at which there is localizationon the hub. For this purpose, we study the asymptotic properties of the total mass, E ( t ) = h ( t ) + (cid:80) i m i ( t ). We calculate the Lyapunov exponents of the solution [15, 16]. Here, weperform the analysis for the case of δ = 0 . The Lyapunov exponents describe the growth rates ofthe solution and its moments. They indicate the level of complexity in the solution’s landscape.The first moment Lyapunov exponent of the solution is as follows: γ = lim t,N →∞ t ln ( (cid:104) E ( t ) (cid:105) ) = (cid:40) f + σ f − J out + σ σ σ < α ∆ σ σ ≥ α , (9)where ∆ σ = σ − σ is the variance difference, see B for details of the proof together withcorrections for finite network size. Interestingly, what determines the growth on average is thedifference between the variance, ∆ σ of the stochastic noises of the hub and that of the non-hubnodes. To understand Eq. (9), let us look at the limit where the stochastic noise in the non-hubnodes has significantly higher variance compared to the hub. This is equivalent to the presence oflarge fluctuations in the non-hub pores with respect to the extracellular space, i.e., ∆ σ ≈ − σ .In this case, comparing to the stability points in Eq. (5) in the regime α ≤ −
1, there is a highconcentration of magnetization on the hub and a few non-hub nodes which contribute most tothe total growth rate. For α > −
1, in the detailed balance limit the non-hub nodes and thehub contribute equally to the total growth. On the other hand, in the limit of ∆ σ ≈ σ , thesystem exhibits three regimes: for α > − < α ≤ α < −
1, thereis concentration of the magnetization on the hub and/or several non-hub nodes. This analysisis verified by calculating the value of the sample Lyapunov exponent, which provide the samplegrowth rate.We define the sample Lyapunov exponent, (cid:101) γ , as the limit of the logarithm of the total massof the solution divided by time as N, t → ∞ . Knowing the dynamics of the average mass ratio,5ee Eq. (5), we can calculate the sample growth rate of the mass ratio exactly. The resultingsample/quenched Lyapunov exponent is (cid:101) γ = lim t,N →∞ t ln ( h + m ) = f + σ , (10)where we denote m = (cid:80) i m i ; see C for details of the proof. Note that, in order to take the limitone needs to specify at which rate the number of nodes growth with time. We show that when thenumber of nodes growth at least exponentially with time, then the sample Lypunov is as in Eq.(10). This value is independent of the initial conditions and is bounded from above by the firstmoment Lyapunov exponent, γ . Localization of the solution is defined as the regime at whichstrict inequality hold, (cid:101) γ < γ [15]. The gap ∆ γ = γ − (cid:101) γ between these two exponents i.e., thedifference between the expression in Eq. (9) and Eq. (10), characterizes the localization regime inthe system. Combining the transition in the values of the exponent (cid:101) γ with the stability analysisof the steady-state solutions of Eq. (5), we identify three regimes: a regime of full equality, for α > ∆ σ σ , in which the mass is spread equally between all the non-hub nodes and the hub, asecond regime for − < α ≤ ∆ σ σ , in which the mass is localized mainly on the hub, and a thirdregime for α ≤ −
1, in which the mass is localized on the hub and a few non-hub nodes.
In this section, we consider the effect of H hub nodes, h i , connected to all the nodes in thesystem, and a set of N non-hub independent nodes, m i , connected only to the hubs nodes.Figure 2 illustrate of this topology. This kind of topology appears in many applications, forinstance, in the economic setting in the presence of more than one central bank/company. Ina porous structure, it can describe different extra-cellular regions interacting with intracellularpores. It is also applied in analyzing the dynamics of control systems [4], and in machine learningalgorithms. The equations of the system now read as follows: dh i dt = J in HN (cid:88) i m i − J out H h i + f i h i + δH − H (cid:88) j (cid:54) = i h j − δh i + σ i g hi h i , (11) dm i dt = J out N H (cid:88) i h i − J in N m i + q i m i + ν i g mi m i . (12)Note that, as before, the total interaction rates between the hubs and the non-hub nodes, J in , and J out , are defined to be finite in the limit of infinite N and H . The processes g hi ( t ) and g mi ( t ),are Gaussian processes with mean zero and variance one. Similar to the one-hub topology, we cannow pass to the relative mass parameters by dividing by the average hubs mass, see D for moredetails. The results of the previous sections are recovered when H = 1. Note that the equationsfor the relative mass are decoupled in the case H = 1 and also in the limit of H very large. Inthe presence of a finite small number of hub nodes, one can show that the total variance dependson the hubs value. Therefore, having finite number of hubs decreases the value of the parameter α , and causes more equality in the system and reduces the localization. Moreover, in the simplecase where all the hubs have the same statistics, such that the stochastic noise, g hi ( t ) = g ( t ),does not depend on the hub location, H hubs are equivalent to one hub with the total effectivenet flux J out . The limit of one non-hub node with J out = J in = δ , N = 1, and an infinitenumber of hubs, H → ∞ , is the mean field model with exponential growth of the total mass [1].6 in 𝑁𝐻 𝐽 out 𝑁𝐻 𝑚 𝑖 (𝑡) ℎ (𝑡) ℎ (𝑡) ℎ (𝑡) 𝛿 Figure 2: An illustration of the multi-hub topology for M = 3, and, N = 7.Note that when both the number of non-hub nodes and the number of hubs are very large, i.e., N → ∞ and H → ∞ , the average mass ratio grows exponentially. The exponent depends onthe average and variance difference of the stochastic multiplicative noises. The phase transitionin this case is equivalent to the results in Sec. 2, i.e., there are three phases: localization on thenon-hub nodes, localization on the hubs and equal spreading over all the nodes. Note that in thislimit, the non-hub nodes play the same role as the hubs, since they are connected to infinitelymany hub nodes. The interaction among the hubs, δ reduces localization, but doesn’t effectthe growth rate. The analysis above affects mainly the sample Lyapunov exponent. The firstmoment Lyapunov exponent remains the same for any H and N , since the average equationsdo not change. This shows that the phase transitions predicted in Sec. 2.2 are general and canbe observed with small modifications to a system of multiple hubs. In addition, the transitionbetween a logistic growth in the relative mass to an exponential growth is a function of the ratiobetween the N and H . In the context of MR measurements the model in Eq. (1) is a generalization of the K¨argermodel [17], which accounts for random changes in the diffusivity due to restricted diffusion or anon-homogeneous magnetic field. This model was already analyzed on a general network, wherethe importance of the spectral dimension as a measurable parameter is stated [8]. The hubtopology (see Figure 1) is a simplified version of this model. We show that in this topologyunder the assumption that all the non-hub pores have similar properties, the average equationsof the model are those for the K¨arger model for two compartments [17], see B. The parametersare then as follows: f = − q D ex and f = − q D in , σ = − q σ , σ = − q σ , such that theparameter α ( q ) = 2 J out − ∆ fσ = − J out + q ( D in − D ex ) q ( σ + σ ) , is controlled by the gradient of the appliedmagnetic field, incorporated in the value of q . For example consider the basic Stejskal-Tannersequence [23], which is composed of two gradients pulses of the magnetic field with magnitude G in opposite direction and with duration δ . The pulses are separated by a diffusion time ∆. For q one takes the wave vector, q = δγG π , where the parameter, γ , is the gyro-magnetic ratio. Theaverage equations for the magnetization in the ( x, y ) − plane, and with a magnetic field gradientin the ˆ z direction, read d (cid:104) h (cid:105) dt = J in N (cid:104) m (cid:105) − J out (cid:104) h (cid:105) − q D ex (cid:104) h (cid:105) , (13)7a) N -2.5-2-1.5-1-0.50 -+ (b) N -2.5-2-1.5-1-0.50 q=0[1/cm] q=50[1/cm]q=100[1/cm]q=150[1/cm]q=200[1/cm]q=250[1/cm]q=300[1/cm] Figure 3: (a) The eigenvalue as a function of the number of independent pores for different valuesof the wave vector q (b) The eigenvalue difference as a function of the number of independentpores for different values of the wave vector q . d (cid:104) m (cid:105) dt = J out N (cid:104) h (cid:105) − (cid:18) J in N + q D in + q σ (cid:19) (cid:104) m (cid:105) . (14)Note that, the wave vector q turns on the stochastic dynamics. We denote δD = D ex − D in , and δσ = σ − σ . The multiplicative white noise accounts for the random diffusivity changes of themedium due to restricted diffusion. Based on the analysis in Sec. 2.2, one can find the transitionpoint in terms of the wave vector q . The transition point is at q c = (cid:113) J δD − δσ . The averagesignal reveals only the transition at q c . Note that, the signal decay is also affected by the noisevariance of the extracellular space. Taking typical values such as, J out = = 0 . ], D ex = 2 e − cm sec ], D in = 0 . e − cm sec ]. Then the critical value is q c ≈ (cid:114) J out D ex − D in + σ = (cid:113) . . e − ≈ ] = 0 . µ m ]. A larger variance in the non-hub pores will set the transition fora lower value of q , value. Whereas larger variance in the diffusion in the extracellular space willset the transition to a higher value of q value. The decay of the signal has a bi-exponential formas predicted by the K¨arger model: γ = lim t,N →∞ t ln ( (cid:104) E ( t ) (cid:105) ) = (cid:40) − q (cid:16) D in + σ (cid:17) − q D ex − J out q > q c q ≤ q c . (15)The stochastic model is a natural generalization of the K¨arger model. Using this generalization,we are able to explore and analyze the behavior of the model in the presence of complex topolog-ical structures, as well as the effect of changes in the apparent diffusivity as a result of stochasticnoise. Note that this transition appears also in the presence of any interaction δ among thenon-hub nodes. Figure 3 presents the behavior of the eigenvalues as a function of N . We have presented a stochastic model that describes diffusion on a graph with an additionalmultiplicative stochastic noise. We analyze this model on a directed graph with one hub node8hat is connected to a large number of non-hub nodes. We derive a non-linear equation for theaverage mass ratio between the non-hub nodes and the hub. This equation describes logisticgrowth. It has two phases one in which the overall mass is mainly concentrated on the hub, andthe other in which there is a “detailed balance” such that the steady-state depends on the ratiobetween the exchange rate between the hub and the non-hub nodes. We show that this modelis completely solvable in the large N limit. In addition, we identify the phase transitions of themodel in terms of the two parameters α and ∆ σ σ . We show that in order for localization phaseoccur the number of non-hub nodes needs to grow at least exponentially with time. Surprisingly,in the limit of large number of independent nodes the transition points do not depend on theamount of resource given by the hub, provided that it is finite and non-zero. We generalize thisanalysis to a system of multiple hubs. We show that in the limit of infinitely many hubs thegrowth of the system becomes exponential.The model has numerous applications. We introduce an application of this model in thecontext of MR measurements of complex structures. Our results in this context may providea theoretical framework that may help interpret and propose new MR experiments to identifythe concentration phases that we see theoretically. This may have impact on the predictionof the underlying measured geometry. Our results and analysis can also be of interest in otherapplications, for example, in predicting economic growth, and in analyzing the stability of controlsystems. Acknowledgment
We would like to thank Prof. Ofer Pasternak for proposing the idea for the paper.
A Transformations of Variable
In this section, we derive the relative magnetization equations, Eq. (4). We use Itˆo’s formula inorder to perform the change of variables df ( t, m ) = ∂f∂t dt + (cid:88) i ∂f∂m i dm i + 12 (cid:88) i,j ∂f ∂m i ∂m j (cid:2) B (cid:3) ij dt = ∂f∂t + (cid:88) i ∂f∂m i A i + 12 (cid:88) i,j ∂f ∂m i ∂m j (cid:2) B (cid:3) ij dt + (cid:88) ij ∂f∂m i B ij dW j , where A and B are the coefficients of the stochastic equations Eq. (2) and Eq. (3), respectively,defined as follows: A i = J i h − J i m i − δ (cid:80) j L ij m j + σ i m i + f i m i , A = (cid:80) j J j m j − (cid:80) j J j h + σ h + f h = J in N (cid:80) j m j − J out h + σ h + f h , and B ij = δ ij σ out m i , and B i = B i = δ i σ m .9 M i dt = J out N − M i (cid:18) J in N + J in M − J out − σ f (cid:19) − δ (cid:88) j L ij M j + (cid:113) σ + σ ξ i M i = J out N − M i (cid:18) J in N + J in M − J out − σ f (cid:19) − δ (cid:88) j L ij M j + σξ i M i = − M i σ (cid:16) βM − α − − ˜ δ (cid:17) − δM ( t ) + σξ i M i . We introduce the dimensionless parameters α = 2 J out − ∆ fσ , and β = J in σ . The equations arewritten under the assumption that the interaction among the nodes and the hub is describedby J i = J out N , and J j = J in N , respectively. We also assume, for simplicity, that all the non-hub compartments obey the statistics σ i = σ out and f i = f . The last transition is under theassumption of mean-field interaction δ among the non-hub nodes. In the Itˆo form in the limit of N → ∞ , we havelim N →∞ N (cid:88) i σ i g i ( t ) M i = 0 . Note that one cannot take the limit N → ∞ in Eq. (3), since the variables m i depend on thestochastic noise of the hub, g . Taking the sum and letting N → ∞ , we arrive to the followingdeterministic equation describing the growth of the average relative mass of the N non-hub nodesas a function of time: dM ( t ) dt = M ( t ) (cid:18) J out + σ − ∆ f (cid:19) − J in M ( t ) . (16)The steady-state solutions of this non-linear equation are: M → − ∆ f − J out − σ J in = αβ as N → ∞ , and M → N → ∞ . The equation is also valid when δ (cid:54) = 0. A.1 Finite- N corrections In this subsection, we consider finite N corrections to the average equation. With this effectaccumulated for the equation reads dM ( t ) dt = εa + εbM ( t ) + cM ( t ) + bM ( t ) , (17)where ε = N . We denote a = σ J out , b = − σ β , and c = σ ( α + 1). The equation has a Riccatiform. Taking the first-order correction in ε , M ( t ) = M ( t ) + εM ( t ), we get dM ( t ) dt = cM ( t ) + bM ( t ) . (18)The solution for M ( t ) is the logistic function as before. Next, dM ( t ) dt = a + bM ( t ) + ( c + 2 bM ( t )) M ( t ) . (19)10he solution for M ( t ) is given by M ( t ) = M (0)exp( (cid:90) t ( c + 2 bM ( s )) ds ) + (cid:90) t exp( (cid:90) ts ( c + 2 bM ( s )) ds )( a + bM ( s )) ds = M (0)exp( ct + 2 b (cid:90) t M ( s ) ds ) + (cid:90) t exp( c ( t − s ) + 2 b (cid:90) ts M ( τ ) dτ )( a + bM ( s )) ds. Substituting the expression of the solution to M ( t ), one can show that lim t →∞ M ( t ) is finite,meaning that the correction of order N to Eq. (16) is negligible in the large N limit. B Moments Lyapunov Exponent
In this section, we derive the moment Lyapunov exponent Eq. (9). For this purpose, we presentthe first moment equations of Eq. (2) and Eq. (3) for a finite number of non-hub nodes N and δ = 0. These equations can be derived using the Fokker-Plank equation or alternatively theFeynman-Kac formula [8, 24, 16, 15]. They read d (cid:104) h (cid:105) dt = J in N (cid:104) m (cid:105) − J out (cid:104) h (cid:105) + f (cid:104) h (cid:105) + σ (cid:104) h (cid:105) , (20) d (cid:104) m i (cid:105) dt = J out N (cid:104) h (cid:105) − J in N (cid:104) m i (cid:105) + f (cid:104) m i (cid:105) + σ (cid:104) m i (cid:105) . (21)Here it is assumed that all the non-hub nodes are independent and with the same dynamics anddenoting the total mass, (cid:104) m (cid:105) = N (cid:104) m i (cid:105) . Summing over N in Eq. (21), we get d (cid:104) m (cid:105) dt = J out (cid:104) h (cid:105) − J in N (cid:104) m (cid:105) + f (cid:104) m (cid:105) + σ (cid:104) m (cid:105) . (22)Eqs. (22) and (20) can be written in vector form as d a dt = A a ( t ) = (cid:32) f − J out + σ J in N J out f − J in N + σ (cid:33) a ( t ) . This is a simple system of linear equations, and the dynamics it describes is determined by theeigenvalues of the matrix A . The resulted eigenvalues are than, λ = f + σ − J in N (cid:32) J out (∆ f + ∆ σ − J in N − J out ) (cid:33) , (23)and, λ = f − J out + σ J out J in N (∆ f + ∆ σ − J in N − J out ) . (24)11he solution is then given by (cid:18) (cid:104) m (cid:105) ( t ) (cid:104) m (cid:105) ( t ) (cid:19) = Av e λ t + Bv e λ t = A (cid:18) λ − dc (cid:19) e λ t + B (cid:18) λ − dc (cid:19) e λ t = A (cid:18) ∆ f − J out + ∆ σ J out (cid:19) e λ t + B (cid:18) J out (cid:19) e λ t . Plugging this expression into the equation ofthe first moment Lyapunov exponent we get γ = lim t →∞ lim N →∞ t ln ( (cid:104) h (cid:105) + (cid:104) m (cid:105) )= lim t →∞ t ln e λ t (cid:18) A (cid:18) ∆ f − J + ∆ σ J out (cid:19) e ( λ − λ ) t + BJ out (cid:19) = (cid:40) f + σ f − J out + σ σ + ∆ f − J out < ∆ σ + ∆ f − J out ≥ (cid:40) f + σ f − J out + σ α > ∆ σ σ α ≤ ∆ σ σ . C Sample Lyapunov Exponent
Here, we calculate the sample Lyapunov exponent of the mass of the system, defined as (cid:101) γ = lim t,N →∞ t ln (cid:32) h + (cid:88) i m i (cid:33) . To prove the formula (10) in the main text, we consider at the system of equations dm i dt = J out N h − J in N m i + f m i + σ m i + σ out g i ( t ) m i ,dh dt = J in N m + (cid:18) f + σ − J out (cid:19) h + σ g ( t ) h . Using a transformation of variable given by the Itˆo formula12 ln ( m + h ) dt = ∂f∂t + (cid:88) i ∂f∂m i A i + ∂f∂h A + 12 ∂f ∂ h (cid:2) B (cid:3) + 12 (cid:88) i ∂f ∂ m i (cid:2) B (cid:3) ii + ∂f∂h B dW dt + (cid:88) i ∂f∂m i B ii dW i dt = 1 N M + 1 (cid:32) J out − J in M + f M + σ M + σ out (cid:88) i g i ( t ) M i (cid:33) + 1 N M + 1 (cid:18) J in M + f + σ − J out (cid:19) − σ + σ (cid:80) i M i (cid:0) N M + 1 (cid:1) + 1 N M + 1 σ g ( t ) + σ out (cid:88) i g i ( t ) M i N M + 1= 1
N M + 1 σ g ( t ) + σ out (cid:88) i g i ( t ) M i N M + 1+ 1
N M + 1 (cid:18)(cid:18) f + σ (cid:19) N M + (cid:18) f + σ (cid:19)(cid:19) − σ + σ (cid:80) i M i (cid:0) N M + 1 (cid:1) . Integrating with respect to time and dividing by t , we have1 t ln ( m ( t ) + h ( t )) = 1 t (cid:90) t d ln ( m ( τ ) + h ( τ )) dτ dτ. Letting the limit of t, N → ∞ , we get (cid:101) γ = lim t,N →∞ t ln ( m + h ) = lim t,N →∞ t (cid:90) t d ln ( m ( τ ) + h ( τ )) dτ dτ = lim t,N →∞ t (cid:90) t (cid:34) σ g ( τ ) N M + 1 − σ (cid:0) N M + 1 (cid:1) + σ out (cid:88) i g i ( τ ) M i N M + 1 − σ (cid:80) i M i (cid:0) N M + 1 (cid:1) (cid:35) dτ + lim t,N →∞ t (cid:90) t N M + 1 (cid:18)(cid:18) f + σ (cid:19) N M + (cid:18) f + σ (cid:19)(cid:19) dτ. (25)Here we used the deterministic law of M given that higher corrections in N are negligible, seeA.1 for more details. Using the stationarity property of the Gaussian processes g ( t ) and g i ( t ),and the ergodic theorem: (cid:101) γ = lim t,N →∞ t (cid:90) t (cid:34) (cid:0) N M + 1 (cid:1) (cid:18)(cid:18) f + σ (cid:19) N M + (cid:18) f + σ (cid:19)(cid:19)(cid:35) dτ = lim t,N →∞ (cid:18) f + σ (cid:19) t (cid:90) t MM + N dτ + (cid:18) f + σ (cid:19) t (cid:90) t dτN M + 1 . (26)This formula is obtained under the assumption that the fluctuations of the noise are finite13Novikov condition), i.e.,lim t,N →∞ t (cid:90) t (cid:80) i M i (cid:0) N M + 1 (cid:1) dτ < ∞ (27)and, lim t,N →∞ t (cid:90) t (cid:0) N M + 1 (cid:1) dτ < ∞ . (28)In order to calculate the above integrals one needs to specify how the number of nodes in thegraphs growth with time. We show that if the number of nodes grows exponentially in t , i.e., N ∼ e εt , then there exists a localization phase. In addition, the fluctuation, conditions 27 and28, are satisfied. These conditions are verified in C.1 below. Since the fluctuations are finite, weare left with calculating integrals of the logistic function M ( t ). We use the following propertiesof the logistic function: M ( t ) = 1 A + Be − ξt = e ξt Ae ξt + B (29)and (cid:90) t M dτ = 1 ξA log (cid:0) Ae ξt + B (cid:1) . (30)In our case, A = βα +1 , B = M − βα +1 , and ξ = σ ( α +1)2 . Using Eqs. (29), and (30) it is easy toshow thatlim t,N →∞ t (cid:90) t MM + N dτ = lim t →∞ t (cid:90) tT MM + e − ε ( τ − T ) dτ = lim t →∞ t (cid:90) tT
11 + Ae εT e − ετ + Be εT e − ( ξ + ε ) τ dτ = 1 , (31)where T is some finite time. In the same manner, one can calculate the second term in Eq. (26):lim t,N →∞ t (cid:90) t dτN M + 1 = lim t →∞ t (cid:90) t A + Be − ξτ dτN + A + Be − ξτ = lim t →∞ t (cid:90) t Adτe ε ( τ − T ) + A + Be − ξτ + lim t →∞ t (cid:90) t Be − ξτ e ε ( τ − T ) + A + Be − ξτ dτ = lim t →∞ t (cid:90) t Be − ( ξ + ε ) τ e − εT + Be − ( ξ + ε ) τ dτ = lim t →∞ − t ( ξ + ε ) log (cid:16) e − εT + Be − ( ξ + ε ) t (cid:17) = 0 , (32)provided that max {− ξ, } ≤ ε . Note that, in order to have finite fluctuations, σ ≤ (cid:15) , see C.1.Substituting the results in Eq. (31) and Eq. (32) we obtain (cid:101) γ = f + σ . (33)Therefore, for any exponential growth rate σ max {− α +12 , } ≤ ε , the fluctuations and sampleLyapunov are finite. 14 .1 Noise fluctuations In this subsection, we prove that the fluctuation are finite when N ∼ e εt , i.e., the conditions in(27) and (28) are satisfied. That condition (28) is satisfied, is readily seen since the functionunder the integral is bounded between zero and one, and so the integral itself is also finite:lim t,N →∞ t (cid:90) t dτ (cid:0) N M + 1 (cid:1) < ∞ In order to calculate the fluctuation of the non-hub nodes, i.e., establish (27), we approximatethe limit, lim N →∞ N (cid:80) i M i → (cid:104) M i (cid:105) ; since M i are i.i.d. random variables. The second momentcan be calculated using the Fokker-Plank equation Eq. (7): ∂P ( M, t ) ∂t = − σ ∂ (cid:0)(cid:0)(cid:0) α + 1 − βM (cid:1) M (cid:1) P ( M, t ) (cid:1) ∂M + σ ∂ (cid:0) M P ( M, t ) (cid:1) ∂ M .
Averaging over the second moment yields: d (cid:104) M (cid:105) dt = σ (cid:0) (cid:0) α + 2 − βM (cid:1) (cid:104) M (cid:105) (cid:1) . The solution for ˜ δ = 0 is (cid:104) M ( t ) (cid:105) = (cid:104) M (cid:105) (0)exp (cid:18) σ t ( α + 2) − σ β (cid:90) t M dτ (cid:19) = (cid:104) M (cid:105) (0) e σ t (cid:18) β ( α +1) + β ( α +1 β − M ) M ( α +1) e − σ α +1) t (cid:19) = (cid:104) M (cid:105) (0) e σ t M ( t ) . The fluctuations of the non-hub nodes, i.e., the last term in Eq. (25) are then, given bylim t,N →∞ t (cid:90) t (cid:104) M i (cid:105) N (cid:0) M + N (cid:1) dτ = lim t,N →∞ t (cid:90) t (cid:104) M i (cid:105) N (cid:0) M + N (cid:1) dτ = (cid:104) M (cid:105) (0) lim t,N →∞ t (cid:90) t e σ t M N (cid:0) M + N (cid:1) dτ. Substituting here the logistic function M ( t ) (Eq. (29)), we getlim t,N →∞ t (cid:90) t e σ τ M N (cid:0) M + N (cid:1) dτ = lim t,N →∞ t (cid:90) t e σ τ N + 2 ( A + Be − ξτ ) + N ( A + Be − ξτ ) dτ = lim t →∞ t (cid:90) t e σ τ e ε ( τ − T ) + 2 ( A + Be − ξτ ) + e − ε ( τ − T ) ( A + Be − ξτ ) dτ = lim t →∞ t ( σ − ε ) e ( σ − ε ) ( t − T )+ εT dτ = 0 . Therefore, the fluctuation are finite for any ε ≥ σ .15 Multiple Hubs Derivation