Lévy walk dynamics explain gamma burst patterns in primate cerebral cortex

Lévy walks describe patterns of intermittent motion with variable step sizes. In complex biological systems, Lévy walks (non-Brownian, superdiffusive random walks) are associated with behaviors such as search patterns of animals foraging for food. Here we show that Lévy walks also describe patterns of oscillatory activity in primate cerebral cortex. We used a combination of empirical observation and modeling to investigate high-frequency (gamma band) local field potential activity in visual motion-processing cortical area MT of marmoset monkeys. We found that gamma activity is organized as localized burst patterns that propagate across the cortical surface with Lévy walk dynamics. Lévy walks are fundamentally different from either global synchronization, or regular propagating waves, because they include large steps that enable activity patterns to move rapidly over cortical modules. The presence of Lévy walk dynamics therefore represents a previously undiscovered mode of brain activity, and implies a novel way for the cortex to compute. We apply a biophysically realistic circuit model to explain that the Lévy walk dynamics arise from critical-state transitions between asynchronous and localized propagating wave states, and that these dynamics yield optimal spatial sampling of the cortical sheet. We hypothesise that Lévy walk dynamics could help the cortex to efficiently process variable inputs, and to find links in patterns of activity among sparsely spiking populations of neurons.

L évy walks are anomalous diffusive forms of random motion, which have been widely observed in natural systems including movements of humans and other animals [1][2][3][4][5][6] . In a Lévy walk, clusters of short-step sizes are occasionally interspersed with longer movements. A Lévy walk therefore represents a superdiffusive process, and can travel much further from its starting position than a Brownian walk of the same duration ( Fig. 1). This characteristic enables Lévy walks to efficiently implement complex tasks such as searching for sparsely distributed resources 4,7 . Recent evidence showed that Lévy walk dynamics can be intrinsically generated by neural circuits in Drosophila larvae 8 , but the question whether Lévy walks are a characteristic of vertebrate brain activity remains open.
Lévy walks dwell near one location for a while and then intermittently switch to new locations, leading to a "heavy-tail" probability distribution of movement step lengths (Fig. 1b). This heavy-tail property naturally generates a bursting characteristic in time and space (Fig. 1a). In the vertebrate cerebral cortex, temporal bursts of high-frequency (gamma-band) brain oscillations  have been reported, under both awake and anaesthetized conditions [9][10][11][12][13][14] . Thus the temporal properties of gamma bursts raise the possibility that they may be characterized by Lévy walks, but the spatiotemporal organization of gamma bursts remains unexplored.
Here, we first present an empirical analysis of gamma bursts in a visual motion-processing area (area MT) of marmoset cerebral cortex, under anesthesia and in absence of patterned visual stimulus (resting state). Under these conditions we find that gamma bursts typically appear as localized patterns, which propagate across the cortical surface. The propagation dynamics can be explained by Lévy walks. The occasional large displacements characteristic of Lévy walks yield wave-like propagations, and the intervening short-step clusters give rise to transient local patches of gamma synchrony. The Lévy walk dynamics of gamma burst patterns thus represent a novel mode of neural population activity.
To explain our experimental findings, we employ a spatiallyextended cortical circuit model of spiking neurons 15 to show that the Lévy walk dynamic of gamma bursts is produced by transitions between asynchronous and regular wave cortical states. We show that this mechanism allows functional interactions in the cortex across many spatial scales, without large increases in spiking rate.
The high metabolic expense of spiking activity in the brain 16,17 means that cortical spikes are a sparsely distributed and scarce resource. We propose that Lévy walk dynamics of gamma bursts enable efficient functional linking of sparsely distributed spikes across the cortical surface, and rapid state switching to process rapid changes in the locations and timing of inputs to cortical circuits.

