From Logistic Growth to Exponential Growth in a Population Dynamical Model

TL;DR

This paper analyzes a stochastic population model with resource dynamics, deriving an exact solution and phase diagram that reveal three regimes of resource distribution, with applications in porous media analysis.

Contribution

It provides an exact solution to a resource distribution model with multiplicative noise and characterizes three distinct dynamical regimes.

Findings

Derives a deterministic equation for the average mass ratio showing logistic growth.

Identifies three regimes of resource distribution among components.

Transition points are independent of resource amount in large systems.

Abstract

Dynamics among central sources (hubs) providing a resource and large number of components enjoying and contributing to this resource describes many real life situations. Modeling, controlling, and balancing this dynamics is a general problem that arises in many scientific disciplines. We analyze a stochastic dynamical system exhibiting this dynamics with a multiplicative noise. We show that this model can be solved exactly by passing to variables that describe the mass ratio between the components and the hub. We derive a deterministic 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 calculating the sample and moment Lyapunov exponent of the system. The first regime describes full balance between the non-hub components and the hub, in the second regime the entire resource is concentrated mainly in the hub, and in the third regime the resource is localized on a few non-hub components and the hub. Surprisingly, in the limit of large number of components the transition values do not depend on the amount of resource given by the hub. This model has interesting application in the context of analysis of porous media using Magnetic Resonance (MR) techniques.

Equations

$$\frac{d{m}_i(t)}{dt}=\sum _j{J}_{ij}{m}_j(t)-\sum _j{J}_{ij}{m}_i(t)+g_i(t)m_i(t),$$

$$\frac{d{h}_0}{dt}=\frac{J_{\mathrm{in}}}{N}\sum _j{m}_j-J_{\mathrm{out}}{h}_0+f_0{h}_0+\frac{{\sigma }_0^2}{2}{h}_0+{\sigma }_0{g}_0{h}_0,$$

$$\frac{d{m}_i}{dt}=\frac{J_{\mathrm{out}}}{N}{h}_0-\frac{J_{\mathrm{in}}}{N}{m}_i+f{m}_i+\frac{{\sigma }_{\mathrm{out}}^2}{2}{m}_i-\delta \sum _j{L}_{ij}{m}_j+{\sigma }_{\mathrm{out}}{g}_i{m}_i.$$

$$\frac{d{M}_i}{dt}=\frac{J_{\mathrm{out}}}{N}-\left(\frac{J_{\mathrm{in}}}{N}+\Delta f-J_{\mathrm{out}}-\frac{{\sigma }^2}{2}+J_{\mathrm{in}}\overline{M}(t)\right)M_i-\delta \sum _j{L}_{ij}{M}_j+\sigma {\xi }_i{M}_i.$$

$$\frac{d\overline{M}}{dt}=\frac{{\sigma }^2\beta }{2}\overline{M}\left(\frac{\alpha +1}{\beta }-\overline{M}\right),$$

$$\overline{M}(t)=\frac{{\overline{M}}_0(\alpha +1)}{\beta \left({\overline{M}}_0+\left(\frac{\alpha +1}{\beta }-{\overline{M}}_0\right)e^{-\frac{{\sigma }^2\left(\alpha +1\right)t}{2}}\right)}.$$

$$\frac{\partial P}{\partial t}=-\frac{{\sigma }^2}{2}\frac{\partial \left(\left(\left(\alpha +1-\beta \overline{M}-\tilde{\delta }\right)M+\tilde{\delta }\overline{M}(t)-\tilde{\delta }\right)P\right)}{\partial M}+\frac{{\sigma }^2}{2}\frac{{\partial }^2\left(M^{2}P\right)}{{\partial }^{2}M}.$$

$$P_{\mathrm{eq}}(M)=A\exp \left(-\frac{\tilde{\delta }{\overline{M}}_{\mathrm{eq}}}{M}\right)M^{-\mu ({\overline{M}}_{\mathrm{eq}})}.$$

$$\gamma_{1}=\underset{t,N\rightarrow\infty}{\mathrm{lim}}\,\frac{1}{t}\mathrm{ln}\left(\langle E(t)\rangle\right)=\begin{cases}\begin{array}[]{c}f+\frac{\sigma_{\mathrm{out}}^{2}}{2}\\ f_{0}-J_{\mathrm{out}}+\frac{\sigma_{0}^{2}}{2}\end{array}&\begin{array}[]{c}\frac{\Delta\sigma^{2}}{\sigma^{2}}<\alpha\\ \frac{\Delta\sigma^{2}}{\sigma^{2}}\geq\alpha\end{array}\end{cases},$$

$$\tilde{\gamma }=\underset{t,N\to \infty }{\lim }\frac{1}{t}\ln \left(h_0+\overline{m}\right)=f+\frac{{\sigma }_{\mathrm{out}}^2}{2},$$

$$\frac{d{h}_i}{dt}=\frac{J_{\mathrm{in}}}{HN}\sum _i{m}_i-\frac{J_{\mathrm{out}}}{H}{h}_i+f_i{h}_i+\frac{\delta }{H-1}\sum _{j\ne i}^H{h}_j-\delta h_i+{\sigma }_i{g}_i^h{h}_i,$$

$$\frac{d{m}_i}{dt}=\frac{J_{\mathrm{out}}}{NH}\sum _i{h}_i-\frac{J_{\mathrm{in}}}{N}{m}_i+q_i{m}_i+{\nu }_i{g}_i^m{m}_i.$$

$$\frac{d⟨h_0⟩}{dt}=\frac{J_{\mathrm{in}}}{N}⟨\overline{m}⟩-J_{\mathrm{out}}⟨h_0⟩-q^2{D}_{\mathrm{ex}}⟨h_0⟩,$$

