Papers in Physics, vol. 7, art. 070017 (2015) Received: 18 October 2015, Accepted: 10 November 2015 Edited by: E. Mizraji Reviewed by: J. Lin, Department of Physics, Washington College, Maryland, USA. Licence: Creative Commons Attribution 3.0 DOI: http://dx.doi.org/10.4279/PIP.070017 www.papersinphysics.org ISSN 1852-4249 How we move is universal: Scaling in the average shape of human activity Dante R. Chialvo,1 Ana Maŕıa Gonzalez Torrado,2 Ewa Gudowska-Nowak,3 Jeremi K. Ochab,4 Pedro Montoya,2 Maciej A. Nowak,3, 4 Enzo Tagliazucchi5 Human motor activity is constrained by the rhythmicity of the 24 hours circadian cycle, including the usual 12-15 hours sleep-wake cycle. However, activity fluctuations also ap- pear over a wide range of temporal scales, from days to a few seconds, resulting from the concatenation of a myriad of individual smaller motor events. Furthermore, individuals present different propensity to wakefulness and thus to motor activity throughout the cir- cadian cycle. Are activity fluctuations across temporal scales intrinsically different, or is there a universal description encompassing them? Is this description also universal across individuals, considering the aforementioned variability? Here we establish the presence of universality in motor activity fluctuations based on the empirical study of a month of con- tinuous wristwatch accelerometer recordings. We study the scaling of average fluctuations across temporal scales and determine a universal law characterized by critical exponents α, τ and 1/µ. Results are highly reminiscent of the universality described for the average shape of avalanches in systems exhibiting crackling noise. Beyond its theoretical relevance, the present results can be important for developing objective markers of healthy as well as pathological human motor behavior. 1 Consejo Nacional de Investigaciones Cient́ıficas y Tec- nológicas (CONICET), Rivadavia 1917, Buenos Aires, Argentina. 2 Institut Universitari d’Investigacions en Ciències de la Salut (IUNICS) & Universitat de les Illes Balears (UIB), Palma de Mallorca, Spain. 3 M. Kac Complex Systems Research Center and M. Smoluchowski Institute of Physics, Jagiellonian Univer- sity, Kraków, Poland. 4 Biocomplexity Department, Ma lopolska Center of Biotechnology, Jagiellonian University, Kraków, Poland. 5 Institute for Medical Psychology, Christian Albrechts University, Kiel, Germany. I. Introduction The most obvious periodicity of human (as well as animal) motor activity is the circadian twenty four hours modulation. However, smaller fluctuations are evident on a wide range of temporal scales, from days to a few seconds. Data shows that the activity evolves in bursts of all sizes and durations which are known to be scale-invariant [1–8] regardless of the origins and intended consequences of such activity. Despite the variety of results, the mechanisms un- derlying the scale-invariant behavior of motor ac- tivity remain to be elucidated. Considering the in- termittent nature of human motor activity - com- prising brief activity excursions separated by peri- ods of quiescence - a natural approach would be to study the average shape of the events, following re- cent results [9–12] which show that for a large class of processes, the average shape is a scaling function 070017-1 Papers in Physics, vol. 7, art. 070017 (2015) / D. Chialvo et al. determined mostly by the temporal correlations of the process and its nonlinearities [13]. In the present work, long time series of human motor activity are analyzed, recorded via wrist- watch accelerometer, lasting approximately one month. We establish first the presence of truncated scale-invariance in the distribution of the durations of the events as well as in its power spectral den- sity, as described previously in similar type of data. Afterwards, we uncover the average shape of the bursts of activity and derive the scaling function and its associated exponents. Finally, we discuss the origins of such scaling and some possible appli- cations. II. Materials and methods The recordings analyzed were part of a larger study and included six healthy, non-smokers, drug-free volunteers (mean age 50.1 years, S.D. = 6.8). The study was approved by the Bioethics Commission of the University of Isles Baleares (Spain). Par- ticipants were informed about the procedures and goals of the study, and provided their written con- sent. After determining their handedness, each subject was provided with a wristwatch-sized activ- ity recorder (Actiwatch from Mini- Mitter Co., OR, USA) measuring acceleration changes in the fore- arm in any plane. Each data point of activity corre- sponded to the number of zero crossings in acceler- ation larger than 0.01 G (sampled at 32 Hz and in- tegrated over a 30-second window length). Records of several thousands of data points were kept in the device’s internal memory until being downloaded to a personal computer every week. Subjects wore the device in their non-dominant arm continuously for up to several weeks (mean 28.1 days, S.D.= 4.). Af- ter careful visual inspection of the data to exclude sets with gaps (due to subject non-compliance), a combined total of 280 days of data was available for further analysis. III. Results For ease of presentation, we will use recordings from a single subject to describe the main results. Nev- ertheless, results are robust as well as similar for the entire group of subjects in the study. A typical recording is presented in Fig. 1. Panel A shows 6 8 10 12 14 16 18 20 22 0 2 4 6 Mon12/14 Tue12/15 Wed12/16 Thu12/17 Fri12/18 Sat12/19 Sun12/20 Mon12/21 Tue12/22 Wed12/23 Thu12/24 Fri12/25 Sat12/26 Sun12/27 Mon12/28 Tue12/29 Wed12/30 Thu12/31 Fri01/01 Sat01/02 Sun01/03 Mon01/04 Tue01/05 Wed01/06 Thu01/07 Fri01/08 Sat01/09 Sun01/10 Mon01/11 Tue01/12 Hour D at e 6 8 10 12 14 16 18 20 22 0 2 4 6 0 0.5 1 A ct iv ity (a .u ) 0 7 14 21 28 −5 0 5 Time (days) I ( n ) / S .D . −5 0 5 −4 −3 −2 −1 0 I /S.D. lo g 1 0 P ( I ) / P ( 0 ) −6 −5 −4 −3 −2 −6 −4 −2 0 log10f (sec − 1) lo g 1 0 S ( f ) B C D A 1 day 1 hour 1 min. Figure 1: Example data set, distribution of suc- cessive increments and their spectral power. Panel A: Time series of activity x(n) recorded continu- ously from a subject during a month. Individual traces correspond to consecutive days. The top subpanel depicts daily activity averaged over the entire month. Panel B: Time series of successive increments I(n) = x(n + 1) −x(n) (normalized by its SD) for the same data. Panel C: Probability density distribution of the time series of succes- sive increments I(n) (continuous line), exhibiting exponential tails (compare with the dotted line, a Gaussian of the same variance). Panel D: Power spectral density (black line) of the time series of successive increments I(n) of panel B. This is scale invariant S(f) ∼ fγ with γ = 0.9 (dashed line). In contrast, for the randomly shuffled increments, the serial correlations vanish and a flat spectral density is obtained (red). a full month of continuously recorded activity from this subject, who is particularly regular in her daily routines. The subject wakes up with the alarm clock at 6:45 a.m. on week days and has lunch fol- lowed by a short nap each day (between 2:00 p.m. and 4:00 p.m. Panel B displays the time series of 070017-2 Papers in Physics, vol. 7, art. 070017 (2015) / D. Chialvo et al. 102 103 104 T 102 103 104 105 100 101 102 103 104 105 S,Wt,T 10-3 10-2 10-1 100 C C D F S Wt T A B α' τ' µ + 1 Figure 2: Scaling of activity events in a single subject (same dataset as in Fig. 1). Panel A: The complementary cumulative distribution func- tion (CCDF) for event durations (T) and sizes (S) obeys power-laws with exponents α′ = 0.70 and τ′ = 0.44, respectively (dashed lines). Note that here the densities are cumulative, thus the expo- nents of the respective PDFs are α = α′ + 1 and τ = τ′ + 1. The waiting time between events falls exponentially. Panel B: The average size of a given duration is well described (for small T) by 〈S〉(T) ∼ Tµ+1 with µ + 1 = 1.59 (blue dashed line) comparable with results obtained from fitting within the scaling region (red filled symbols) giving µ + 1 = 1.61. the successive increments of the signal x(n), defined as I(n) = x(n + 1) −x(n). The large-scale statistical features of the time se- ries presented in Fig. 1 are already well known. The density distribution of the successive incre- ments i(n) is non-Gaussian, as can be appreci- ated by a joint plot with a Gaussian distribution of the same variance (Fig. 1, Panel C). It is known that the power spectrum of the activity decays as S(f) ∼ fβ [1, 2]. Because this type of processes are likely to be non-stationary, it is best to esti- mate the exponents of the spectral density by doing the calculations over the time series of successive increments, whose density distribution is station- ary. For instance, for Brownian motion (which is summed white noise), the power spectrum decays S(f) ∼ fβ with β = −2 and for white noise β = 0; the summed time series has an exponent +2 larger than the non-summed time series. As discussed in [14], this can be generalized for all self-affine pro- cesses: summing a self-affine time series shifts the theoretical power-spectral density exponent by +2, and the reverse process is also true: the differences in consecutive values (the “first differences”) of a Brownian motion result in white noise, thus tak- ing the first differences shifts the theoretical power- spectral density exponent β by −2. In our case, the exponent obtained for the time series of successive increments I(n) was γ = 0.9. Thus, the exponent of the raw data is β = γ − 2 = −1.1 [14]. For com- parison, the spectral densities of the actual signal and of a surrogate obtained after randomly shuf- fling the increments are jointly displayed in Panel D of Fig. 1. To further study the time series from the perspec- tive of individual bursts of activity, we introduce the definition of an event. We consider the time series of activity x(n) and select a threshold value U to be vanishingly small. An event is defined by the consecutive points starting when x(n) > U and ending when x(n) < U. This is equivalent to the definition of avalanches in other contexts [9, 15]. In the following part, we will be concerned primarily with the statistics of event lifetimes T, as well as of their average size S and shape. In all subjects, we found that the distributions of event durations and sizes (defined by the area, i.e., the integral of the signal corresponding to the individual events) can be well described, for relatively small values, by a power-law (Fig. 2, Panel A). In contrast, the dis- tribution of waiting times between events demon- strated an exponential decay. In addition to the scale invariance, we found that the longer an event lasted, the stronger the motor activity executed by the subject. The plot of average event size 〈S〉 as a function of duration T follows a power-law (for small values of T) described by 〈S〉(T) = Tµ+1 with µ + 1 = 1.59. The exponents in this power-law are 070017-3 Papers in Physics, vol. 7, art. 070017 (2015) / D. Chialvo et al. 0 0.2 0.4 0.6 0.8 1 t/T 0 5 10 15 / T µ 0.2 0.4 0.6 µ 0.1 0.15 0.2 ε 0 0.2 0.4 0.6 0.8 1 t/T 0 0.2 0.4 0.6 0.8 1 < x> / T µ / max (4(t/T)(1-t/T))0.48 (4(t/T)(1-(t/T))0.59 0 600 1200 1800 Time(sec) 0 150 300 A ct iv ity < x( t) > T A B C Figure 3: Collapse of events of different duration into a single functional form. Panel A: Three ex- amples of typical events of duration T=480, 960 and 1920 sec.. Panel B: The heterogeneous events shown in Panel A can be collapsed onto the av- erage shape (dashed black line) by normalizing t to t/T and 〈x(t)〉 to 〈x(t)〉/Tµ. The inset shows the cumulative variance for a range of µ. Panel C: The average event shape, i.e., fshape(t/T), re- covered from six data sets (thin lines). The best fit using an inverted parabola is shown as a red dashed line (µ = 0.49) as well as the one expected from the critical exponent µ = 0.59 as a dot-dashed blue line. robust across subjects and to changes of threshold over a reasonable range of values. This type of scaling is well known in the statisti- cal mechanics of critical phenomena [15]. Examples range from earthquakes [16] to active transport pro- cesses in cells [17], crackling noise [11], the statis- tics of Barkhausen noise in permalloy thin films [10] and plastic deformation of metals [18]. In all these cases, the distributions obey universal functional forms: f(S) ∼ S−τ, (1) f(T) ∼ T−α, (2) 〈S〉(T) ∼ T1/σνz, (3) where f denotes the probability density functions of the size of the event S and its duration T, and 〈S〉(T) is the expected size for a given duration. The parameters τ, α and 1/σνz are the critical ex- ponents of the system and are expected to be inde- pendent of the details, being related to each other by the scaling relation: α− 1 τ − 1 = 1 σνz . (4) We found that the empirical exponents very closely fulfill the expression above. Using the fitting approach introduced by Clauset [19] in the scaling regions depicted in Panel A of Fig. 2, we found τ = 1.44 and α = 1.70 . Thus, from Eq. (4) a value of 1/σνz = µ + 1 = 1.59 is expected. The experimental data points are very close to this the- oretical expectation (dashed line), especially for the relatively small T values within the scaling region of Panel A (where a linear fit estimates µ + 1 = 1.61), while those for relatively larger T values (corre- sponding to the cutoff of the distributions) are a bit apart, probably due to undersampling. After repeating this analysis for all subjects in our sam- ple, the average exponents were all within 5% of the reported values. From scaling arguments, it is expected that the average shape of an event of duration T 〈x(T,t)〉 scales as : 〈x(T,t)〉 = Tµfshape(t/T). (5) Thus, the shapes of events of different durations T rescaled by µ should collapse on a single scaling function given by fshape(t/T). Note that µ corre- sponds in this context to the wandering exponent (i.e., the mean squared displacement) of the activ- ity [13, 20]. Examples of this collapse are presented in Pan- els A and B of Fig. 3. Considering the number of 070017-4 Papers in Physics, vol. 7, art. 070017 (2015) / D. Chialvo et al. 102 103 T 102 103 104 0 2000 4000 S,Wt,T 10-4 10-3 10-2 10-1 100 C C D F S Wt T A B Figure 4: Scaling is absent in a null model result- ing from defining events after randomly reordering the time series x(n). Panel A: Density distribu- tions (CCDF) for event duration, size, and waiting time. All the distributions are exponential (note the logarithmic-linear scale). Panel B: The ex- pected average size for a given duration in the null model is a linear function of T (the dashed line rep- resents the fit with slope 1), therefore, µ = 0 and there is no collapse. events here averaged (in the order of N ∼ 102), the data collapse is quite satisfactory, while the value of the exponent (µ = 0.48) does not exactly match the one predicted in Eq. (4), µ = 0.59 (likely a conse- quence of insufficient sampling). To determine the generality of our results, we extended this analy- sis to six other data sets. For each data set, the value of µ was first determined. Subsequently, the x(T,t) obtained from the events were rescaled with Tµ and their average computed. To account for individual differences in mean activity, shape func- tions were normalized by their mean value. The results for the six datasets are presented in Panel C of Fig. 3. They can be accurately described by an inverted parabola, as in other systems previously studied using this method. The best fit disagrees with the empirical functions near their peak, the latter being flatter, likely an effect related to satu- ration observed in long events. Finally, we turn to discuss simple null models. We consider two extreme cases, in both of them the raw time series are randomly shuffled to remove serial correlations. In the first case, we remove all temporal correlations by randomly reordering x(n), thus attaining a flat power spectral density. After repeating the above analysis in this surrogate data set, it becomes clear (as shown in Fig. 4) that the scale invariance is absent in all the statistics un- der study: size S, waiting time Wt and duration T of events (note that the distributions are here plotted using a logarithmic-linear scale). Results in Panel B show that µ + 1 = 1, thus µ = 0, imply- ing that there is not collapse, because with Tµ = 1 in Eq. (5), the amplitude of the individual events remains invariable. To consider the second case, we need first to reorder randomly the time series of increments I(n) and then proceed to integrate the increments. Since each increment is now a random variable, the power spectral density for this surro- gate process obeys fβ with β = −2 , and as shown analytically by Baldassari et al. [13], for this case µ = 1/2 and the scaling function is a semicircle. Please note that the fluctuations of human activity described here differ from a simple auto-regressive process: indeed successive increments I(n) are anti- correlated and the power spectral density corre- sponds to non-trivial power law correlations (i.e., β 6= −2). IV. Discussion The present findings can be summarized by six styl- ized facts describing bursts of human activity: I) the spectral density of the time series of activity x(t) obeys a power law, with exponent β ∼ 1; II) successive increments I(n) are anti-correlated with a spectral density obeying a power law with exponent γ ∼ 1, which corresponds to a spec- tral density for the raw data fβ with β ∼ −1; III) the PDF of the increments I(n) is definitely non-gaussian; IV) the PDF of duration and sizes of events obeys truncated power laws with expo- nents 1 < τ < 2 and 1 < α < 2; V) the aver- 070017-5 Papers in Physics, vol. 7, art. 070017 (2015) / D. Chialvo et al. age size of the events scales with its lifetime T as 〈S〉(T) ∼ Tµ, where µ + 1 = (α − 1)/(τ − 1); VI) the time series of individual events can be appropri- ately rescaled via a transformation of its duration T and amplitude x(t) onto a unique functional shape: 〈x(T,t)〉 = Tµfshape(t/T). We are aware that these observations are novel only for human activity, because similar statistical regularities of avalanching activity are well known for a large variety of inanimate systems [9–12]. The rescaling of the average shape is not surprising be- cause, placed in the appropriate context, it can be traced back to Mandelbrot’s study of the fractal properties of self-affine functions [21]. A curve or a time series are said to be self-affine if a transfor- mation can be found, such that rescaling their x,y coordinates by k and kµ, respectively, and the vari- ance in y is preserved (with µ = 1 corresponding to self-similarity). In that sense, the successful col- lapse of the events shape is a trivial consequence of the overall self-affinity of the x(t) time series. Thus, it is clear that the existence of the scal- ing uncovered here is not informative per se of the type of mechanism behind: scale-invariance can be constructed via different processes, ranging from critical phenomena [15] to simple stochastic auto- regressive dynamics [13, 20]. What is then the mechanism by which the above six facts are gen- erated? It seems that this question cannot be easily an- swered by the type of experiments reported here. Fluctuations of this type could have either an in- trinsic (i.e., brain-born) origin but also could be the reflection of a collective phenomena (including humans and its environment). In either case, the correlations observed seem to reject the case of in- dependent random events starting and stopping hu- man actions, because neither the distribution of the increments I(n), nor the exponents match the case of a random walk. In terms of brain-born process, it is hard to accept some of the implications of the scaling function in the activity shape. The average parabolic shape means that the very beginning of the motion activity contains information about how long the activity will last, in the same sense that the initial trajectory of a projectile predicts when and where it will land. This proposal is hardly realistic, because there is hardly a reasonable physiological argument in support of any motor planning for the length of time we are observing (∼ 103 secs). In terms of collective processes, the results here sug- gest that the interaction with other humans could determine when and where, on the average, we start and stop moving. Despite our current relative ignorance, a possi- bility that sounds interesting is to determine in children, as they grow, if their behavioral product of parental (and otherwise) education are reflected in the shape of their individual scaling function. This seems reasonable given the fact that “tireless running around” is almost a definition of early age well-being, which gives way to less hectic activity as children mature. In the same line of thoughts, if changes in the scaling function can be quantita- tively traced to behavioral changes, one could also consider to explore applications of these techniques to monitor eventual progress in the treatment of hyperactivity disturbances such as in the subjects affected by the Attention Deficit Hyperactivity Dis- order syndrome. The converse, i.e., cases in which the average activity diminish, as in elderly subjects shall be also explored. Further experiments and analysis should shed light on these possibilities. In, the meantime, the present results provide a guide and six important constraints for the models that should best capture the physics (and biology) of the process. Acknowledgements - Work supported by Na- tional Science Center of Poland (ncn.gov.pl, grant DEC-2011/02/A/ST1/00119); State Secretary for Research and Development (grants PSI2010-19372 and PSI2013-48260) from Spain and by CONICET from Argentina.) [1] T Nakamura, K Kiyono, K Yoshiuchi, R Naka- hara, Z R Struzik, Y Yamamoto, Universal scaling law in human behavioral organization, Phys. Rev. Lett. 99, 138103 (2007). [2] T Nakamura, et al., Of mice and men - uni- versality and breakdown of behavioral organi- zation, PLoS ONE 3, e2050 (2008). [3] K Hu, P C Ivanov, Z Chen, M F Hilton, H E Stanley, S A Shea, Non-random fluctuations and multi-scale dynamics regulation of human activity, Physica A 337, 307 (2004). 070017-6 Papers in Physics, vol. 7, art. 070017 (2015) / D. Chialvo et al. [4] L A N Amaral, D J B Soares, L R da Silva, L S Lucena, M Saito, H Kumano, N Aoyagi, Y Ya- mamoto, Power law temporal auto-correlations in day-long records of human physical activ- ity and their alteration with disease, Europhys. Lett. 66, 448 (2004). [5] C Anteneodo, D R Chialvo, Unravelling the fluctuations of animal motor activity, Chaos 19, 033123 (2009). [6] K Christensen, D Papavassiliou, A de Figueiredo, N R Franks, A B Sendova-Franks, Universality in ant behaviour, J. R. Soc. Inter- face 12, 20140985 (2014). [7] A Proekt, J Banavar, A Maritan, D Pfaff, Scale invariance in the dynamics of sponta- neous behavior, Proc Natl Acad Sci USA 109, 10564 (2012). [8] J K Ochab, et al., Scale-free fluctuations in be- havioral performance: Delineating changes in spontaneous behavior of humans with induced sleep deficiency, PLoS ONE 9, e107542 (2014). [9] L Laurson, X Illa, S Santucci, K T Tallak- stad, K J Maloy, M J Alava, Evolution of the average avalanche shape with the universality class, Nature Comm. 4, 2927 (2013). [10] S Papanikolaou, F Bohn, R L Sommer, G Durin, S Zapperi, J P Sethna, Universality beyond power laws and the average avalanche shape, Nature Phys. 7, 316 (2011). [11] J P Sethna, K A Dahmen, C R Myers, Crack- ling noise, Nature 410, 242 (2001). [12] N Friedman, S Ito, B A W Brinkman, M Shi- mono, R E L DeVille, K A Dahmen, J M Beggs, T C Butler, Universal critical dynam- ics in high resolution neuronal avalanche data, Phys. Rev. Lett. 108, 208102 (2012). [13] A Baldassarri, F Colaiori, C Castellano, Av- erage shape of a fluctuation: Universality in excursions of stochastic processes, Phys. Rev. Lett. 90, 060601 (2003). [14] B D Malamud, D L Turcotte, Self-affine time series: I. Generation and analyses, Adv. Geo- phys. 40, 1 (1999). [15] P Bak. How nature works. The science of self- organized criticality, Copernicus, New York (1996). [16] B Gutenberg, C F Richter, Magnitude and en- ergy of earthquakes, Ann. Geofis. 9, 1 (1956). [17] B Wang, J Kuo, S Granick, Burst of active transport in living cells, Phys. Rev. Lett. 111, 208102 (2013). [18] L Laurson, M J Alava, 1/f noise and avalanche scaling in plastic deformation, Phys. Rev. E 74, 066106 (2006). [19] A Clauset, C R Shalizi, M E J. Newman, Power-law distributions in empirical data, SIAM Rev. 51, 661 (2009). [20] F Colaiori, A Baldassarri, C Castellano, Aver- age trajectory of returning walks, Phys. Rev. E 69, 041105 (2004). [21] B B Mandelbrot, Self-affine fractals and frac- tal dimension, Physica Scripta 32, 257 (1985). 070017-7