Results
Spatiotemporal dynamics of gamma bursts. We recorded local field potentials (LFPs) from the middle temporal (MT) cortical area of four sufentanil-anesthetized adult marmosets using multielectrode arrays (10 × 10 electrodes; see Methods), while the animal viewed a spatially uniform monitor screen (Fig. 2a). As expected 18 , the time-averaged power spectrum of the LFP showed bias to low temporal frequencies. Inspection of the LFP recordings at finer temporal scales, however, showed distinct bursts of power in the gamma band (30-90 Hz; Fig. 2d). These bursts were identifiable in the raw LFP traces and in the time-frequency spectrograms, computed by either bandpass filtering or wavelet decomposition ( Fig. 2b-d). The burst properties are consistent with those described for spontaneous activity in mouse barrel cortex and macaque visual cortex 13,14,19 , insofar that bursts at individual electrodes exhibited variable duration and peak frequency (Fig. 2d).
The gamma bursts observed at one electrode were usually accompanied by similar bursts at nearby electrodes, and we found that the gamma bursts were organized as spatially localized patterns. These localized patterns emerged at seemingly random   locations and moved across the recorded area with complex spatiotemporal dynamics; Fig. 3a shows such patterns starting at different timepoints. The localized patterns dwelled near one location for a while and then quickly shifted to another location in an intermittent way. To characterize these movement dynamics, we developed a method to track the gamma bursts (see Methods). At each time moment t, we calculated the center of mass (CoM) position, X t ¼ x t ; y t À Á , of the pattern (see Methods). The positions of the CoM at sequential timepoints provide the movement trajectory of a burst pattern (two examples are shown in Fig. 3a). These movement trajectories show superdiffusive dynamics and heavy-tailed, power-law distributions of step lengths, which are two key properties of Lévy walks.
Superdiffusive dynamics of gamma bursts. To determine whether these trajectories of gamma burst patterns were consistent with the superdiffusive property of Lévy walks, we calculated the , as a function of time increment τ. The MSD was linear on a log-log scale (Fig. 3b), indicating that that MSD is a power function of time increment τ, such that r τ ð Þ / τ β . Brownian motion is characterized by MSD with β = 1; β > 1 indicates a superdiffusive process and β < 1 indicates a subdiffusive process 20 . We estimated the value of β by minimizing the square error between the observations and this model using a Levenberg-Marquardt algorithm (see Methods); 2648 trajectories analyzed in all 4 animals showed a good fit to the model with the residuals < 0.1. The typical trajectory shown in Fig. 3b is a superdiffusive process with β = 1.24. The gamma bursts exhibited variable propagation dynamics: some patterns were very mobile while others were more stationary, with the diffusion exponent β ranging from 0.3 to 1.9 for all well-fit trajectories (for example, case MA026 in Fig. 3c). Most trajectories appeared superdiffusive: the mean diffusion exponent was β ¼ 1:42 (s.d = 0.24, n = 773, case MA026). The diffusion exponents for other animals and the aggregated data showed similar distributions with the mean diffusion exponents > 1.0 as well ( Supplementary Fig. S1). This superdiffusion suggests that long-range correlations underlie the seemingly random motion of gamma burst patterns 20 .
Distributions of movement step lengths. We found that the step lengths that form the trajectories of gamma bursts were also consistent with Lévy walk.
Step length was defined as the distance travelled between two points that were preceded and followed by either a pause or a change in direction 3,21,22 . We used two methods to detect a change in movement direction: a onedimensional (1D) model [22][23][24] and an angle model 3,21 . As in Rhee et al. 3 , we then characterized the distribution of step length with a complementary cumulative probability distribution (CCPD), which is better than a simple probability distribution at revealing the tails of distributions. In the 1D model, we projected the 2D trajectories onto 1D and defined a turning point as a reversal of movement direction 22,24 . The 1D method is unbiased because it is not dependent on choice of criterion turning angle, while preserving the distribution properties of the original 2D trajectory 24 . We found that the resultant CCPDs of the step lengths of all four animals showed heavy tails ( Fig. 3d-g, and the aggregated data shown in Supplementary Fig. S2). In the angle model, we  identified turning points as sequential time steps when the change in movement angle was larger than a criterion turning angle, θ (Fig. 3h). We varied θ in 20 increments from 20 to 120 and found that the angle method yielded the similar heavy-tailed distributions as in the 1D model ( Supplementary Fig. S3).
The measured trajectories of the gamma burst patterns are naturally bounded by the recording area, which restricts the distribution of movement step-lengths, such that only truncated Lévy walks are experimentally plausible. We therefore tested the CCPD of movement step lengths against a truncated Lévy distribution (that is, a truncated power-law model). We first used Maximum Likelihood Estimation (MLE) to fit five models to the data, including a truncated power-law, exponential, normal, lognormal (a heavy-tailed distribution) and gamma functions (these functions are detailed in Table 1). We then used Akaike Information Criteria weights (AIC) for model comparison and selection 3,25 (see Methods). We examined distributions of step lengths for individual animals and calculated the AIC weights; Table S1 shows the AIC weights in the 1D model and the angle model with different parameters. In all cases, the truncated power-law distribution had the largest AIC weight among the tested distributions, indicating that it provided the best characterization of the step-length distributions. The exponents of the truncated power-law distribution were all in the Lévy range 1<λ ≤ 3, with λ = 1.32, 1.33, 1.  Table S1). The explanatory power of the truncated Lévy distribution was not sensitive to the threshold used for gamma burst pattern detection (1.5 s.d. to 2.5 s.d.) ( Supplementary Fig. S4), the minimal step length, or the spatial filtering (Supplementary Note, Fig. S5). These results therefore indicate that the movements of gamma burst patterns were consistent with Lévy walks. We also calculated the movement speeds of the gamma patterns and found that they were very variable with a broad distribution. The average speed was 0.376 m/s (s.d. = 0.309, n = 2.4 × 10 6 movements), consistent with the propagation speed of gamma patterns found in the visual cortex of rabbit 26 .
We performed parallel analyses on surrogate datasets. One set was generated by randomizing the phases of the entire LFP recordings in the Fourier domain-the Fourier amplitudes were retained to ensure the same autocorrelation as the original data 27 . The same analysis from the raw data to the gamma burst patterns was applied. The numbers of bursts at each electrode were similar in the experimental data (mean: 1317 over all electrodes and all animals, s.d.: 177, n = 4 animals) and the surrogate data (mean: 612, s.d.: 19, n = 4 animals). However, propagating burst patterns were rarely observed in the surrogate data (n p < 10 in all recordings), where n p is the number of patterns, whereas propagating bursts were common in the original data (mean: 1785, s.d.: 190, n = 4 animals, p < 10 À3 ). This result rules out the possibility that the propagating gamma burst patterns are random events. The second surrogate analysis was done after the gamma patterns were detected. We shuffled the CoM position within each detected trajectory (mean duration: 79.13 ms, s.d.: 107.31 ms, n = 7027 for four animals) and then applied the procedure described above to study the step length distribution; the distribution of step length of the surrogate data follows a normal distribution (Fig. 3i). These results indicate that the Lévy walk dynamics is an inherent property of gamma burst patterns.
Gamma burst patterns in model circuits. To understand what circuit mechanisms could underlie the Lévy walk dynamics of gamma burst patterns, we implemented a biophysically plausible, spiking circuit model with excitatory and inhibitory neurons that capture the known anatomy and physiological of cortical circuits 15 . It incorporates distance-dependent synaptic connectivity and balanced excitation and inhibition 15 (Methods). The model exhibits a rich repertoire of dynamic activity, ranging from asynchronous to propagating wave states (Fig. 4). If the inhibition-to-excitation (I-E) ratio ξ is small (ξ < ξ c ; ξ c = 3.4), localized propagating waves emerge and propagate across the neural circuit with relatively smooth and regular moving trajectories (Fig. 4b). On the other hand, if the I-E ratio is large (ξ > ξ c ), the circuit exhibits an asynchronous state without any structured patterns in spiking activity (Fig. 4d). When the circuit is close to the transition region between the asynchronous state and the wave state (ξ % ξ c ), localized spiking activity patterns emerge intermittently and exhibit complex spatiotemporal dynamics ( Fig. 4c) consistent with the Lévy walk dynamics of gamma bursts found in our experimental data.
where ΓðaÞ is the Gamma function  Fig. 5a). Individual neurons in the model fire variably with mean coefficient of variation of interspike intervals 0.96; these variable spikes are phase locked to gamma bursts (Fig. 5b), consistent with empirical observations 9,14 . Further, the amplitudes of both excitatory and inhibitory postsynaptic currents onto individual neurons exhibit large cycle-to-cycle fluctuations (Fig. 5c, left). Despite these fluctuations, excitation and inhibition in the circuit model are tightly balanced (Fig. 5c, left), as measured in experimental studies 28 . The I-E ratio is~3.4; it is interesting to note that this value and the presence of large fluctuations are quantitatively consistent with experimental studies of gamma oscillations 29 . The fluctuations form a skewed distribution with a heavy-tail (Fig. 5c, right). We used the MLE method to fit the Lévy stable, normal and exponential functions to the distribution. The probability density function of Lévy stable where γ is the scale parameter, c α = sin(πα/2)Γ(α)/π, Γ is the Gamma function, α s is the skewness parameter and α is the stability parameter with the range 0 < α < 2 30 . The AIC weights showed that the Lévy stable distribution provided the best characterization. The stability parameter α of the fitted Lévy stable distribution is 1.216 (confidence interval α 2 [1.214, 1.218]). Such large fluctuations only occur in the transition regime of the circuit. In the asynchronous state, the synaptic inputs stabilize and have Gaussian-like statistics. In the propagating wave state, the inputs become semi-periodic and regular. The synaptic fluctuations of our model resemble the large fluctuations typically observed in complex physical systems near phase transitions of different states 31 .
Consistent with empirical observations 29 , fluctuations in amplitudes of both excitatory and inhibitory currents in our model are correlated with the amplitudes of LFP gamma oscillation on a cycle-by-cycle basis (excitatory: r = 0.43 ± 0.27, p < 10 −3 ; inhibitory: r = 0.55 ± 0.15, p < 10 −3 ). The amplitude of each gamma oscillation cycle is strongly correlated with the latency to the subsequent cycle (r = 0.76 ± 0.01, p < 10 −3 ) with higher amplitude leading to longer latency of the next cycle; this is a characteristic feature of gamma oscillations in rat hippocampus 29 and awake macaque area V1 32 . The duration of each gamma cycle is also inversely correlated with the average spike rate (r = −0.54 ± 0.01, p < 10 −3 ), as found by Spyropoulos et al. 32 .
We next tracked and analyzed the trajectories of the localized activity patterns produced by the model, using the same methods we used to analyze the experimental data (Fig. 6a). Figure 6b shows the MSD of a typical burst pattern, which is a power function of time increment τ, d τ ð Þ / τ β , with β = 1.56. The distribution of diffusion exponents in the model is likewise comparable to that found in our experimental data (Fig. 6c), with the diffusion exponent β ranging from 0.48 to 1.98 (mean β = 1.25). We detected the turning points of movement trajectories as described above and used these points to define the distributions of step lengths, again characterizing the CCPDs using MLE to derive AIC weights and determine the best model (see Supplementary Table S2). We found that the CCPDs for both the 1D and angle methods were best characterized by a truncated power law distribution with exponents in the Lévy range 1<λ ≤ 3 ( Fig. 6d and Fig. 6e), as in the experimental data.
Shifting the model network away from the dynamic regime and further into the asynchronous regime led to a loss of localized activity patterns (Fig. 4d). By contrast, shifting the network further into the wave regime resulted in a regular propagation of activity patterns (the red curve in Fig. 4b indicates a smooth trajectory) lacking the spatiotemporal dynamics characteristic of the transition regime. It is also interesting to note that in the critical transition regime, the pattern propagation speed was highly variable (Supplementary Fig. S6) as found in our data. However, the mean of the propagation speeds was 11.2 μm/ms, which was much smaller than that found in our data. This discrepancy is likely because the total number of neurons of our network is much smaller than real neural circuits. In the excitation-dominated state of our model, the localized wave patterns propagated more slowly (~6.82 μm/ms) than the patterns in the critical transition regime and did not exhibit significant variability.
Metabolic efficiency of activity patterns sampling cortical space. To test for a functional advantage of the activity pattern with Lévy walks in the transition regime, we used the neural circuit model to calculate a spike energy budget while the pattern moves over the cortical surface. To quantify the sampling capacity of the spiking pattern, we uniformly divided the 2-dimensional network into N t ¼ 100 small areas (i.e., squares) and calculated the number of different squares (N s ) visited by the CoM within 3 s. The sampling rate is To calculate neural energy expenditure we used the method of Lévy and Baxter 33 : in each time window 4τ ¼ 50 ms, each neuron costs r energy units due to leak currents and each spike costs one extra unit of energy, so the energy expenditure is: where m is the number of spikes during the time interval, n is the total number of neurons in the network. Based on Eq. (1), we then obtained the total energy expenditures E total . The sampling efficiency of the activity pattern is defined as the sampling rate S r divided by the total energy expenditure, η ¼ S r E total . The searching efficiency indicates that the larger the value η, the more metabolically efficient it is for the pattern to sample the 2D cortical space. In our calculation we used r = 0.01; 4t ¼ 50 ms and N t ¼ 100 to calculate the sampling efficiency; other choices near these values yielded similar results. In our circuit model, we found that as the I-E ratio (ξ) decreases, the sampling rate S r increases (Supplementary Fig. S7) and that the propagating wave regime (State I) provides higher overall S r than other regimes. The energy expenditure E total also increases as ξ decreases (Supplementary Fig. S7). However, as shown in Fig. 4a (green line) the sampling efficiency η (η ¼ S r E total ), which is the trade-off between the sampling rate and the energy expenditure, is maximal when the network is in the critical transition region (State II, where the I-E ratio is~3.4). The sampling efficiency in State II is significantly larger than the values of η in either State I or State III (p < 10 À3 , Kolmogorov-Smirnov test). This result thus indicates that the activity pattern with Lévy walk dynamics are an energy-efficient way to dynamically sample the cortical space.

