Fractional calculus ties the microscopic and macroscopic scales of complex network dynamics

A two-state master equation based decision making model has been shown to generate phase transitions, to be topologically complex and to manifest temporal complexity through an inverse power-law probability distribution function in the switching times between the two critical states of consensus. These properties are entailed by the fundamental assumption that the network elements in the decision making model imperfectly imitate one another. The process of subordination establishes that a single network element can be described by a fractional master equation whose analytic solution yields the observed inverse power-law probability distribution obtained by numerical integration of the two-state master equation to a high degree of accuracy.


I. INTRODUCTION
The extraordinary advancements made in the physical sciences at the beginning of twentieth century originated from the effort to understand the simplest, that is, the most fundamental, elements of physical phenomena. The separation of matter and light into its basic constituents: electrons, protons, atoms, and molecules in the first case and photons in the second, enabled scientists to explain the puzzling properties of solid materials, such as their sound transmission and electrical properties, as well as, light emission and reflection characteristics. The discipline of statistical mechanics demonstrated that the behavior of systems composed of millions of individual particles can be captured with simple laws, involving only their average properties [44]. At that time it seemed as if the Aristotelian approach of reductionism, where complex phenomena consisted of nothing more than the sum of essential elements, was the prescription underlying the correct method of scientific discovery.
The French mathematician Poincaré was probably the first to rigorously demonstrate the failure of this attractive but overly optimistic method by extending Newton's law of universal gravitation to a system consisting of three celestial bodies [47]. His mathematical treatment of the three-body problem demonstrated that unlike the two-body problem of the earth and sun, a planet's orbit in the three-body system need not be periodic. Poincaré proved that the long held belief that planetary motion could be built up from the superposition of simple cycles was false. This unanticipated behavior emerges from nonlinear dynamics and has come to be known as chaos theory.
Another perplexing phenomena that violates the principle of superposition is critical phase transitions in magnetic materials. If criticality were truly nothing more than collective behavior resulting from the superposition of some basic building blocks, why then does the change in an external parameter, for example temperature, induce such a dramatic shift in its macroscopic behavior?
More recently the identification of emergent phenomena across multiple disciplines, from the swarming of insects [68], the schooling of fish [30] and the flocking of birds [10] observed in animal groups by naturalists; to the spatiotemporal activity of the brain [6,11,19] observed by neurophysiologists; to the collective and cooperative behavior observed in social groups studied by psychologists and sociologists; all demonstrate collective behavior reminiscent of particle dynamics near the critical phase transitions studied by physicists [54]. Each of these disciplines has demonstrated the need to investigate the dynamics of complex networks across scales in order to develop a deeper understanding of how largescale behavior emerges from microscale dynamics and the sensitivity of the observed behavior to those dynamics.
Of particular interest to us here are the biological fields in which we observe a need for a system wide approach [49]. The recent discoveries in biology were propelled by the successes of molecular biology and genetics that have made available genomic blueprints of numerous organisms, which are complemented by extensive experimental data describing cell functions. At the same time however the realization came that biological function emerges out of the interaction of numerous molecular components, such as depicted in Figure 1 for the human brain, making the detailed knowledge of specific components at any level of organization insufficient to capture macroscopic functionality. There is probably no better example of this limitation then the study of the neurological systems, whose goal is to understand, predict and ultimately modify (in order to heal) brain function. Ongoing initiatives of the Human Connectome Project [26], the Human Brain Project [27] or the Allen Brain Atlas [28] illustrate the fact that the system wide approach and integration of data from across different spatial and temporal scales has become the norm in modern scientific disciplines.
Despite experimental developments, the ability of science to make theoretical predictions of the behavior of complex networks is still in its infancy. The adoption of methods from non-equilibrium statistical physics have demonstrated limitations, resulting from the fact that liv- ing systems, in contrast to inert physical materials, are extremely heterogeneous, non-generic, highly specialized and operate far from an equilibrium state [18]. Herein we demonstrate that what was for a very long time a niche branch of mathematics, the fractional calculus, might very well be able to span the gap between the inert materials of physics and the living networks of biology.
Although developed along side the classical calculus, fractional differential equations have only recently been shown to be a convenient way to describe the dynamics of complex phenomena characterized by long-term memory and spatial heterogeneity [36,46,64]. Fractional differential equations were demonstrated to capture the time evolution of fractal processes, such as in anomalous diffusion, viscoelasticity and turbulent fluid flow, as reviewed by West and Grigolini [65]. In spite of the success of the mathematical descriptions of such processes there has been a lack of identification and interpretation of mechanisms that entail fractional dynamic equations in the context of complex networks. Herein we provide an explanation for one source of a fractional differential equation that describes the dynamics of a complex network using a fractional master equation.
Here the utility of the fractional calculus is demonstrated by capturing the dynamics of the individual elements of a complex network from the information quantifying that network's global behavior. The phase transitions observed in complex social and physiologic networks suggest the wisdom of using a generic model from the Ising universality class to characterize network dy-namics. Using such a model it is then possible to demonstrate that the individual trajectory response to the collective motion of the network is described by a linear fractional differential equation. The solution to this fractional equation is obtained through a subordination procedure without the necessity of linearizing the underlying dynamics, that is, the solution retains the influence of the nonlinear network dynamics on the individual. Moreover the solutions to the fractional equations of motion suggest a new direction for designing control mechanisms for complex networks.
In Section II we consider the dynamics of a complex network described by a two-state master equation. The decision making model (DMM), defined by the two-state master equation, undergoes phase transitions at a critical value of the control parameter [58]. It is understood that at criticality the short-range local interactions between the two-state elements generate long-range global correlations, thereby producing effective long-range interactions. Consequently at criticality there is global cooperation among the network elements.
An individual disconnected from the network is assumed to choose randomly between two states with an exponential distribution of decision times and a given average decision time. When coupled to the other individuals of the network, the global distribution for the time intervals between decisions is determined to be inverse power-law [66]. In this latter case the power-law index of the survival probability is a measure of the complexity of the underlying dynamics determining whether that process is non-stationary and non-ergodic [57,58]. In Section III the DMM network dynamics is incorporated into that of an individual element through a process known as subordination. In order to formalize the subordination process we introduce the concept of subjective time to distinguish between clock time that determines the activities of the network and operational or subjective time that determines the activities of the individual.
The subordination process results in the two-state master equation of the DMM being replaced by a fractional master equation for the individual whose solution is shown to be a Mittag-Leffler function in Section III. This predicted behavior of the single element dynamics is compared with the numerical results from the DMM implemented on a two-dimensional lattice and found to be in excellent agreement. In Section IV we draw some conclusions.