$$\frac{d⟨\overline{m}⟩}{dt}=\frac{J_{\mathrm{out}}}{N}⟨h_0⟩-\left(\frac{J_{\mathrm{in}}}{N}+q^2{D}_{\mathrm{in}}+\frac{q^2{\sigma }_{\mathrm{in}}^2}{2}\right)⟨\overline{m}⟩.$$

$$\gamma_{1}=\underset{t,N\rightarrow\infty}{\mathrm{lim}}\,\frac{1}{t}\mathrm{ln}\left(\langle E(t)\rangle\right)=\begin{cases}\begin{array}[]{c}-q^{2}\left(D_{\mathrm{in}}+\frac{\sigma_{\mathrm{in}}^{2}}{2}\right)\\ -q^{2}D_{\mathrm{ex}}-J_{\mathrm{out}}\end{array}&\begin{array}[]{c}q>q_{\mathrm{c}}\\ q\leq q_{\mathrm{c}}\end{array}\end{cases}.$$

$$\begin{gathered} df(t,\mathbf{m})=\frac{\partial f}{\partial t}dt+\sum _i\frac{\partial f}{\partial m_i}d{m}_i+\frac{1}{2}\sum _{i,j}\frac{\partial f^2}{\partial m_i\partial m_j}{\left[B^2\right]}_{ij}dt \\ =\left(\frac{\partial f}{\partial t}+\sum _i\frac{\partial f}{\partial m_i}{A}_i+\frac{1}{2}\sum _{i,j}\frac{\partial f^2}{\partial m_i\partial m_j}{\left[B^2\right]}_{ij}\right)dt+\sum _{ij}\frac{\partial f}{\partial m_i}{B}_{ij}d{W}_j, \end{gathered}$$

$$\begin{gathered} \frac{d{M}_i}{dt}=\frac{J_{\mathrm{out}}}{N}-M_i\left(\frac{J_{\mathrm{in}}}{N}+J_{\mathrm{in}}\overline{M}-J_{\mathrm{out}}-\frac{{\sigma }^2}{2}+\Delta f\right)-\delta \sum _j{L}_{ij}{M}_j+\sqrt{{\sigma }_{\mathrm{out}}^2+{\sigma }_0^2}{\xi }_i{M}_i \\ =\frac{J_{\mathrm{out}}}{N}-M_i\left(\frac{J_{\mathrm{in}}}{N}+J_{\mathrm{in}}\overline{M}-J_{\mathrm{out}}-\frac{{\sigma }^2}{2}+\Delta f\right)-\delta \sum _j{L}_{ij}{M}_j+\sigma {\xi }_i{M}_i \\ =-M_i\frac{{\sigma }^2}{2}\left(\beta \overline{M}-\alpha -1-\tilde{\delta }\right)-\delta \overline{M}(t)+\sigma {\xi }_i{M}_i. \end{gathered}$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\sum _i{\sigma }_i{g}_i(t)M_i=0.$$

$$\frac{d\overline{M}(t)}{dt}=\overline{M}(t)\left(J_{\mathrm{out}}+\frac{{\sigma }^2}{2}-\Delta f\right)-J_{\mathrm{in}}\overline{M}(t{)}^2.$$

$$\frac{d\overline{M}(t)}{dt}=\epsilon a+\epsilon b\overline{M}(t)+c\overline{M}(t)+b\overline{M}(t{)}^2,$$

$$\frac{d{\overline{M}}_0(t)}{dt}=c{\overline{M}}_0(t)+b{\overline{M}}_0(t{)}^2.$$

$$\frac{d{\overline{M}}_1(t)}{dt}=a+b{\overline{M}}_0(t)+(c+2b{\overline{M}}_0(t)){\overline{M}}_1(t).$$

$$\begin{gathered} {\overline{M}}_1(t)={\overline{M}}_1(0)\exp ({\int }_0^t(c+2b{\overline{M}}_0(s))ds)+{\int }_0^t\exp ({\int }_s^t(c+2b{\overline{M}}_0(s))ds)(a+b{\overline{M}}_0(s))ds \\ ={\overline{M}}_1(0)\exp (ct+2b{\int }_0^t{\overline{M}}_0(s)ds)+{\int }_0^t\exp (c(t-s)+2b{\int }_s^t{\overline{M}}_0(\tau )d\tau )(a+b{\overline{M}}_0(s))ds. \end{gathered}$$

$$\frac{d⟨h_0⟩}{dt}=\frac{J_{\mathrm{in}}}{N}⟨\overline{m}⟩-J_{\mathrm{out}}⟨h_0⟩+f_0⟨h_0⟩+\frac{{\sigma }_0^2}{2}⟨h_0⟩,$$

$$\frac{d⟨m_i⟩}{dt}=\frac{J_{\mathrm{out}}}{N}⟨h_0⟩-\frac{J_{\mathrm{in}}}{N}⟨m_i⟩+f⟨m_i⟩+\frac{{\sigma }_{\mathrm{out}}^2}{2}⟨m_i⟩.$$

$$\frac{d⟨\overline{m}⟩}{dt}=J_{\mathrm{out}}⟨h_0⟩-\frac{J_{\mathrm{in}}}{N}⟨\overline{m}⟩+f⟨\overline{m}⟩+\frac{{\sigma }_{\mathrm{out}}^2}{2}⟨\overline{m}⟩.$$

$$\frac{d\boldsymbol{a}}{dt}=A\boldsymbol{a}(t)=\left(\begin{array}[]{cc}f_{0}-J_{\mathrm{out}}+\frac{\sigma_{0}^{2}}{2}&\frac{J_{\mathrm{in}}}{N}\\ J_{\mathrm{out}}&f-\frac{J_{\mathrm{in}}}{N}+\frac{\sigma_{\mathrm{out}}^{2}}{2}\end{array}\right)\boldsymbol{a}(t).$$