Discussion
We find that the spatiotemporal organization of gamma bursts in primate cortex can be characterized as Lévy walks. Our circuit model explains how the Lévy walk dynamics can arise from intrinsic, near-critical transitions between cortical states. Thus, gamma bursts exhibit more complex spatiotemporal dynamics than expected by conventional views of stable steady states [34][35][36][37] or filtered Gaussian noise 38,39 . In the following we first relate the Lévy walk dynamics revealed here to previous studies of gamma bursts, then consider the origin of these dynamics in synaptic activity as illustrated by our circuit model. Finally, we consider implications for cortical function.
Spatial and temporal properties of gamma bursts. The Lévy walk property means that localized gamma burst patterns hover around one location for a while and then move or jump to another location in an intermittent manner. These spatiotemporal dynamics can reconcile and extend some previous findings on gamma oscillations. For instance, spatially localized gamma oscillations were reported by Freeman and Barrie 26 and Sirota et al. 40 but their spatiotemporal evolution was not studied. The Lévy walk characteristic of gamma bursts means they occasionally exhibit long propagation trajectories; this property would yield a wave-like propagation as previously described [41][42][43] . At the same time, the intermittent dwelling of gamma bursts gives rise to gamma synchrony (correlation) of neurons in the local cortical area 34,35 ; we show here that such gamma synchrony is localized in space but is transient in time.
Other work has shown that the magnitude of gamma oscillation is modulated by slower oscillations such as theta oscillation 44,45 . Our findings imply that such cross-frequency coupling would also exhibit complex spatiotemporal dynamics. If propagating patterns at different frequencies are coupled across multiple spatial and temporal scales, the resulting cascades of pattern interactions could be comparable to cascaded interactions in coherent structures such as vortices in turbulent fluids 46 .
Overall, our results show that the simple concept of temporal correlation, which has guided our understanding of gamma oscillations during the past two decades 34,47 , can usefully be expanded to take spatial variables into account; such expansion can provide a new perspective for understanding gamma oscillations. Modeling Lévy walk dynamics of gamma bursts. Our modeling indicates that Lévy walk dynamics of gamma burst patterns depends on critical transitions between asynchronous states and propagating wave states. Lévy walks arising from criticality have been mainly found in models of low-dimensional dynamical systems. For example, one recent study 48 showed Lévy walks can by generated by coupled chaotic oscillators near the transition between synchronous and asynchronous states. In ref. 49 , bursty time-series of coupled Lotka-Volterra (L-V) models and spike potentials recorded from rat hippocampus were used to generate superdiffusive Lévy processes. Complex neural dynamics have been characterized as critical avalanches 50 , but it is unclear whether these avalanches could be organized as localized activity patterns with Lévy movements in space. In contrast, our model is neural circuit-based and incorporates properties such as distancedependent synaptic coupling 15 , synchronized excitatory and inhibitory synaptic fluctuations (as measured empirically by Atallah and Scanziani 29 ), and localized propagating patterns.
In classical neural circuit models, in contrast to our model, neurons fire in regular fashion, so these classical models are unable to capture variable spike rates and non-stationary bursts 36,37,51 . Population gamma oscillation and variable spike rates can be reconciled in the sparse synchrony regime of randomly coupled circuits 47,52 . But synaptic inputs in sparse synchrony models are characterized as Gaussian noise 53 , which is incompatible with the heavy-tailed distributions that we observe experimentally and in our circuit model. Heavy-tailed, non-Gaussian properties often emerge from complex non-equilibrium systems, suggesting that methods for analyzing such complex systems (for example, fractional Fokker-Plank formalisms 54 ) could be a promising direction to pursue for formal analysis of the Lévy walk dynamics.
Relevance to brain function. What could be the relevance of our results for understanding information processing in the brain? First, it is natural to ask whether the results we found here under anesthesia would apply to the waking state. In the waking state, gamma oscillations are associated with functions including sensory processing 59 , cognition, memory 35,44,60 and attention 61,62 . In the anesthetized state, and at rest in absence of patterned sensory stimulus (resting state: Fox and Raichle 63 ) the amplitude and coherence of gamma oscillations is decreased but their dynamic properties are otherwise preserved 9,64 . Likewise the main change in subcortical visual evoked spike activity under sufentanil anesthesia is reduced response amplitude without other changes in receptive field properties 65 . Finally, resting-state fluctuations (in terms of up and down changes over time) in blood-oxygenlevel-dependent fMRI signals are preserved at anesthetic levels that produce profound loss of consciousness (0.8-1.5% isoflurane) 66 . These examples of preserved elementary brain functions under anesthesia predict that basic properties of Lévy walk dynamics in gamma bursts will be preserved in the waking state. Such experiments are beyond the scope of the present study, but it would be straightforward to apply our circuit model to recordings from waking brains.
Lévy walk movements have been observed of Drosophila larva neural circuits 8 and in other biological systems and complex physical systems 1,6,20 . For instance, Lévy walks are essential for optimally transporting energy in turbulent fluids 67 , for animals to optimally search for spatially distributed food 7 , and for lymphocytus T-cells to efficiently find target pathogens in brain explants 4 . In the context of normal brain function we propose here two hypotheses, which are not mutually exclusive, as follows.
Our first hypothesis relates to the fact that gamma bursts have been widely observed in attentional and cognitive tasks. The Lévy walk dynamics of gamma bursts may help an attentional "spotlight" to focus on one location for a while and then switch/jump to another location in an intermittent manner. Relatedly, the scale-free property of Lévy walks could provide an efficient mechanism for sampling natural environments, which are inherently scale-free 68,69 . It has been found that saccadic shifts when viewing natural scenes follow Lévy motion 70 . Attention and saccadic eye movement are closely related, as visual spatial attention determines the end point of saccades 71 . Taken together, these observations suggest the Lévy motion property of gamma bursts may improve the efficiency of attention sampling. This hypothesis could be tested by neural population decoding 72 , and comparing the statistical properties of spatial attention to predictions of Gaussian and Lévy walk dynamics. Based on our results (Fig. 3, Fig. 6) these predictions would diverge for temporal scales above 100 ms and spatial scales represented by 100-1000 µm on the surface of area MT, at least in marmoset monkeys. The question whether the same scale factors hold in larger brains or unanaesthetized brains remains to be tested.
Our second hypothesis relates to the metabolic expense of spiking activity in the cerebral cortex 16,17 and the associated evidence for sparse spiking as a basis for information processing in the brain [73][74][75] . Our modeling (Fig. 4) showed that activity patterns with Lévy walks dynamics are a metabolically efficient way to visit spatially distributed locations on the cortical sheet; that is, the ratio of the locations visited to the energy consumed is maximal in the transition regime where Lévy walks dynamics emerge. In this context, the superdiffusion dynamics of gamma bursts would enable widely separated regions of cortex to be activated within short time steps. By analogy, cortical activity patterns could be seen as efficient "foraging" for spikes which are separated in space but close together in time. This hypothesis could be tested by comparing Gaussian to Lévy walk statistics of gamma burst activity, at the millisecond timescales relevant to spike-dependent plasticity 76 and related phenomena.