II. DECISION MAKING MODEL (DMM)
The DMM realized on a complex networks represents the dynamics of the probability for an individual to be in either of the two states: yes or no, up or down, on or off. The model is based on the cooperative interaction of N elements, each of which is described by the two-state master equation [57,58] d dt p The quantity p (i) j (t) is the probability of the element i=1,...,N in the network being in the state j = 1, 2 at time t and the probability is normalized at all times such that The network dynamics are determined by the choice of the functional form of the transition rates in the twostate master equation (Eqs. (1) and (2)). Each probability p The symbol N (i) denotes the total number of elements linked to the i -th element and N implies that the transition rates become erratic functions of time.

A. When every element is interconnected
In the ATA coupling case when the total number of elements within the network becomes infinite (N −→ ∞) the fluctuation frequencies collapse into probabilities according to the law of large numbers. In physics this replacement goes by the name of the mean field approximation, in which case the transition rates in the master equation (1) and (2) are written as The formal manipulation of the master equation even in this simplified case in made a little easier if we introduce a new variable defined as the difference in the probabilities Subtracting Eq.(2) from Eq.(1) after some algebra yields the highly nonlinear rate equation for the difference variable where the nonlinearity enters through the transition rate dependence on the difference variable in the mean field approximation. By inserting Eqs.(10) and (11) into Eq.(9) we obtain and the network dynamics are determined by the potential function V (Π), which is a symmetric double well potential with the explicit form Note that the network is not described by a Hamiltonian and yet the global dynamics in indeed described by an effective Hamiltonian, that being the double well potential given by Eq. (13) and depicted in Figure 2.
The cooperative behavior of the infinitely large ATA coupled network described by Eq.(12) is that of an overdamped particle hopping from one potential minimum to the other, whose position is Π within the potential Eq.(13). For K < 1, half of the nodes are in one state and half are in the other because there is only a single broad minimum in the potential. At the critical value of the control parameter K = K C = 1.0 a bifurcation occurs and the potential develops two wells separated by a barrier as discussed by Turalska et al. [57]. The height of the barrier increases with the value of the control parameter. It is now convenient to define the stochastic global variable where s i (t) is the state of element i at time t. The variability of the global variable characteristic of the entire network, capturing the cooperation between units at any moment of time. It is interesting that at the critical value of the control parameter the ATA version of the DMM undergoes a phase transition. Note that the amplitude of ξ(t) depends on the value of the control parameter K. When K = 0, all elements in the network are independent Poisson processes; thereby an average taken at any moment of time over all of them yields zero. Once the value of the coupling becomes nonzero, K > 0, single elements are less and less independent, resulting in nonzero averages. The quantity K c is the critical value of the control parameter K, at which point a phase transition to a global majority state occurs.
In numerical calculations we use the time average ξ eq = |ξ (t)| as a measure of this global majority. More precisely, after an initial 10 6 time steps, the average is taken over the same number of the consecutive time steps of the model. In Figure 3 the thin line indicates the ATA phase transition as measured by ξ eq . The other phase transitions indicated are for the Ising model (dashed line) and the DMM model on a two-dimensional lattice as discussed in Section II B and elsewhere [58].
Real network are not ATA coupled since interactions typically have finite range and elements are spatially separated. Thus, the ATA approximation may be valid for small networks but certainly breaks down for large systems. Moreover, real-world networks have finite numbers of elements. It is therefore useful to examine how strongly the mean field solutions are violated when we relax these constraints. The stability condition can be violated in at least two different ways. The first way is by reducing the number of elements N to a finite value. The second way is by restricting the number of links so the network no longer has ATA coupling. In real-world networks both sources of equilibrium disruption are expected to occur. For the time being we retain the ATA coupling within the networks and consider the number of elements N to be finite. In this latter case we can no longer make the mean field approximation and the dynamic picture stemming from the above master equation is radically changed.
If the number of elements is still very large, but finite, we consider the mean-field approximation to be nearly valid and replace the average Eq. (14) with the stochastic quantity where f (t) is a small amplitude random fluctuation. After inserting Eq.(15) into (12) and some straight forward algebra we obtain the stochastic differential equation [8,24] to lowest-order in the strength of the fluctuations: The additive fluctuations ε (t) have amplitudes that are computationally determined to be on the order of 1/ √ N . Note that the double-well potential in the mean field approximation persists in the present description even though we have relaxed the mean field approximation to a finite number of network elements. The random fluctuations resulting from the finite size of the network induces transitions between the two states of the potential well. Consequently, for a network with a finite but large number of elements the phase synchronization of Eq. (12) is not stable and the stochastic Langevin equation Eq. (16) determines the dynamics of the network. Furthermore, the fluctuations can drive the particle from one well of the potential to the other when its amplitude is sufficient to traverse the barrier separating the wells. However, here the fluctuations arise from the finite number of elements in the network rather than from non-existent thermal excitations and are consequently suppressed as the network size increases.
Although Eq. (16) is written in the continuous time representation, in practice the numerical calculations of DMM correspond to the adoption of a finite integration time step ∆t = 1. Note that the stochastic rate equation Eq.(16) replaces Eq. (12) in the case of a finite N , and that Eq. (12) is recovered in the ideal case N = ∞. We incorporate the ATA coupling condition with a finite number of elements by numerically integrating the the master equation for each element in the network and then calculating the number of elements in each of the two states. In Figure 4 the fluctuating global variable ξ(t) is depicted as a function of time, under differing conditions. Notice that with increasing N the fluctuation ξ(t) become more distinctly dichotomous-like, with an increasingly sharp transition from the 'up' to the 'down' state. This pattern corresponds to the entire network keeping a decision for a longer and longer time as the size of the network increases. The condition of a decision lasting forever is reached in the ideal case N = ∞. The global variable fluctuates between the two minima of the double-well potential as described by Eq.(16) for K = 1.05 > K C . The single element follows the fluctuations of the global variable, switching back and forth from the condition where the upper state is preferred statistically to that where the lower sate is preferred statistically. The complete properties of the DMM on an ATA network are explored by Turalska et al. [57,58].