$${\lambda }_1=f+\frac{{\sigma }_{\mathrm{out}}^2}{2}-\frac{2J_{\mathrm{in}}}{N}\left(1+\frac{J_{\mathrm{out}}}{(\Delta f+\frac{\Delta {\sigma }^2}{2}-\frac{J_{\mathrm{in}}}{N}-J_{\mathrm{out}})}\right),$$

$${\lambda }_2=f_0-J_{\mathrm{out}}+\frac{{\sigma }_0^2}{2}+\frac{2J_{\mathrm{out}}\frac{J_{\mathrm{in}}}{N}}{(\Delta f+\frac{\Delta {\sigma }^2}{2}-\frac{J_{\mathrm{in}}}{N}-J_{\mathrm{out}})}.$$

$$\gamma_{1}=\underset{t\rightarrow\infty}{\mathrm{lim}}\,\underset{N\rightarrow\infty}{\mathrm{lim}}\,\frac{1}{t}\mathrm{ln}\left(\langle h_{0}\rangle+\langle\overline{m}\rangle\right)\\ =\mathrm{\underset{t\rightarrow\infty}{lim}}\,\frac{1}{t}\mathrm{ln}\,e^{\lambda_{2}t}\left(A\left(\Delta f-J_{0}+\frac{\Delta\sigma^{2}}{2}+J_{\mathrm{out}}\right)e^{\left(\lambda_{1}-\lambda_{2}\right)t}+BJ_{\mathrm{out}}\right)\\ =\begin{cases}\begin{array}[]{c}f+\frac{\sigma_{\mathrm{out}}^{2}}{2}\\ f_{0}-J_{\mathrm{out}}+\frac{\sigma_{0}^{2}}{2}\end{array}&\begin{array}[]{c}\frac{\Delta\sigma^{2}}{2}+\Delta f-J_{\mathrm{out}}<0\\ \frac{\Delta\sigma^{2}}{2}+\Delta f-J_{\mathrm{out}}\geq 0\end{array}\end{cases}\\ =\begin{cases}\begin{array}[]{c}f+\frac{\sigma_{\mathrm{out}}^{2}}{2}\\ f_{0}-J_{\mathrm{out}}+\frac{\sigma_{0}^{2}}{2}\end{array}&\begin{array}[]{c}\alpha>\frac{\Delta\sigma^{2}}{\sigma^{2}}\\ \alpha\leq\frac{\Delta\sigma^{2}}{\sigma^{2}}\end{array}\end{cases}.$$

Full text

1 Introduction

Population dynamics on large scale networks has attracted a lot of attention due to its wide occurrence in many disciplines, such as social sciences [1, 2], physics [3] and biology, communication and control theory [4]. This dynamics is mainly affected by the topology of the network as well as some internal stochastic noise. In many applications there are only a few nodes playing a major role in the dynamical process, distributing and carrying most of the resources [5, 6]. For example, 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 context of 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 water molecules can be utilized to extract information about the network of cells (neurons) in the brain. The concept of self-diffusion of molecules in a network of pores was already introduced in Refs. [8, 9, 10]. The main challenge is how to determine the topology of the network based on the MR measurements [8].

Our model consists of a system of interacting sites on a graph $\mathcal{G}$ with $N$ vertices and $E$ edges between them. We are interested in the stochastic dynamics of some characteristic property \{m_{i}(t)\}_{i\in\mathcal{G},t\text{\geq}0}. The property $m_{i}(t)$ is linked to a physical measurable quantity in the real world and the graph is the underlying geometry/topology in which the property lives. The topology is a complex network of sites. The model is described by the following family of stochastic differential equations in the Stratonovich form on the graph $\mathcal{G}$:

[TABLE]