Materials and methods
Multielectrode array (10 × 10 electrodes, 1.5 mm length, electrode spacing 400 µm, Blackrock Microsystems) recordings were made from the MT cortical area of four adult male marmosets (Callithrix jacchus). Procedures conformed to the Australian National Health and Medical Research Council Code of Practice for the use and care of animals, and were approved by institutional committee at the University of Sydney. Details of the preparation are given elsewhere 77 . Anesthesia and analgesia were maintained by intravenous sufentanil infusion (6-30 μg Á kg À1 h À1 ) and inspired 70:30 mix of N 2 O and carbogen. The eyes were held open with microretractors and protected with oxygen-permeable contact lenses. Refraction was optimized (with optometric trial lenses) by maximizing the response of the first recorded units to high spatial frequency drifting gratings. Stimuli were viewed through natural pupil, which typically had diameter close to 2 mm. The animal viewed a uniform gray field presented on a cathode-ray-tube monitor (Sony G500, refreshed at 100 Hz, viewing distance 45 cm, mean luminance 45-55 cd m −2 ). Recording surface insertion depth was targeted to 1 mm. The LFP sampling frequency was 1024 Hz. Recording duration ranged from 5 to 25 min in total: the first 5 min of recording were analyzed for each animal (n = 4). The identification numbers of the four animals are MA026, MA027, MY144 and MY147. The four corner electrodes were excluded for all the animals. For animal MY144, electrodes [3,8] and [5,4] were excluded due to bad signal quality. The missing channels were interpolated by a Gaussian filter with filter width of 3 electrodes and standard deviation of 0.6.
Detection of gamma burst patterns. Electric line noise at 50 Hz was removed by applying a second-order Butterworth band-stop filter to the raw LFP signals. Zerophase forward and reverse filters were used to preserve features in the filtered waveform. Following Lundqvist et al. 12 , bandpass filtering and/or wavelet analysis were applied to calculate gamma-band signal amplitude; the methods yielded similar results and comparable burst extraction outcomes (Fig. 2d). For bandpass filtering analysis, a 4th order Butterworth bandpass filter was applied to the LFP signals, the first and final 1 s of the record were discarded to reduce boundary effects, and amplitude was extracted using the Hilbert transform. For wavelet analysis, the Morse wavelet transform was applied to the LFPs with symmetry parameter 3 and the time-bandwidth product 12 .
To characterize gamma burst patterns, a smoothing Gaussian filter (s.d. 0.6) was applied to the resized (factor 4, bi-cubic interpolation) electrode matrix of gamma amplitudes. Other smoothing and resizing parameters yielded comparable results. Neural data are then represented as z(x, y, t); z is the recorded signal at the position (x, y) at time moment t. Following Lundqvist et al. 12 , a burst event at (x, y) is defined where the gamma amplitude exceeds (by 2.5 s.d.) the trial mean value; other threshold values yielded similar results. We defined burst structures as a spatial connected set of burst points, with connectedness defined in terms of the six orthogonal neighbours in the Cartesian mesh of the recorded sites ("bwconncomp" routine in MATLAB 2016b; MathWorks, Natick, MA). If the burst structures were connected in the consecutive time framework, these structures were identified as gamma burst patterns. Burst patterns with duration >60 ms and spatial extent (maximum in x or y direction >1.2 mm) were analyzed. 7027 patterns were detected for the four animals (1866, 1595, 1615, and 195l patterns for the animal MA026, MA027, MY144 and MY147, respectively).