B. Nearest neighbor coupling
In this section we consider the topology of a simple two-dimensional lattice and confine the coupling between elements to its four nearest neighbors thereby setting N = 4 in the transition rates of the two-state master equation. Similarly to the ATA case, the fluctuations of the global variable ξ(t) show pronounced transition as a function of the coupling parameter K. As seen in Figure 5b, the global variable shifts from a configuration dominated by randomness to an organized state once the control parameter is increased above the critical value K C . For values of the control parameter K corresponding to the disorganized phase K < K C , single elements of the lattice are only weakly influenced by the decisions of the neighbours. Thus, the fluctuations of the global order parameter ξ(t) are characterized by small amplitude and very fast oscillations about the zero-axis. For K > K C , the interaction between individuals give rise to a majority or a consensus state, during which a significant number of agents adopts the same opinion at the same time.
At the same time the global behavior is undergoing a phase transition, the presence of the lattice apparently exerts only very subtle influence over the behavior of single individuals. The latter influence can be observed as a change in the interval timing for a single element as the control parameter is increased (Fig. 5a). Note that if attention is concentrated on a single network element when the network is in a consensus state that individual would still appear to make transitions according to an exponential distribution as exhibited in Figure 5c. The strict exponential is indicated by the black dotted curve. The single particle survival probabilities do not look too much different, the light gray dashed curve with the subcritical value K = 1.5 < K c is very close to the exponential. The remaining single particle curves, whether critical K c ≈ 1.70 or supercritical K > K c appear to be exponential on this graph.
To characterize the changes in temporal properties of the microscopic and macroscopic variables we evaluate the survival probability function Ψ(τ ) of time intervals τ between consecutive events defined as changes of the state or crossing of the zero-axis, for the single element and the global variable, respectively. These calculations unveil modest deviation of the survival probability for a single individual from the exponential form Ψ(τ ) = exp[−g 0 τ ]. The strict exponential is indicated in Figure  5c by the black dashed curve. The single particle survival probabilities do not look too much different, the blue curve with the subcritical value K = 1.5 < K C is very close to the exponential. The remaining single particle curves, whether critical K = K C ≈ 1.70 or supercritical K > K C appear to be very nearly exponential on this graph. The difference in the behavior of the individual from that in the non-interacting state would be that she tends to be more reluctant to change her mind.
The deviation of the individual survival probability from the exponential form in Figure 5c appears to be modest when compared with the dramatically greater deviation of the survival probability of the global variable from the exponential as depicted in Figure 5d. The average network behavior differs markedly as the control parameter is increased from the subcritical through the supercritical regions. However the influence of the global variable on the behavior of the individual does not appear to induce a significant change. For the individual the change is however a subtle yet profound difference and is a direct result of the imitation mechanism, that is the ERH. So if the individual survival probability is not exponential, what is it? To answer this question we turn our attention to describing an alternate construction of the dynamics of the individual elements.