with the initial conditions $m_{i}(0)=m_{0}$. The term $g_{i}(t)$ is a multiplicative white noise, such that $\langle g_{i}(t)\rangle=f_{i}$, and $\langle g_{i}(t)g_{j}(t^{\prime})\rangle=\sigma_{i}^{2}\delta_{ij}\delta(t-t^{\prime})$. We choose the Stratonovich form, since its 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. The model consists of two parts: an interacting part, where the interaction strength depends on the location on the graph, and a non-interacting part, where each component follows a stochastic noise with different variance $\sigma_{i}^{2}$. The first part causes spreading, while the second pushes towards concentration (a.k.a localization or condensation). The model was already analyzed in the mean field 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 on the lattice is known in the mathematical literature as the time-dependent Parabolic Anderson Model (PAM) [15, 16]. The phase diagram of the model in this case depends on the dimension of the lattice. On a general network, phase transitions depend on the spectral dimension of the network [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, the MR signal measured is assumed to be composed of two contributions: one coming from hindered diffusion in the extracellular space and the other from restricted diffusion in the intracellular space [17, 18, 19, 20]. The hub node represents the magnetization in the extracellular space (e.g., water), $h_{0}(t)$, and the non-hub nodes represents $N$ independent intracellular pores with magnetization, $m_{i}(t)$. The motion of molecules between these regions changes the value of the magnetization 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 $\sigma_{i}^{2}$. In the economic context, the system describes the dynamics of the 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. The stochastic noise represents the bank/state and the agent’s investments in the stock market and housing [1, 7]. Analysis of the dynamics of the sums of the money held by the agents and the bank/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 of nodes growth at least exponentially with time there exist a localization phase. We also generalize our analysis to multiple number of hubs.

Our main result is a full phase-diagram of the model. We show that this model can be described by a stochastic equation for the mass ratio between each of the non-hub nodes and the hub, and a deterministic non-linear equation for the average relative mass of all the non-hub nodes with respect to the hub. To identify the phases of the system, we calculate the sample and moment Lyapunov exponents and identify a gap between them. The phase transitions are characterized by one parameter. This parameter takes into account the exchange rate between the non-hub nodes and the hub and the variances of the multiplicative noises.

2 Hub Topology

The basic hub topology is composed of an infinite number of nodes, $\{m_{i}\}$, interacting at constant rate with a hub node, $h_{0}$, such that, $J_{i0}=\frac{J_{\mathrm{out}}}{N}$, and $J_{0j}=\frac{J_{0}}{N}$, respectively. Our normalization is such that, the overall interaction between the nodes and the hub is finite in the limit of infinite number of nodes. The interaction among the non-hub nodes is characterized by the parameter $\delta$; when $\delta=0$, any interaction (transfer of mass) between the non-hub nodes is done only through the hub. The topology of the interaction between the non-hub nodes is defined by a Laplacian matrix, $L$, satisfying $\sum_{i}L_{ij}=0$. Figure 1 presents an illustration of such a system for $\delta=0$.

We assume that the stochastic noise acting on the non-hub nodes has the same variance and average for all nodes, $\sigma_{i}=\sigma_{\mathrm{out}}$, $f_{i}=f$. Eq. (1) in the Itô form reduces to the following system of stochastic equations:

[TABLE]

[TABLE]

The Itô form will be useful latter on when one takes the limit $N\rightarrow\infty$. It is instructive to pass to the following normalized variables by introducing the mass ratio between the non-hub nodes and the hub node, $M_{i}=\frac{m_{i}}{h_{0}}$. This leads to the following system of $N$ equations in the Itô form (See A for more details):

[TABLE]

We introduce the average mass ratio $\overline{M}=\frac{1}{N}\sum_{i}M_{i}$, the effective variance $\frac{\sigma^{2}}{2}=\frac{\sigma_{0}^{2}+\sigma_{\mathrm{out}}^{2}}{2}$, and the average difference $\Delta f=f_{0}-f$. Here, $\xi_{i}(t)$ is a Gaussian process with mean zero and variance one. The transformation of variables introduces a non-linear term, which accounts for the interaction between any two nodes through the hub. In the limit $N\rightarrow\infty$, averaging over all the nodes in Eq. (4) yields the following deterministic law for the average mass ratio:

[TABLE]

where we use the dimensionless parameters $\alpha=2\frac{J_{\mathrm{out}}-\Delta f}{\sigma^{2}}$ and $\beta=\frac{2J_{\mathrm{in}}}{\sigma^{2}}$. Since, Eq. (5) it is a non-linear equation, there are two steady-state solutions: $\overline{M}_{1\mathrm{eq}}{\rightarrow}\frac{J_{\mathrm{out}}-\Delta f+\frac{\sigma^{2}}{2}}{J_{\mathrm{in}}}=\frac{\alpha+1}{\beta}$ as $N\rightarrow\infty$, and $\overline{M}_{2\mathrm{eq}}{\rightarrow}0$ as $N\rightarrow\infty$. Note that, in the case of $\delta=0$ the system geometry can be viewed as a directed tree structure with one level and infinitely many leaves. Therefore, the presence of a the non-linear term is caused by to the indirect interaction between the non-hub nodes following the tree topology of the system [12]. Convergence to each one of these fixed points depend on the initial condition of the system and on the stability of these points. Stability analysis of these two points shows that the point $\overline{M}_{1\mathrm{eq}}$ is stable when $\alpha+1>0$, while the point $\overline{M}_{2\mathrm{eq}}=0$, is stable when $\alpha+1\leq 0$. For example, in the context of MR measurements in porous media, this system can model a complex structure measured from a single voxel in the MR image. The value $\overline{M}_{\mathrm{eq}}$ describes the steady-state average magnetization ratio, between the extracellular space and the intracellular space. The first fixed point $\overline{M}_{1\mathrm{eq}}$ is reached when the steady-state magnetization ratio is equal to the amount of molecules leaving the non-hub pores reduced by the magnetic field effects divided by the amount of molecules leaving the extracellular space. The second fixed point $\overline{M}_{2\mathrm{eq}}$ represents the case in which on average most of the contribution to the magnetization in a single voxel comes mainly from the extracellular space. Eq. (5) shows logistic growth and is a version of Lotka-Volterra equation, which describes many social phenomena in nature [21, 22]. It can be solved exactly: for an initial condition $\overline{M}(0)=\overline{M}_{0}$, the solution is

[TABLE]

2.1 Equilibrium Distribution

Given Eqs. (4) and (5) for the relative mass between the non-hub nodes and the hub, one can 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 is governed by a deterministic non-linear equation, in the limit $N\rightarrow\infty$ the system reduces to a set of stochastic independent equations for the relative mass in each node, $M_{i}$. Therefore, we can 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-Plank equation:

[TABLE]

To study the effect of interaction between nodes, we introduce a small mean field interaction between the non-hub nodes, $\delta$, and define the dimensionless interaction rate, $\tilde{\delta}=\frac{2\delta}{\sigma^{2}}$. By equating the left hand side to of Eq. (7) zero one can find the steady-state distribution. The solution shows a Pareto power-law behavior:

[TABLE]

The power is a function of the steady-state average relative mass $\overline{M}_{\mathrm{eq}}$, i.e., steady-state solution of Eq. (5): $\mu(\overline{M}_{\mathrm{eq}})=\beta\overline{M}_{\mathrm{eq}}+1-\alpha+2\tilde{\delta}$. Substituting the value of the average mass ratio, we find that $\mu(\frac{\alpha+1}{\beta})|_{\alpha+1\geq 0}=2+2\tilde{\delta}$, and $\mu(0)|_{\alpha+1<0}\rightarrow 1-\alpha+2\tilde{\delta}$. This shows that the system has 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 $\mu$ is greater than $2$ when the system collapses to the state $\overline{M}_{\mathrm{eq}}=\frac{\alpha+1}{\beta}$, showing equality among the non-hub nodes and the hub. For, $\overline{M}_{\mathrm{eq}}=0$, the mass is localized on a few nodes within the set of non-hub nodes. The power is, $\mu<2$, as long as $2\tilde{\delta}<1+\alpha$. Therefore, adding interaction between non-hub nodes reduces localization, as expected.

2.2 Balance and Localization

The analysis above reveals the localization regime within the non-hub nodes when the influence of the hub is renormalized. In this section, we analyze the regime at which there is localization on the hub. For this purpose, we study the asymptotic properties of the total mass, $E(t)=h_{0}(t)+\sum_{i}m_{i}(t)$. We calculate the Lyapunov exponents of the solution [15, 16]. Here, we perform the analysis for the case of $\delta=0.$ The Lyapunov exponents describe the growth rates of the 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:

[TABLE]

where $\Delta\sigma^{2}=\sigma_{0}^{2}-\sigma_{\mathrm{out}}^{2}$ is the variance difference, see B for details of the proof together with corrections for finite network size. Interestingly, what determines the growth on average is the difference between the variance, $\Delta\sigma^{2}$ of the stochastic noises of the hub and that of the non-hub nodes. To understand Eq. (9), let us look at the limit where the stochastic noise in the non-hub nodes has significantly higher variance compared to the hub. This is equivalent to the presence of large fluctuations in the non-hub pores with respect to the extracellular space, i.e., $\Delta\sigma^{2}\approx-\sigma_{\mathrm{out}}^{2}$. In this case, comparing to the stability points in Eq. (5) in the regime $\alpha\leq-1$, there is a high concentration of magnetization on the hub and a few non-hub nodes which contribute most to the total growth rate. For $\alpha>-1$, in the detailed balance limit the non-hub nodes and the hub contribute equally to the total growth. On the other hand, in the limit of $\Delta\sigma^{2}\approx\sigma_{0}^{2}$, the system exhibits three regimes: for $\alpha>1$ the magnetization spreads equally among the nodes, for $-1<\alpha\leq 1$ the magnetization is mainly concentrated on the hub, and for $\alpha<-1$, there is concentration of the magnetization on the hub and/or several non-hub nodes. This analysis is verified by calculating the value of the sample Lyapunov exponent, which provide the sample growth rate.

We define the sample Lyapunov exponent, $\widetilde{\gamma}$, as the limit of the logarithm of the total mass of the solution divided by time as $N,t\rightarrow\infty$. Knowing the dynamics of the average mass ratio, see Eq. (5), we can calculate the sample growth rate of the mass ratio exactly. The resulting sample/quenched Lyapunov exponent is

[TABLE]

where we denote $\overline{m}=\sum_{i}m_{i}$; see C for details of the proof. Note that, in order to take the limit one needs to specify at which rate the number of nodes growth with time. We show that when the number 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 first moment Lyapunov exponent, $\gamma_{1}$. Localization of the solution is defined as the regime at which strict inequality hold, $\widetilde{\gamma}<\gamma_{1}$ [15]. The gap $\Delta\gamma=\gamma_{1}-\widetilde{\gamma}$ between these two exponents i.e., the difference between the expression in Eq. (9) and Eq. (10), characterizes the localization regime in the system. Combining the transition in the values of the exponent $\widetilde{\gamma}$ with the stability analysis of the steady-state solutions of Eq. (5), we identify three regimes: a regime of full equality, for $\alpha>\frac{\Delta\sigma^{2}}{\sigma^{2}}$, in which the mass is spread equally between all the non-hub nodes and the hub, a second regime for $-1<\alpha\leq\frac{\Delta\sigma^{2}}{\sigma^{2}},$ in which the mass is localized mainly on the hub, and a third regime for $\alpha\leq-1$, in which the mass is localized on the hub and a few non-hub nodes.

3 Multiple Hubs Topology

In this section, we consider the effect of $H$ hub nodes, $h_{i},$ connected to all the nodes in the system, 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, for instance, in the economic setting in the presence of more than one central bank/company. In a porous structure, it can describe different extra-cellular regions interacting with intracellular pores. It is also applied in analyzing the dynamics of control systems [4], and in machine learning algorithms. The equations of the system now read as follows:

[TABLE]

[TABLE]

Note that, as before, the total interaction rates between the hubs and the non-hub nodes, $J_{\mathrm{in}},$ and $J_{\mathrm{out}}$, are defined to be finite in the limit of infinite $N$ and $H$. The processes $g^{h}_{i}(t)$ and $g^{m}_{i}(t)$, are Gaussian processes with mean zero and variance one. Similar to the one-hub topology, we can now pass to the relative mass parameters by dividing by the average hubs mass, see D for more details. The results of the previous sections are recovered when $H=1$. Note that the equations for the relative mass are decoupled in the case $H=1$ and also in the limit of $H$ very large. In the presence of a finite small number of hub nodes, one can show that the total variance depends on the hubs value. Therefore, having finite number of hubs decreases the value of the parameter $\alpha$, and causes more equality in the system and reduces the localization. Moreover, in the simple case where all the hubs have the same statistics, such that the stochastic noise, $g^{h}_{i}(t)=g(t)$, does not depend on the hub location, $H$ hubs are equivalent to one hub with the total effective net flux $J_{\mathrm{out}}$. The limit of one non-hub node with $J_{\mathrm{out}}=J_{\mathrm{in}}=\delta$, $N=1$, and an infinite number of hubs, $H\rightarrow\infty$, is the mean field model with exponential growth of the total mass [1]. Note that when both the number of non-hub nodes and the number of hubs are very large, i.e., $N\rightarrow\infty$ and $H\rightarrow\infty$, the average mass ratio grows exponentially. The exponent depends on the average and variance difference of the stochastic multiplicative noises. The phase transition in this case is equivalent to the results in Sec. 2, i.e., there are three phases: localization on the non-hub nodes, localization on the hubs and equal spreading over all the nodes. Note that in this limit, the non-hub nodes play the same role as the hubs, since they are connected to infinitely many hub nodes. The interaction among the hubs, $\delta$ reduces localization, but doesn’t effect the growth rate. The analysis above affects mainly the sample Lyapunov exponent. The first moment Lyapunov exponent remains the same for any $H$ and $N$, since the average equations do not change. This shows that the phase transitions predicted in Sec. 2.2 are general and can be observed with small modifications to a system of multiple hubs. In addition, the transition between a logistic growth in the relative mass to an exponential growth is a function of the ratio between the $N$ and $H$.

4 A Note on MRI

In the context of MR measurements the model in Eq. (1) is a generalization of the Kärger model [17], which accounts for random changes in the diffusivity due to restricted diffusion or a non-homogeneous magnetic field. This model was already analyzed on a general network, where the importance of the spectral dimension as a measurable parameter is stated [8]. The hub topology (see Figure 1) is a simplified version of this model. We show that in this topology under the assumption that all the non-hub pores have similar properties, the average equations of the model are those for the Kärger model for two compartments [17], see B. The parameters are then as follows: $f_{0}=-q^{2}D_{\mathrm{ex}}$ and $f=-q^{2}D_{\mathrm{in}}$, $\sigma_{0}^{2}=-q^{2}\sigma_{\mathrm{ex}}^{2}$, $\sigma_{\mathrm{out}}^{2}=-q^{2}\sigma_{\mathrm{in}}^{2}$, such that the parameter $\alpha(q^{2})=2\frac{J_{\mathrm{out}}-\Delta f}{\sigma^{2}}=-2\frac{J_{\mathrm{out}}+q^{2}\left(D_{\mathrm{in}}-D_{\mathrm{ex}}\right)}{q^{2}\left(\sigma_{\mathrm{ex}}^{2}+\sigma_{\mathrm{in}}^{2}\right)}$, is controlled by the gradient of the applied magnetic field, incorporated in the value of $q$. For example consider the basic Stejskal-Tanner sequence [23], which is composed of two gradients pulses of the magnetic field with magnitude $G$ in opposite direction and with duration $\delta$. The pulses are separated by a diffusion time $\Delta$. For $q$ one takes the wave vector, $q=\frac{\delta\gamma G}{2\pi}$, where the parameter, $\gamma$, is the gyro-magnetic ratio. The average equations for the magnetization in the $(x,y)-$plane, and with a magnetic field gradient in the $\hat{z}$ direction, read

[TABLE]

[TABLE]

Note that, the wave vector $q$ turns on the stochastic dynamics. We denote $\delta D=D_{\mathrm{ex}}-D_{\mathrm{in}}$, and $\delta\sigma^{2}=\sigma_{\mathrm{ex}}^{2}-\sigma_{\mathrm{in}}^{2}$. The multiplicative white noise accounts for the random diffusivity changes of the medium due to restricted diffusion. Based on the analysis in Sec. 2.2, one can find the transition point in terms of the wave vector $q$. The transition point is at $q_{\mathrm{c}}=\sqrt{\frac{J_{0}}{\delta D-\frac{\delta\sigma^{2}}{2}}}$. The average signal reveals only the transition at $q_{\mathrm{c}}$. Note that, the signal decay is also affected by the noise variance of the extracellular space. Taking typical values such as, $J_{\mathrm{out}}=\frac{1}{1800\mathrm{msec}}=0.5556[\frac{1}{\mathrm{sec}}]$, $D_{\mathrm{ex}}=2e-5[\frac{\mathrm{cm}^{2}}{\mathrm{sec}}]$, $D_{\mathrm{in}}=0.1e-5[\frac{cm^{2}}{\mathrm{sec}}]$. Then the critical value is $q_{\mathrm{c}}\approx\sqrt{\frac{J_{\mathrm{out}}}{D_{\mathrm{ex}}-D_{\mathrm{in}}+\frac{\sigma_{\mathrm{in}}^{2}}{2}}}=\sqrt{\frac{0.5556}{1.9e-5}}\approx 171[\frac{1}{\mathrm{cm}}]=0.1[\frac{1}{\mathrm{\mu m}}]$. A larger variance in the non-hub pores will set the transition for a lower value of $q$, value. Whereas larger variance in the diffusion in the extracellular space will set the transition to a higher value of $q$ value. The decay of the signal has a bi-exponential form as predicted by the Kärger model:

[TABLE]

The stochastic model is a natural generalization of the Kärger model. Using this generalization, we are able to explore and analyze the behavior of the model in the presence of complex topological structures, as well as the effect of changes in the apparent diffusivity as a result of stochastic noise. Note that this transition appears also in the presence of any interaction $\delta$ among the non-hub nodes. Figure 3 presents the behavior of the eigenvalues as a function of $N$.

5 Discussion and Conclusion

We have presented a stochastic model that describes diffusion on a graph with an additional multiplicative stochastic noise. We analyze this model on a directed graph with one hub node that is connected to a large number of non-hub nodes. We derive a non-linear equation for the average mass ratio between the non-hub nodes and the hub. This equation describes logistic growth. It has two phases one in which the overall mass is mainly concentrated on the hub, and the other in which there is a “detailed balance” such that the steady-state depends on the ratio between the exchange rate between the hub and the non-hub nodes. We show that this model is completely solvable in the large $N$ limit. In addition, we identify the phase transitions of the model in terms of the two parameters $\alpha$ and $\frac{\Delta\sigma^{2}}{\sigma^{2}}$. We show that in order for localization phase occur 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 the amount of resource given by the hub, provided that it is finite and non-zero. We generalize this analysis to a system of multiple hubs. We show that in the limit of infinitely many hubs the growth of the system becomes exponential.

The model has numerous applications. We introduce an application of this model in the context of MR measurements of complex structures. Our results in this context may provide a theoretical framework that may help interpret and propose new MR experiments to identify the concentration phases that we see theoretically. This may have impact on the prediction of the underlying measured geometry. Our results and analysis can also be of interest in other applications, for example, in predicting economic growth, and in analyzing the stability of control systems.

Acknowledgment

We would like to thank Prof. Ofer Pasternak for proposing the idea for the paper.

Appendix A Transformations of Variable

In this section, we derive the relative magnetization equations, Eq. (4). We use Itô’s formula in order to perform the change of variables

[TABLE]

where $A$ and $B$ are the coefficients of the stochastic equations Eq. (2) and Eq. (3), respectively, defined as follows: $A_{i}=J_{i0}h_{0}-J_{0i}m_{i}-\delta\,\underset{j}{\sum}L_{ij}m_{j}+\frac{\sigma_{i}^{2}}{2}m_{i}+f_{i}m_{i}$, $A_{0}=\underset{j}{\sum}\,J_{0j}m_{j}-\underset{j}{\sum}\,J_{j0}h_{0}+\frac{\sigma_{0}^{2}}{2}h_{0}+f_{0}h_{0}=\frac{J_{\mathrm{in}}}{N}\underset{j}{\sum}\,m_{j}-J_{\mathrm{out}}h_{0}+\frac{\sigma_{0}^{2}}{2}h_{0}+f_{0}h_{0}$, and $B_{ij}=\delta_{ij}\sigma_{\mathrm{out}}m_{i}$, and $B_{0i}=B_{i0}=\delta_{i0}\sigma_{0}m_{0}$.

[TABLE]

We introduce the dimensionless parameters $\alpha=2\frac{J_{\mathrm{out}}-\Delta f}{\sigma^{2}}$, and $\beta=\frac{2J_{\mathrm{in}}}{\sigma^{2}}$. The equations are written under the assumption that the interaction among the nodes and the hub is described by $J_{i0}=\frac{J_{\mathrm{out}}}{N}$, and $J_{0j}=\frac{J_{\mathrm{in}}}{N}$, respectively. We also assume, for simplicity, that all the non-hub compartments obey the statistics $\sigma_{i}=\sigma_{\mathrm{out}}$ and $f_{i}=f$. The last transition is under the assumption of mean-field interaction $\delta$ among the non-hub nodes. In the Itô form in the limit of $N\rightarrow\infty$, we have

[TABLE]

Note that one cannot take the limit $N\rightarrow\infty$ in Eq. (3), since the variables $m_{i}$ depend on the stochastic noise of the hub, $g_{0}$. Taking the sum and letting $N\rightarrow\infty$, we arrive to the following deterministic equation describing the growth of the average relative mass of the $N$ non-hub nodes as a function of time:

[TABLE]

The steady-state solutions of this non-linear equation are: $\overline{M}_{1\mathrm{eq}}{\rightarrow}-\frac{\Delta f-J_{\mathrm{out}}-\frac{\sigma^{2}}{2}}{J_{\mathrm{in}}}=\frac{1+\alpha}{\beta}$ as ${N\rightarrow\infty}$, and $\overline{M}_{2\mathrm{eq}}{\rightarrow}0$ as ${N\rightarrow\infty}$. The equation is also valid when $\delta\neq 0$.

A.1 Finite-$N$ corrections

In this subsection, we consider finite $N$ corrections to the average equation. With this effect accumulated for the equation reads

[TABLE]

where $\varepsilon=\frac{1}{N}$. We denote $a=\frac{\sigma^{2}}{2}J_{\mathrm{out}}$, $b=-\frac{\sigma^{2}}{2}\beta$, and $c=\frac{\sigma^{2}}{2}(\alpha+1)$. The equation has a Riccati form. Taking the first-order correction in $\varepsilon$, $\overline{M}(t)=\overline{M}_{0}(t)+\varepsilon\overline{M}_{1}(t)$, we get

[TABLE]

The solution for $\overline{M}_{0}(t)$ is the logistic function as before. Next,

[TABLE]

The solution for $\overline{M}_{1}(t)$ is given by

[TABLE]

Substituting the expression of the solution to $\overline{M}_{0}(t)$, one can show that $\underset{t\rightarrow\infty}{\mathrm{lim}}\,\overline{M}_{1}(t)$ is finite, meaning that the correction of order $\frac{1}{N}$ to Eq. (16) is negligible in the large $N$ limit.

Appendix B Moments Lyapunov Exponent

In this section, we derive the moment Lyapunov exponent Eq. (9). For this purpose, we present the first moment equations of Eq. (2) and Eq. (3) for a finite number of non-hub nodes $N$ and $\delta=0$. These equations can be derived using the Fokker-Plank equation or alternatively the Feynman-Kac formula [8, 24, 16, 15]. They read

[TABLE]

[TABLE]

Here it is assumed that all the non-hub nodes are independent and with the same dynamics and denoting the total mass, $\langle\overline{m}\rangle=N\langle m_{i}\rangle$. Summing over $N$ in Eq. (21), we get

[TABLE]

Eqs. (22) and (20) can be written in vector form as

[TABLE]

This is a simple system of linear equations, and the dynamics it describes is determined by the eigenvalues of the matrix $A$. The resulted eigenvalues are than,

[TABLE]

and,

[TABLE]

The solution is then given by \left(\begin{array}[]{c}\langle m_{0}\rangle(t)\\ \langle\overline{m}\rangle(t)\end{array}\right)=Av_{1}e^{\lambda_{1}t}+Bv_{2}e^{\lambda_{2}t}=A\left(\begin{array}[]{c}\lambda_{1}-d\\ c\end{array}\right)e^{\lambda_{1}t}+B\left(\begin{array}[]{c}\lambda_{2}-d\\ c\end{array}\right)e^{\lambda_{2}t}=A\left(\begin{array}[]{c}\Delta f-J_{\mathrm{out}}+\frac{\Delta\sigma^{2}}{2}\\ J_{\mathrm{out}}\end{array}\right)e^{\lambda_{1}t}+B\left(\begin{array}[]{c}0\\ J_{\mathrm{out}}\end{array}\right)e^{\lambda_{2}t}. Plugging this expression into the equation of the first moment Lyapunov exponent we get

[TABLE]

Appendix C Sample Lyapunov Exponent

Here, we calculate the sample Lyapunov exponent of the mass of the system, defined as

[TABLE]

To prove the formula (10) in the main text, we consider at the system of equations

[TABLE]

[TABLE]

Using a transformation of variable given by the Itô formula

[TABLE]

Integrating with respect to time and dividing by $t$, we have

[TABLE]

Letting the limit of $t,N\rightarrow\infty$, we get

[TABLE]

Here we used the deterministic law of $\overline{M}$ given that higher corrections in $N$ are negligible, see A.1 for more details. Using the stationarity property of the Gaussian processes $g_{0}(t)$ and $g_{i}(t)$, and the ergodic theorem:

[TABLE]

This formula is obtained under the assumption that the fluctuations of the noise are finite (Novikov condition), i.e.,

[TABLE]

and,

[TABLE]

In order to calculate the above integrals one needs to specify how the number of nodes in the graphs growth with time. We show that if the number of nodes grows exponentially in $t$, i.e., $N\sim e^{\varepsilon t}$, then there exists a localization phase. In addition, the fluctuation, conditions 27 and 28, are satisfied. These conditions are verified in C.1 below. Since the fluctuations are finite, we are left with calculating integrals of the logistic function $\overline{M}(t)$. We use the following properties of the logistic function:

[TABLE]

and

[TABLE]

In our case, $A=\frac{\beta}{\alpha+1}$, $B=\frac{1}{\overline{M}_{0}}-\frac{\beta}{\alpha+1}$, and $\xi=\frac{\sigma^{2}\left(\alpha+1\right)}{2}$. Using Eqs. (29), and (30) it is easy to show that

[TABLE]

where $T$ is some finite time. In the same manner, one can calculate the second term in Eq. (26):

[TABLE]

provided that $\mathrm{max}\{-\xi,0\}\leq\varepsilon$. Note that, in order to have finite fluctuations, $\sigma^{2}\leq\epsilon$, see C.1. Substituting the results in Eq. (31) and Eq. (32) we obtain

[TABLE]

Therefore, for any exponential growth rate $\sigma^{2}\mathrm{max}\{-\frac{\alpha+1}{2},1\}\leq\varepsilon$, the fluctuations and sample Lyapunov are finite.

C.1 Noise fluctuations

In this subsection, we prove that the fluctuation are finite when $N\sim e^{\varepsilon t}$, i.e., the conditions in (27) and (28) are satisfied. That condition (28) is satisfied, is readily seen since the function under the integral is bounded between zero and one, and so the integral itself is also finite:

[TABLE]

In order to calculate the fluctuation of the non-hub nodes, i.e., establish (27), we approximate the limit, $\underset{N\rightarrow\infty}{\mathrm{lim}}\frac{1}{N}\sum_{i}M_{i}^{2}\rightarrow\langle M_{i}^{2}\rangle$; since $M_{i}$ are i.i.d. random variables. The second moment can be calculated using the Fokker-Plank equation Eq. (7):

[TABLE]

Averaging over the second moment yields:

[TABLE]

The solution for $\tilde{\delta}=0$ is

[TABLE]

The fluctuations of the non-hub nodes, i.e., the last term in Eq. (25) are then, given by

[TABLE]

Substituting here the logistic function $\overline{M}(t)$ (Eq. (29)), we get

[TABLE]

Therefore, the fluctuation are finite for any $\varepsilon\geq\sigma^{2}$.

Appendix D Multiple Hubs Derivation

In this section, we derive the normalized equation in the multiple hubs models of Sec. 3. Taking normalized variables $M_{i}=\frac{m_{i}}{\overline{h}}$, $x_{i}=\frac{h_{i}}{\overline{h}}$, so that, $\overline{M}=\frac{1}{N}\sum_{i}\frac{m_{i}}{\overline{h}}$, and $\frac{1}{H}\sum_{j}x_{j}=1$, Eq. (11), and Eq. (12) are transformed into the following set of equations:

[TABLE]

[TABLE]

Set

[TABLE]

and

[TABLE]

so that $\frac{1}{H}\sum_{i}^{H}\sigma_{i}\left(x\right)\phi_{i}x_{i}=0$, and $\nu_{i}^{2}\left(x\right)-\nu_{i}^{2}=\sigma_{i}^{2}\left(x\right)-\sigma_{i}^{2}$, and $\Delta f=f-q$. In the limit of $H,N\rightarrow\infty$, the variances are constants, $\nu^{2}\left(x\right)=\nu^{2}+\sigma^{2}$, and $\sigma^{2}\left(x\right)=2\sigma^{2}$. In this limit,one can average Eq. (34), since the variables are decoupled.