Analysis of step lengths.
To study the trajectory of the gamma burst patterns, the centre of mass X t ¼ x t ; y t À Á of the detected patterns at each time-step was calculated as the instantaneous location of the patterns, where X i t x t ; y t À Á is the location of the i-th burst within a detected gamma burst pattern at time t and N t is the total number of bursts within the pattern. To demonstrate the robustness of our results, we also tested our results by defining the amplitude weighted centre of mass (X w t ¼ Next, the time series of locations was converted into a series of movement steps. Here, two methods were used to define movement steps. The first method used is the 1D model. In this model, the 2D pattern trajectories were projected to the x direction. Turns were identified as a reversal in direction in the 1D movement. Step length was then defined as the distance between the turns. The second is the 2D turning angle model. Given two points at X t ð Þ and X t þ 1 ð Þ, we calculated the distance between these two points and considered this distance a step if the angle formed by X t ð Þ, X t þ 1 ð Þand the next point X t þ 2 ð Þis greater than the characteristic angle (θ as shown in Fig. 3h). If the condition was not met, we calculated the distance between X t ð Þ and X t þ 2 ð Þand determined whether these points met our criteria. This process was continued until a step length was found that met the criteria described above. In addition, the distance between the location where one pattern disappears and a new pattern emerges was also considered as a step. Different turning angles were used to interrogate the 2D turning angle models ( Supplementary Fig. S3).
Maximum likelihood estimation (MLE) and model selection. We used the MLE method 78 to fit truncated power law, exponential, normal, log-normal and gamma distributions. Briefly, MLE finds estimated parameter that maximizes the likelihood function in all the models. The probability density functions for all the models are shown in Table 1. The step lengths smaller than 20 µm were firstly discarded as they are below the minimal resolution. The remaining step lengths were used to find the best fit for all the models using MLE respectively. When fitting the truncated power law, an iterative process was used to derive a best fit value for the parameter a and b (the lower bound and the upper bound of the truncated power law distribution).
After the best fit for all the distributions were found, the log-likelihoods (LLH) for all the distributions were calculated 78 . Then, AIC weights were employed for model selection 3,25 , where L ð Þ is the likelihood function and K is the number of estimable parameters (the value of θ) in the approximating model (probability distribution).
As AIC values contain arbitrary constants and are greatly affected by the sample size, they do not represent an absolute metric and cannot be used directly. The following transformation makes the result an interpretable metric 3 : where AIC j is the AIC value of the jth model and AIC min is the minimum of different AIC values. The Akaike weights w j are useful as the weight of evidence 3 , where R is the size of a set of the approximating models (distribution). In addition to wAIC, a Vuong test of the truncated power law was used to check if the competing distribution should be rejected.
Neural circuit model of excitory and inhibitory neurons. Our circuit model was described in Gu et al. 15 . Briefly, N E excitatory and N I inhibitory model neurons are uniformly located in a 2D square plane and connected by chemical synapses with probability p αβ (α is the post-synaptic neuron population, and β is the pre-synaptic population); p αβ is an exponential distance-dependent connection probability function between connected neurons. The subthreshold membrane potential V α i of neuron i in population α follows where the membrane capacitance C = 0.25 nS, the leak conductance g L = 16.7 nS, and V L = −70 mV is the reversal potential for the leak current. I α i;K t ð Þ is the potassium current, I α i;rec t ð Þ is the recurrent synaptic current received by the neuron and I α i;ext t ð Þ is the external current. When the membrane potential reaches the threshold V th ¼ À50 mV, a spike is emitted and the membrane potential is reset to the potential V rt ¼ À60 mV for an absolute refractory period τ f ¼ 4 ms. The potassium current is given by is the active potassium conductance and V K = −85 mV. The dynamics of the potassium conductance are described by where t α j;k is the time of the k th spike emitted by neuron i from population α, 4g K = 10 nS and τ K = 80 ms. Because spike frequency adaptation has been primarily observed in cortical pyramidal neurons, we only include such adaptation for excitatory neurons in our model.
The recurrent synaptic current I α i;rec t ð Þ in Eq. (6) is: where g αβ i t ð Þ is the conductance of the recurrent current from the pre-synaptic population β. The excitatory and inhibitory reversal potential are V E rev ¼ 0 mV and V I rev ¼ À80 mV, respectively. The conductance g αβ i t ð Þ is given by where connection topology a αβ ij and strength J αβ ij capture the numbers of connections and coupling weights respectively, as detailed in Gu et al. 15 . The nondimensional gating s αβ ij t ð Þ describes the synaptic dynamics, where τ β d and τ β r are decay and rise time constants respectively, t β j;k is the time point of the k th spike of neuron j from population β, and d αβ ij is the conduction delay drawn from a uniform distribution between 0 and 4 ms. We implemented a vital neurophysiological feature of mammal primary cortex in our model 15 , that is, there is an identical ratio between the excitatory post-synaptic currents and inhibitory post-synaptic currents with respect to the whole excitatory population. To model this in our circuit, we consider the I-E ratio ik =∑ j J EE ij , where K EI i;in denotes the number of connections (in-degree) received by neuron i from the inhibitory population. To equalize the I-E ratio ξ i across the neurons to a desired network-wide ratio, that is, <ξ i > ¼ ξ, the J EI ik values for neuron i are sampled from a Gaussian distribution with a mean equal to ξ∑ j J EE ij =K EI i;in and a standard deviation that is 25% of the mean. The I-E ratio ξ is varied as a system parameter to explore the spatiotemporal dynamics of our cortical circuit model. Model parameters are the same as in Gu et al. 15 , except for τ I d ¼ 5 ms. N E ¼ 63 63 ¼ 3969 excitatory neurons and N I ¼ 1000 inhibitory neurons are modeled. Each neuron is described by a 1-dimensional differential equation (i.e., the leaky integrate-and-fire model). The circuit model is composed of a large number of leaky integrate-and-fire neurons (i.e., 4969); it is thus a high-dimensional dynamical system.
Neighboring excitatory neurons (one grid unit) are separated by 7.4 μm. Simulations used the forward Euler method with a time-step of 0.1 ms. Initial membrane potentials were uniformly distributed between V rt = −60 mV and V th = −50 mV. A typical trial covered 100 s, with first 500 ms excluded. Variance (SEM) was calculated from 100 trials unless otherwise stated. All these parameters were calibrated based on the biophysics properties of neurons. Near-criticality dynamics in the transition region between the asynchronous state and the propagating wave state was characterized by the branching parameter of spikes and the susceptibility 15 , two typical quantities indexing critical transitions. In the present study, we only changed one parameter, the I-E ratio, to study the spatiotemporal nature (i.e., Lévy walk characteristic) of the emergent activity patterns.
To calculate local field potential signals in our model, we adapted the LFP proxy as used by Mazzoni et al. 79 , in which LFP was approximated as the sum of the absolute values of excitatory and inhibitory currents. This proxy was able to accurately reproduce the observed LFP power spectra. Because LFP signals can be well captured by a weighted spatial sum within a two-dimensional Gaussian window centered on the recording electrode 80 , we applied a spatial Gaussian window to the sum of absolute values of currents: where r e is the coordinate vector of the simulated LFP recording electrode, r i is the coordinate vector of neuron i and I αβ i is the total current received by neuron i in population α from population β. A spatial scale σ L = 8 (56-80 μm) was used, consistent with the upper estimated limit~100 µm 80 . Our results are not sensitive to the spatial scale. The simulated LFP signal was bandpass filtered between 1 and 1000 Hz with a 4th order Butterworth filter to model the experimentally measured broadband LFP. Synaptic currents were estimated as: where I AMPA and I GABA are interior excitatory and inhibitory synaptic currents, I ext is external excitatory current, V I and V E are the inhibitory and excitatory synaptic reversal potentials, and V t ð Þ is the time-varying membrane potential.
Statistics and reproducibility. The results were expressed as mean ± standard deviation (s.d.). The experimental data were analyzed and compared against surrogate data using Wilcoxon rank sum test. Fifty surrogate datasets were generated and compared. The models were selected based on a Vuong test. The p values < 0.05 were considered to be significance.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.