III. SUBORDINATION AND FRACTIONAL DYNAMICS -TWO VIEWS OF TIME
In this section we demonstrate the equivalence between a fractional trajectory that is the solution of a Caputo fractional differential equation, and the ensemble average trajectory that results from a subordination process. We here consider only fractional derivatives of the Caputo type, in part because it requires the least amount of explanation. However for the aficionado we note that an approach using Riemann-Liouville fractional derivatives would be equivalent as long as the initial conditions are properly specified. We begin the discussion with a derivation of the fractional derivative from a subordination argument.
The master equation for a single isolated individual is, with the index suppressed, according to the numerical simulation given by (17) whose discrete solution is Here φ(n) is the difference variable for a single individual chosen from the network at random and as n → ∞ and ∆τ → 0 such that clock time is t = n∆τ we have the apparently trivial result Subordination implies the existence of two different notions of time [50,55]. One is the operational time τ , which is the internal time of a single individual, with an element generating the ordinary dynamics of a nonfractional system. The other notion is experimental time t; the time as measured by the clock of an external observer. Typically in the operational time frame the temporal behavior of an element is regular and evolves exactly according to the ticks of a clock leading to the exponential. Therefore it is assumed that the trajectory of a network's element in operational time is well defined and given by φ(τ ), which is the solution given by Eq. (19).
It is perhaps worthwhile to point out that this notion of two different times was introduced into psychology in the middle nineteenth century and lead to the general Weber-Fechner law. It has been further developed in a contemporary setting to explain the observation of 1/f noise in cognition by discriminating between subjective and objective times, that being operational and chronological time, respectively.
In operational time an element's behavior appears ordinary, but to an experimenter observing the elements from outside the network their temporal behavior appears erratic, evolving in time then abruptly freezing in different states for extended time periods. Because of the random nature of the experimental or chronological time evolution of the elements the subordination process involves an ensemble average over many realizations of the element's dynamics each evolving according to its own internal clock, independent of one another. Making an ensemble average over a large number of realizations results in a smooth average trajectory, which is equivalent to the fractional trajectory.
To find the average behavior we move from the operational time solution to the experimental time solution adopting the subordination interpretation. We assume that the chronological time lies in the interval (n − 1)∆τ ≤ t ≤ n∆τ and obtain It is evident that the trajectory resulting from the subordination process inherently involves an ensemble average. Here we see that Eq. 20 replaces the solution to the single element two-state master equation of the DMM.
Note that ψ n (t) dt is the probability that n events have occurred, the last one in the time interval (t, t + dt). The function Ψ(t) denotes the probability that no event occurs up to time t and is given empirically by numerical calculation in Figure 5d and mathematically by Eq. (27). The occurrence of an event corresponds to activating a decision with (1−g 0 ∆τ ), so that activating n such events transforms the initial condition φ(0) into the product (1 − g 0 ∆τ ) n φ (0). This form of the equation is kept from time t , at which time the last event occurs, up to time t, the time interval t − t being characterized by no event occurring. Of course, the expression Eq.(20) takes into account that the number of possible events may range from the no-event case to that of infinitely many events. The conditions necessary for this result to occur are discussed by Svenkenson et al [55]. To interpret the physical meaning of Eq. (20), consider each tick of the internal clock n of an element measured in experimental time as an event. Since the observation is made in experimental time, the time intervals between events are random. We assume that the waiting times between consecutive events are identically distributed independent random variables. The integral in Eq. (20) is then built up according to renewal theory. After the n-th event the individual changes from state φ(n) to φ(n + 1), where it remains until the action of the next event. The sum over n takes into account the possibility that any number of events could have occurred prior to an observation at experimental time t. The events occur randomly with a waiting-time probability density function (pdf ) ψ(t) and survival probability Ψ(t). The waiting-time pdf is related to the survival probability through Taking advantage of the renewal nature of the events, the waiting-time pdf for the n-th event in a sequence is connected to the previous event by ψ n (t) = t 0 ψ n−1 (t )ψ(t − t )dt (22) At this point it is useful to introduce Laplace variables in our discussion. The Laplace transform of a function f (t) is denoted To find an analytical expression for the behavior in experimental time it is convenient to study the Laplace transform of Eq. (20) where we assume the intervals between successive transitions are independent of one another, it is a renewal process. Consequently, the waiting time pdf for n transitions is the product of n single transition pdf 's: which was used to collapse the convolution of Eq. (22) to the power-law form in Eq. (24).
Consequently the time t is derived from a waiting-time pdf given by that of the network as a whole: and the survival probability is empirically determined from Figure 5d to be The Laplace transform of the survival probability in terms of that for the waiting-time pdf is Inserting these last two expressions into Eq.(24) and evaluating the sum yields whose inverse Laplace transform yields: a generalized master equation.

