Original Research ARTICLE
Fractional Dynamics of Individuals in Complex Networks
- 1Computational and Information Sciences Directorate, US Army Research Laboratory, Adelphi, MD, United States
- 2Information Sciences Directorate, US Army Research Office, Durham, NC, United States
The relation between the behavior of a single element and the global dynamics of its host network is an open problem in the science of complex networks. We demonstrate that for a dynamic network that belongs to the Ising universality class, this problem can be approached analytically through a subordination procedure. The analysis leads to a linear fractional differential equation of motion for the average trajectory of the individual, whose analytic solution for the probability of changing states is a Mittag-Leffler function. Consequently, the analysis provides a linear description of the average dynamics of an individual, without linearization of the complex network dynamics.
The last decade has witnessed the blossoming of two quite different strategies for the mathematical modeling of the complex systems, which are network science [1–3] and fractional calculus [4–6]. The widespread adoption of the network science perspective to study phenomena such as epidemic spreading of diseases , neuronal avalanches , or social dynamics  derives from the fact that these systems are composites of many simpler, interconnected, and dynamically interacting elements. Similarly, popularization of fractional calculus in research that concerns physical processes that are characterized by long-term memory and spatial heterogeneity [10, 11] stems from its particular mathematical formulation, based on a definition of the nonlocal differentiation and integration operators. Therefore, since memory effects and heterogeneity are frequently observed in biological, social, and man-made systems [12, 13], the application of fractional calculus in the domain of complex networks is a natural step toward providing novel analytical tools that are capable of addressing research questions arising in the field.
Despite the simplicity of their basic building blocks, complex systems, such as cooperative animal behavior , the flow of highway traffic , or the cascades of load shedding on power grids , are characterized by rich self-emergent behavior. However, since in most cases, solving a system of coupled nonlinear equations that trace the dynamics of a network composed of N units is not possible, the primary focus of investigations into complex networks has been on their global behavior . This approach follows the path taken by classical statistical physics, with Boltzmann's realization that the description of the state of a gas or a solid state could be only achieved on the scale of the entire system . Analogously, the ability to characterize the global behavior of a complex network comes at a price of not being able to quantify the dynamics of the components that give rise to it. Typically, one attempts to infer the global dynamics by averaging the behavior of single elements within the system, following a bottom-up approach of the mean field theory (see Figure 1).
Figure 1. Typical description of the dynamics arising from the interaction of numerous basic elements over a complex network that focuses on the global behavior of the system (left). Such an approach, however, comes with a price of not being able to quantify dynamics of individual elements within the system. In this paper, we address this problem by adopting statistical properties of the macroscopic dynamics in order to infer the behavior of individual units.
In this paper, we address this issue by posing the inverse question. Rather than inferring the global dynamics by combining the behavior of single elements within the dynamical system, we ask whether it is possible to construct a description of the dynamics of the individual elements, provided information about the network's global behavior. We approach the problem by considering statistical properties of the global variable.
Frequently, the macrovariables observed in complex networks display emergent properties of spatial and/or temporal scale-invariance. These are manifested by, for example, the inverse power scaling of waiting-time probability density functions (PDFs) between events, such as communication instances in human interactions or occurrence of earthquakes. At the same time, the inverse power laws (IPLs) that characterize the emergent macroscopic behavior are reminiscent of particle dynamics near a critical point, where a dynamic system undergoes a phase transition . However, despite the advances made by the renormalization group approach and self-organized criticality theories that have shown how scale-free phenomena emerge at critical points, the issue of determining how the emergent properties influence the microdynamics of individual units of the system remains open.
Herein, we address the problem of quantifying the response of an individual unit to the dynamics of the collective. This is done by taking advantage of the fractional calculus apparatus, whose utility arises from its ability to seamlessly incorporate the IPL statistics into its dynamics. The phase transitions that characterize many complex systems suggest the wisdom of using a generic model from the Ising universality class to characterize system dynamics. It is then possible to demonstrate that the individual trajectory response to the collective dynamics of the system is described by a linear fractional differential equation. This is achieved through a subordination procedure without the necessity of linearizing the underlying dynamics. Following this procedure, it is shown that the analytic solution to the linear fractal differential equation retains the influence of the nonlinear network dynamics on the behavior of the individual. Moreover, the solution to the fractional equation of motion suggests a new direction for designing mechanisms to control the dynamics of complex networks.
In section 2, we sketch out the mathematics of the dynamical decision making model (DMM), introduce renewal events, and subordinate the behavior of the individual to the mean field behavior of the network. In section 2.2, the dynamics of the individual is determined from the subordination theory to be a tempered fractional differential equation. The exact solution to this equation is given by an attenuated Mittag-Leffler function, which is fitted to the numerical solution of the DMM equation. In section 4, we discuss some implications of the high quality convergence of the analytical and numerical results of this complex network.
2. Complex Network Dynamics
As demonstrated by Grinstein et al. , any discrete system, defined by means of local interactions, with symmetric transitions between states and randomness that originate from the presence of a thermal bath or internal causes belongs to the universality class of kinetic Ising models. One such system is the DMM [21–23] and is the one we implement herein. Each individual unit si of the model is a stochastic oscillator and can be found in either of the two states, +1 or −1. The dynamics are defined in terms of the probability of an individual to be in either state, and it is modeled by the coupled two-state master equation,
where I and 1 are the 2 × 2 identity and unit matrices, respectively. The probability of being in one of the two states (+1, −1), , defines a Markovian telegraph noise, with symmetric and constant rate of changing states 0 < g0 < 1.
Positioning N such individuals at the nodes of a complex network introduces coupling between them [21, 22], which, here, is limited to the nearest neighbor interactions. The influence that unit si experiences due to the presence of its neighbors is expressed by a modification of its transition rate
which becomes a time dependent variable. Here, K is the strength of the coupling between nodes, 0 < K < ∞, constant for all nodes in the network. The variable M(i) denotes the degree of the node i, and denotes, respectively, the count of the nearest neighbors in states si(K, t) = 1 and si(K, t) = −1 at time t. As single units si change their states, quantities and fluctuate in time, while their sum is always conserved . In this paper, we consider the case of a regular two-dimensional lattice, where M(i) = 4 and for all the nodes. The single unit in isolation corresponds to the case of K = 0. When the coupling constant K > 0, a unit in state +1(−1) makes a transition to the state −1(+1) faster or slower according to whether or , respectively.
Time-dependent transition rates modify the two-state master (Equation 1) to take the form
where the matrix of rates Gi(t) is defined as
and p(i)(t) is the probability of the element i = 1, 2, …, N in the network at time t and is normalized such that for every i.
Dynamics of an entire network is described by a system of N such coupled equations, resulting in a highly nonlinear system , containing 6N dynamic variables . This number of coupled variables prevents the successful application of analytic methods, as these are usually adopted to solve problems that involve only a few coupled time-dependent differential equations. Instead, extensive numerical calculations are supplemented by an analytic formulation of the evolution of a global variable.
As depicted in Figure 2B, the global behavior of the model, defined by the fluctuations of the mean field variable
shows a pronounced transition as a function of the control parameter K. While in Figure 2A, the single elements appear to be essentially unchanged by their interactions with the rest of the network, the global variable shifts from a configuration dominated by randomness to one in which strong interactions give rise to long-lasting majority states shown in Figure 2B. Note that the origin of the random fluctuation in the DMM is the finite size of the network, which has nothing to do with the thermal fluctuations in the Ising model of magnetization.
Figure 2. Behavior of a discrete, two-state dynamic unit on a two-dimensional lattice. Temporal evolution and corresponding survival probability Ψ(τ) for the transitions between two states for the single unit si(t) of the system, presented on panels (A,C), respectively, are compared with the behavior and statistical properties of the global order parameter ξ(t), showed on panels (B,D). Simulations were performed on a lattice of size N = 50 × 50 nodes, with periodic boundary conditions, for g0 = 0.01 and increasing values of the control parameter K. Blue, red, and green lines correspond to K = 1.50, 1.70, and 1.90, respectively. The critical value of the control parameter is KC ≈ 1.72. Black dashed line on the plots of Ψ(τ) denotes an exponential distribution, with the decay rate g0.
To characterize the changes in the temporal properties of the micro- and macro-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 or the global variable, respectively. These calculations unveil modest deviations of Ψ(τ) for a single individual from the exponential form, Ψ(τ) = exp(−g0τ), that characterizes single non-interacting elements, as shown in Figure 2C. Clearly, the influence of the network on the behavior of the individual does not appear to induce a significant change in the latter. Despite such a modest change in the behavior of the individual, the global variable manifests IPL statistics, as depicted in Figure 2D. Thus, the following question arises: To what extent are individual opinions within a complex network influenced by the network dynamics?
2.1. Renewal Events
Many physical processes, for example earthquakes, radioactive decay, and social processes, such as making a decision, can be viewed as particular events. A characteristic property of an event is that it's onset can be precisely localized in time, even if its occurrence has extended consequences in space. Thus, the dynamics of a process characterized by events is described in terms of the probability of an event occurring, rather than by a more traditional Hamiltonian approach.
The process of event occurrence is characterized by the waiting-time PDF ψ(τ), which specifies the distribution of times between consecutive events. The probability for an event to occur in the short time interval [t, t + dt] is given by
where τ is measured from the occurrence of the previous event. Consequently, one can define the survival probability Ψ(τ) as the probability that no event occurs up to the time since the last event as
As a consequence of this integral, the waiting-time PDF can be written as
and the PDF ψ(τ) is a properly normalized function,
since it is assumed that an event occurs somewhere within the time interval (0, ∞). It is also true that no event occurs at time t = 0, which means that the survival probability Ψ(0) = 1.
A particular class of events can be defined, renewal events, that reset the clock of the system to an initial state instantaneously after their occurrence. After a renewal event takes place, the system evolves in time independently of whatever occurred earlier, having no memory of previous instances in which such an event occurred. Some examples of renewal events found in physics include anomalous diffusion of tagged particles inside living cells, blinking quantum dots, and defects arising in the weak turbulence regime of liquid crystals.
The renewal character of events is captured by the probability of n events occurring as follows. First, one assumes that an event occurs at time t = 0, thus, ψ0(t) = δ(t). Next, the first event occurs at time t > 0, taking place with the probability ψ1(t) = ψ(t). Subsequently, the probability for event n in a sequence to occur at time t is expressed in terms of probabilities of earlier events by the correlation chain condition
Frequently, experimentally observed waiting-time PDFs are exponential, but quite often in complex networks they are IPLs. For the purpose of this paper, we define the waiting-time PDF in terms of the hyperbolic distribution
If the events are generated by an ergodic process, then μ > 2, and the first moment of the hyperbolic PDF is
In the framework of renewal theory, Equation (11) denotes the average time that one would have to wait between successive events. However, when μ < 2, the process is non-ergodic, and the mean value of the distribution diverges. In the non-ergodic case, T becomes a characteristic time scale of the process.
2.2. Subordination of Time
The notion of different clocks associated with different physical systems arises naturally in physics; the linear Lorentz transformation in relativistic physics being probably the most familiar example. Thanks to the recent availability of time-resolved data, biological, and social sciences have also started adopting the notion of multiple clocks, distinguishing between cell-specific and organ-specific clocks in biology and person-specific and group-specific clocks in sociology. Of course, the notion of subjective and objective time dates back to the middle of the nineteenth century with the introduction of the empirical Weber-Fechner law .
However, the striking difference between the clocks of classical physics and natural sciences is that the relations between the latter clocks are nonlinear. While the global activity of an organ, such as the brain or the heart, might be characterized by quite regular, often periodic fluctuations, the activity of single neurons demonstrates burstiness and noisiness. Similarly, in a society, people operate according to their individual schedules, not always being able to perform particular actions in the same global time frame. Thus, owing to the stochastic behavior of one or both clocks, a probabilistic transformation between times is necessary. An example of such a transformation is the subordination procedure.
We begin by defining two clocks. The first clock records a discrete operational time n, which measures the time T(n) of an individual. The second clock records the continuous chronological time t, which measures the time T(t) that a system of individuals have agreed upon. If each advancement of the discrete clock n is thought of as an event, then the relation between the operational time and chronological time can be given by the waiting-time PDF of those events in chronological time ψ(t). Assuming a renewal property for events, as given by chain condition from renewal theory (Equation 9), one can relate operational time to chronological time by
Every advancement of the operational clock is an event, which in the chronological time occurs at time intervals drawn from the renewal waiting-time PDF. Because of this randomness, one needs to sum over all events, and the result is an average over many realizations of the transformation.
As an example, consider the behavior of a two-state operational clock, whose evolution is shown in Figure 3. In operational time, the clock switches back and forth between its two states at equal unit time intervals. In chronological time, however, this regular behavior is significantly distorted. In the figure, the time transformation was taken to be an IPL PDF of waiting times. Thus, a single time step in the operational time corresponds to a time interval being a random number drawn from ψ(t) in chronological time. The long tail of the IPL PDF leads to especially strong distortions of the operational time trajectory, since there exist a non-zero probability of drawing very large time intervals between events. However, since the transformation between the operational and chronological time scales involves a random process, one needs to consider infinitely many trajectories in the chronological time, which leads to the average behavior of the clock in the chronological time denoted in Equation (12) by the bracket.
Figure 3. The upper curve is the regular transition between the two states of the individual in operational time. The lower curve is the subordination of the transition times to an IPL PDF to obtain chronological time.
We note that the time subordination procedure can also be used to model communication delays in the system. However, contrary to frequently used approaches, where individual units of the system are subordinated to model the interaction delay, here, we adopt the statistics of the macroscopic variable to derive the behavior of the interacting individual units. The coupling between units causes them to deviate from the Poisson behavior of an individual non-interacting unit. However, as illustrated in Figure 2, the time scale of interacting units is orders of magnitude that are smaller than the time scale of the macroscopic variable. Thus, we use the statistical properties of the macroscopic variable to provide a first-order estimate of the single unit dynamics. As such, we adopt a top-down approach, which is different from the bottom-up approach adopted for the consideration of communication delays.
3. Complex Network Subordination
To determine the network's influence on the dynamics of the individual, we adapt the subordination argument of the preceding section and relate the time scale of the macro-variable ξ(K, t) to the time scale of the micro-variable si(K, t). The two-state master equation for a single isolated individual in discrete time n in steps of Δτ is
where the notations φ(n) = φ(nΔτ) and φ = p1 − p2 depict the difference in probabilities for the typical individual to assume one of the two states. The solution to this discrete equation is
which, in the limit g0 Δ τ < < 1, becomes an exponential. However, when the individual is a part of a network, the dynamics are not so simple.
Adopting the subordination interpretation, we define the discrete index n as an individual's operational time that is stochastically connected to the chronological time t, in which the global behavior is observed. We assume that the chronological time lies in the interval (n − 1) Δτ ≤ t ≤ n Δτand, consequently, the equation for the average dynamics of the individual probability difference is given by 
Here, the time t in the waiting-time PDF ψ(t) is determined from the derivative of the survival probability. The empirically determined analytic expression for the survival probability is
The dominant behavior of the empirical survival probability is an IPL as indicated in Figure 2D. However, at early times, the probability of not making a transition approaches the constant value of one; at late times, the probability of not making a transition at a given time decays exponentially; it is in the middle range, where the probability is an IPL. The extent of the IPL range of the survival probability is determined by the empirical values of T, μ, and ϵ, and from Figure 2D, the value of ϵ is seen to become smaller as the control parameter K increases. The IPL functional form of the PDF results from the behavior of the survival probability Ψ(τ) of the global variable depicted in Figure 2D, with μ = 3/2.
Using a renewal theory argument, Pramulkkul et al.  show that Equation (15) expressed in terms of Laplace transform variables indicated by for the time-dependent function f(t) has the form
Note that u is replaced by u + ϵ in the Laplace transforms, because the exponential truncation of the empirical survival probability shifts the index on the Laplace transform operation. The asymptotic behavior of an individual in time is determined by considering the waiting-time PDF as u → 0,
so that Equation (17) reduces to
The inverse Laplace transform of Equation (20) yields the tempered rate equation
where the operator is the Caputo fractional derivative for 0 < α = μ − 1 < 1  and
Note that owing to the dichotomous nature of the states, 〈φ(t)〉 is the average opinion of the individual si(K, t).
The solution of the asymptotic fractional master equation (Equation 21) for a randomly chosen unit within the network is given by an exponentially attenuated Mittag-Leffler function (MLF):
and the MLF is defined by the series
The MLF is a stretched exponential at early times and an IPL at late times, with α = μ − 1 being the IPL index in both domains.
3.1. Comparisons With Numerics
We test the above analysis with numerical simulations of the dynamic network on a two-dimensional lattice with nearest-neighbor interactions in all three regions of DMM dynamics: subcritical, critical, and supercritical. The time-dependent average opinion of a randomly chosen individual is presented in Figure 4, where the average is taken over 104 independent realizations of the dynamics in the subcritical, critical, and supercritical regimes.
Figure 4. The probability difference 〈φ(t)〉 estimated as an average over an ensemble of 104 independent realizations of single element trajectories. Each trajectory corresponds to evolution of a randomly selected node within a N = 100 × 100 lattice network, with g0 = 0.01 and the same initial condition si(0) = 1. The parameter values for the numerical data are given in Figure 2 and from left to right K = 1.0 (A), 1.7 (B), 2.5 (C), respectively. The fit of the exponentially truncated MLF to the numerical calculations is summarized in Table 1.
A comparison with the exponential form of 〈φ(t)〉 for an isolated individual indicates that the influence of the network on the individual's dynamics clearly persists for increasingly longer times with increasing values of the control parameter within the network. The parameters μ and λ of Equation (23) obtained through fitting numerical results of Figure 4 with the MLF are summarized in Table 1. It is evident that the influence of the network dynamics on the individual is greatest at long times. The deviation of the analytic solution from the numerical calculation is evident for values of the control parameter at and below the critical value. The analytical prediction is least reliable at extremely long times in the subcritical domain. Consequently, the response of the individual to the group mimics the group's behavior most closely when the control parameter is equal to or greater than the critical value.
Herein, the subordination procedure provides an equivalent description of the average dynamics of a single individual within a complex network, in terms of a linear fractional differential equation. The fractional rate equation is solved exactly, determining the Poisson statistics of the isolated individual becomes attenuated Mittag-Leffler statistics, owing to the interaction of that individual with the other members of a complex dynamic network.
Consequently, an individual's simple random behavior, when isolated, is replaced with behavior that might serve a more adaptive role in social networks. We conjecture that the behavior of the individual is generic, given that the DMM network dynamics belong to the Ising universality class. Members of this universality class share the critical temporal behavior  that drives the subordination process. It is the renewal property of the event statistics, which, through the subordination process, gives rise to the linear fractional master equation for the typical individual's dynamics. The solution to the tempered fractional rate equation manifests the subsequent robust behavior of the individual; it remains to be determined just how robust the behavior of the individual is relative to control signals that might be driving the network.
As pointed out by Liu et al. , the ultimate understanding of complex networks is reflected in the ability to control them. Recent observations of the interconnectedness of infrastructure networks , facilitating the spread of failures  or the tight coupling between banking institutions, posing a danger to the stability of global financial system , demonstrate the importance of developing a systematic approach to influence and/or control the complex networks. The analysis presented here provides an alternative attempt to address this need directly. Subordination suggests a way to impose the conditions of traditional control theory  onto the complex network dynamics by, first, expressing the underlying nonlinear network dynamics in the form of a linear fractional equation of motion. This approach at addressing control will be pursued in a future publication.
BW developed the theoretical formalism and performed the analytic calculations. MT performed numerical simulations. Both authors discussed the results and contributed to the final manuscript.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors acknowledge the US Army Research Laboratory for supporting this research.
1. ^ Adjusted goodness of fit, , is defined as the ratio of the sum of squared residues for the nonlinear fit with the MLF (SSreg) and for the fit to the average value of data points (SStot), where n is the number of data points and κ is the number of free parameters being estimated.
Keywords: fractional calculus, subordination, inverse power law, complex networks, control
Citation: Turalska M and West BJ (2018) Fractional Dynamics of Individuals in Complex Networks. Front. Phys. 6:110. doi: 10.3389/fphy.2018.00110
Received: 25 April 2018; Accepted: 10 September 2018;
Published: 16 October 2018.
Edited by:Dumitru Baleanu, University of Craiova, Romania
Reviewed by:Carla M. A. Pinto, Instituto Superior de Engenharia do Porto (ISEP), Portugal
Yilun Shang, Tongji University, China
Copyright © 2018 Turalska and West. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Malgorzata Turalska, firstname.lastname@example.org