A. Fractional Langevin Equation
The function Φ (t) in the Eq.(30) is a memory kernel containing the information on how the other elements in the network influence the dynamics of the individual element under study. The Laplace transform of the memory kernel is and a complete discussion of its properties is now given in textbooks [65]. Equation (31) so that Eq.(29) reduces to and the rate parameter has the value We now assume that the exact equation for the individual dynamics has both an average and a fluctuating part just as in the mean field treatment of the double well potential. Consequently in terms of the Laplace variables we have the stochastic equation (35) which has the inverse Laplace transform [64] Equation ( and is completely equivalent to that determined using the Riemann-Liouville form of the fractional derivative.
The solution to the fractional Langevin equation is given by where the homogeneous solution to the fractional equation is the Mittag-Leffler function ; θ > 0. (40) and the kernel of the integral is in terms of the Mittag-Leffler function of the second kind .
The dynamics of the individual is determined by the exact Laplace transform equation Eq.(28). However it is notoriously difficult to obtain analytic expressions by direct inversion of the resulting equations due to the complexity of the exact form of the Laplace transform memory kernel. Consequently, the strategy is to consider the asymptotic forms of the solution, which was done by examining the solutions to the fractional Langevin equation given by Eq. (39). We find that the properties of the fluctuations change as the control parameter is varied from the subcritical, critical and supercritical regions.

B. Solution domains
In all three regions of DMM dynamics, subcritical, critical and supercritical, the single elements used in the evaluation of the probability difference ϕ(t) were selected at random among all nodes of the lattice. The calculations were done on a 100 × 100 node two-dimensional lattice, with nearest neighbor interactions. The time-dependent average solution calculated over an ensemble of randomly chosen individuals is depicted in Figure 6, where the average is taken over 10 4 independent realizations of the dynamics. The analytic solution is obtained by averaging Eq.(39) over an ensemble of realizations of the single particle trajectory to obtain the Mittag-Leffler function: From the series form of the Mittag-Leffler function it is evident that for µ = 2 the average probability difference would be an exponential. Consequently the influence of Note that the early time behavior of the Mittag-Leffler function is the indicated stretched exponential. On the Figure 6a, the region where the black dashes of the Mittag-Leffler function fit diverge from the data is the onset of the inverse power-law tail of the Mittag-Leffler function. An exponential truncation of the Mittag-Leffler function would fit the data throughout its domain. The rational for a truncated Mittag-Leffler function will be taken up elsewhere. The fitting of the analytic solution at early times to the DMM numerically generated curves is certainly very good in the subcritical domain with R 2 = 0.9968.
As the critical point is approached from the subcritical region the random influence of fluctuations are diminished as would be expected due to the formation of clusters as the network undergoes a phase transition and encounters critical slowing down. The plunging stretched exponential that was observed in the subcritical region as seen in Figure 8a is replaced with a more gently decreasing function. The time-dependent average solution of a randomly chosen individual in the critical regime is depicted in Figure 8b. It is evident by comparing these data with the curve in Figure 8a that the average solution does not decrease as quickly in time and there is less variability asymptotically in time. This behavior is reflected in the value of the power-law index which is determined to be µ = 1.808 with a quality of fit given by R 2 = 0.9989. Note how well separated the solution is from the exponential function given by the light grey dashed curve. But here again an exponential truncation of the Mittag-Leffler function might provide a better overall fit to the data.
In the supercritical region of the control parameter it is evident from the fit of the analytic solution to the data depicted in Figure 8c that the solution extends far beyond that found in either the subcritical or critical domains with the Mittag-Leffler function solution extending far into the inverse power-law region. Here the power-law index is fitted with the value µ = 1.534 with R 2 = 0.989. The measured inverse power-law index is very close to that obtained for the global survival probability obtained from the numerical calculation of the DMM lattice network.

IV. CONCLUSION
In summary the last few years have witnessed substantial attention focused on the role of criticality [41] to explain the function of complex networks, from flocks of birds [10], to neural networks [20] to the brain [11]. At criticality, the short-range links of Ising-like cooperative models are converted into long-range interactions turning a set of N distinct elements into an organized network behaving as a single individual with extended cognition [7,60]. A complex network at criticality generates 1/f noise [32], which is thought to be of relevance for cognition [62], with the interesting property of maximizing information transport [3,65]. Moreover the network dynamics has a subtle but profound influence on the behavior of each individual within the network.
The numerical solution of the DMM on a 100×100 lattice with elements at each of the nodes and with nearest neighbor interactions gives rise to an inverse power-law survival probability for the global variable introduced in Section II. Using the theory of subordination, that being the time experienced by an individual, with the influence of the network entering into the individual's dynamics through the distribution of critical events, the dynamics of an individual is determined by a fractional Langevin equation.
The explicit form of the fractional Langevin equation depends on whether the network dynamics is in the subcritical, critical or supercritical domains. In all three domains the Mittag-Leffler function solution to the fractional Langevin equation is sufficient to describe the dynamic response of an individual to the other 9, 999 dynamic elements of the network. In the subcritical and critical domains the solutions could be modified to include truncations effects evident in the numerical data.
The lesson to be learned from the combination of computation and analysis presented herein is that complex networks of finite size whose dynamics are members of the Ising universality class, such as described by the DMM, have an analytic not just a numerical description. Instead of confining the dynamic description to that of the macroscopic variable, that being the global or average state of the network, one can also investigate how individual members of the network respond to the influence of the network as a whole. If we consider the fluctuations in the global dynamics to be microscopic, and the potential of the global variable to be macroscopic, then the real-time dynamic description of the individuals is mesoscopic. In general the mesoscopic dynamics can be described by a fractional stochastic differential equation.
Coupling two or more such fractional stochastic equations to model the across-scale coupling within the brain depicted in Figure 1 suggests itself. This is presently an active area of investigation.

V. ACKNOWLEDGEMENT
The authors would like to thank the U.S. Army Research Office for supporting this research. P.G. warmly thanks the ARO and the Welch Foundation for their support through Grants No. W911NF-11-1-0478 and No.B-1577, respectively.