extfieldsofstudy: Biology Medicine

text:
Scaling-Laws of Human Broadcast Communication Enable Distinction between Human, Corporate and Robot Twitter Users

Human behaviour is highly individual by nature, yet statistical structures are emerging which seem to govern the actions of human beings collectively. Here we search for universal statistical laws dictating the timing of human actions in communication decisions. We focus on the distribution of the time interval between messages in human broadcast communication, as documented in Twitter, and study a collection of over 160,000 tweets for three user categories: personal (controlled by one person), managed (typically PR agency controlled) and bot-controlled (automated system). To test our hypothesis, we investigate whether it is possible to differentiate between user types based on tweet timing behaviour, independently of the content in messages. For this purpose, we developed a system to process a large amount of tweets for reality mining and implemented two simple probabilistic inference algorithms: 1. a naive Bayes classifier, which distinguishes between two and three account categories with classification performance of 84.6% and 75.8%, respectively and 2. a prediction algorithm to estimate the time of a user's next tweet with an . Our results show that we can reliably distinguish between the three user categories as well as predict the distribution of a user's inter-message time with reasonable accuracy. More importantly, we identify a characteristic power-law decrease in the tail of inter-message time distribution by human users which is different from that obtained for managed and automated accounts. This result is evidence of a universal law that permeates the timing of human decisions in broadcast communication and extends the findings of several previous studies of peer-to-peer communication.

Introduction
The dynamics of many social, technological and economic phenomena are driven by individual human actions, therefore the quantitative understanding of human behaviour is becoming a central question in neuroscience, economics and social sciences. Human behaviour is highly variable from trial to trial [1], yet appears highly stereotyped (e.g. we can recognise a mime's actions). In neuroscience, the variability in human decision making and motor behaviour has been found to display a characteristic variability structure [2] than can be used to predict the behavioural decisions and performance of individuals [3]. We are investigating here to what extent these computational neuroscience approaches can be applied to analyse human communication decisions on the online social network Twitter, specifically to understand the timing of tweeting. We follow a very simple, easily interpretable approach using non-parametric Bayesian statistics to analyse and then predict the nature of the tweeter, i.e., is the tweeter a genuine individual or somebody or something else. We focus here on the statistical structure of broadcast communication by employing a large anotated dataset collected from Twitter, with the aim of better understanding the timing of human actions in this type of communication and how individual or different they are from each other. While past research has largely focused on using tweets as a representation of collective behaviour [4][5][6], our individual-based approach takes a neuroscience perspective of reality mining [7,8] and uses Twitter data to study users individually and make predictions about them in real life.
Since its creation in 2006, Twitter has become an increasingly popular medium enabling over 500 million active users (Summer 2012) who produce 65 million tweets per day. The popularity of Twitter makes it an important tool for journalism, marketing, political campaigns and social change. It is therefore of immediate interest to be able to determine if the user generating the tweets (irrespective of tweeted content) is 1. a genuine individual, 2. a group of people appearing as one Twitter user (e.g. a corporation or celebrity having a dedicated PR team handling their 'personal' tweets) or 3. an automated system ('bot') that generates tweets. We approach this by creating a non-parametric naive Bayes classifier based on tweeting time. This classification can be very helpful in the recognition and filtering of spammers and malicious accounts, and can therefore assist in understanding the online community and help us recognise who is actually tweeting. The first step in the development of our study was the collection of data from multiple Twitter users through a web crawler. For this purpose, we created the Twitter Reality Miner application, a Twitter crawler which allowed us to retrieve data in an efficient way while conforming to the request limit imposed by the Twitter API. After data collection, we studied tweeting patterns and the probability distributions of timestamps and time intervals in between posts. We then used this information to classify user accounts into three different groups (personal, managed and botcontrolled) and to predict the probability distribution of the time delay before the next tweet of a user was posted.

Related Work
With the growing popularity of online social networks and other means of interaction, recent research has taken advantage of the large amounts of freely available digital data in order to investigate several aspects of human behaviour. The novel field of reality mining, for instance, applies machine-sensed environmental data to the study of human activities in real life. Eagle et al. have used data from 100 mobile phone users in an American educational institution including call logs, Bluetooth devices in proximity, cell tower IDs and phone status, and found that mobile phone usage consistently correlates with users' activities. The authors apply this data to accurately predict real-life friendships and individual-level measures such as job satisfaction [7]. In a different study, the authors extract the mobile phone dataset's principal components and use them to predict user activity on the same day with 79% accuracy [8].
In addition, related studies of human communication behaviour have studied modern e-mail communications and web browsing, as well as Einstein and Darwin's documented correspondence patterns. All these studies consistently find that human communication intervals are governed by an underlying statistical structure, which largely dictates how and when these actions are performed, regardless of the individual characteristics of each person [9]. Barabási and colleagues studied e-mail communication patterns in order to understand how humans prioritise their activities and proposed a priority model which predicts that interevent times should display a heavy-tailed distribution, as found in power-laws [10]. These distributions arise from individuals displaying long periods of inactivity which alternate with bursts of intense activity, a behaviour characteristic of the timing of many human actions, from communication to entertainment and work patterns [11][12][13]. Later, this work was extended to studying Darwin's and Einstein's patterns of correspondence and comparing them with today's e-mail exchanges [14]. The authors found that the probability that a letter would be replied to in t days is well approximated by a power-law, thus following the same scaling laws as current e-mail communication. Dezsö et al. [15] investigated the topology and features of dynamically changing human interaction networks by analysing the visitation patterns of a major news portal. They showed that the timing of the browsing process is not the commonly assumed Poisson process, but instead suggests that heavy tails are a part of a universal scaling law, representing a fundamental pattern of human decision making dynamics.
Here we look at broadcast communication, an aspect of human interaction that has not been studied before in this context. We test and apply our analysis by focusing on identification and classification of specific types of users on Twitter. This classification can be useful for a variety of reasons, from focusing advertisement and political campaigns, to filtering spam and malicious accounts. With a large occurrence of spamming and political campaigning on Twitter, recent research has focused on methods for identifying certain types of behaviour that are characteristic of spammers or propagandists. In [16], Lumezanu et al. aim to understand how Twitter is used to spread propaganda. They studied the Twitter behaviour of propagandists, users who consistently express the same opinion or ideology, and focused their work on hyperadvocates, who show a consistent lack of impartiality in their messages. Four publishing patterns were found to amplify the effect on hyperadvocacy on Twitter. Another example of Twitter account classification can be found in [17]: Chu et al. observe the differences between Twitter accounts controlled by humans, bots, and cyborgs, which refer to either botassisted humans or to human-assisted bots. The authors studied tweeting behaviour, tweet content and account properties in order to characterise the automation feature of Twitter accounts, then used this information to build a classifier for the three account categories. Despite attaining a high correctness rate, the system created has the limitation of heavily relying on processing the contents of tweets in order to identify them as spam, which can be an expensive and time-consuming process.

Data Collection
The first step in the development of our study was the collection of data from multiple Twitter accounts through a web crawler. For this purpose, we created the Twitter Reality Miner, a Twitter crawler application which allowed us to retrieve Twitter data in an efficient way while conforming to the request limit imposed by the Twitter API. After data collection, we studied tweeting patterns and the probability distributions of timestamps and time intervals in between posts. We then used this information to classify user accounts into the three different groups (personal, managed and bot-controlled) and to predict the probability distribution of the time delay before the next tweet of a user is posted.
The application was developed in Python script language with the aid of a variety of third-party libraries (python-twitter, oauth2, httplib, json, psycopg and pyparsing), the details of which are omitted from this paper. The whole TRM application consists of four Python modules: crawler, rateLimiter, databaseAccess and errorReport. The full source code and the data collected (in text and spreadsheet files) can be obtained from the following GitHub repository: www.FaisalLab.com/TRM. Access to Twitter was possible due to the Twitter Application Programming Interface (API), a specification that allows communication between the crawler and Twitter itself. One significant shortcoming in using this API for data retrieval is its restrictive limit policy, which only allows clients to make 150 requests per hour. Even if a client makes calls to the API within the allowed limit, Twitter may throttle the account when too many calls are made repeatedly. For this reason, we created a wrapper module for the API, called rateLimiter, in order to add small time intervals in between requests, thus preventing the account from being ''black listed'' by Twitter. During data collection, the crawler was given a list of screen names of manually selected user accounts to process. Both tweets and retweets in the timelines were collected, up to a total of 800 posts per user account.

Classifying Tweeters
We have developed two classification algorithms with similar implementation: the 2-Classifier distinguishes between personal and managed accounts, while the 3-Classifier distinguishes between personal, managed and bot-controlled accounts. Both our classification systems are based on the maximum a posteriori (MAP) decision rule: where c is a specific class and f is a feature value for a particular sample [18]. According to this rule, a test sample is assigned to the class in which its features yield the largest probability value. Four attempts of classification were made, applying different probability distributions: 1. using the inter-tweet delay marginal distribution (ITD); 2. using the tweet time marginal distribution (TT); 3. using the joint distribution of the two variables assuming independence (JI); and 4. using the joint distribution of the two variables not assuming independence (JNI). We began by applying leave-one-out cross validation to our dataset. In each crossvalidation loop, N{1 sample accounts were grouped into their respective classes, then the probability density function for each class was computed. To classify the left-out account, the feature values (inter-tweet delay, tweet time, or both) of each one of that account's tweets were interpolated into the distribution of each class.
The classification score of a given class for a given account was then computed as the sum of the logarithm of the probabilities obtained for all the sample tweets of that account, when interpolated into the class distribution. For each of the four attempts, the classification score S c (i) of class c for sample i was computed as: where T(i) is the set of tweets for sample account i, and interpolate(t,pdf(c)) is the spline interpolation of the value of t into the probability density function of class c. Once all class scores had been computed, the user was classified into the class with the highest score. Since scores were computed separately for each classification attempt, a different outcome was obtained for each attempt, resulting in four different classification outcomes for each user account. During the cross-validation phase, our best results were obtained when using the joint distribution of the inter-tweet delay and tweet time variables assuming independence (JI), as shown in the Results section below. We have therefore applied this classification system in our next stage: splitting the data into separate training and test sets, using the training set to generate the probability distributions, then classifying the test samples by interpolating their values into the generated distributions [19]. To test the robustness of the algorithms, we varied the size of the training dataset between 5% and 70% of the user accounts (while the remaining accounts were using for testing). In each of these set ups, we repeated the experiment 10 times, each time reshuffling the samples among each class.

Predictive Model for Tweet Time Distribution
Our next step was to create a probabilistic model to predict when a user's next tweet would be posted, based on the inter-tweet delay distribution of that user's account class. Again, we started by applying leave-one-out cross validation to our dataset, which comprised 67 accounts from each class, resulting in 201 accounts in total. At each of the N iterations, N{1 accounts were used to generate the model, while the left-out account was used to validate the model. This allowed us to maximise the number of samples used in the model generation.
In our first predictive model, we used the inter-tweet delay distribution of each class in order to generate a corresponding cumulative distribution function (CDF). The CDF of the intertweet delay t describes the probability that a tweet will occur given that t seconds have passed since the last tweet. The actual (observed) inter-tweet delay of the tweet we wanted to predict (among the left-out sample's tweets) was then used to compute a step function as follows: where t is the actual inter-tweet delay of the left-out sample's tweet, which we aimed to predict. This step function represents the observed cumulative probability of a tweet occurring t seconds after the previous tweet: because the tweet occurred exactly after t seconds, this probability is 0 before t, and 1 after t. For each tweet of the sample user account, a different step function was computed. In order to evaluate the predictive model, each step function was compared to the class CDF using the coefficient of determination R 2 . The R 2 between each step function (observed data) and the class CDF (predictive model) was calculated as 1{SS err =SS tot , where SS err is the sum of squares of residuals, and SS err is the total sum of squares.
As an illustrative example, the prediction for 5 sample tweets in the personal accounts class is demonstrated in Figure 1. Figure 1 (a) shows the CDF computed for the personal accounts class using N{1 accounts (in red), as well as the step functions computed for 5 tweets of the left-out account (in blue). In order to evaluate how well the CDF fits each step function, we show in Figure 1 (b) a 3dimensional histogram where the axis on the left of the plane corresponds to the value of the CDF obtained for the inter-tweet delay (predicted value), and the axis on the right corresponds to the value of the step function obtained for the same delay (actual value, which is either 0 or 1). A perfect predictive model would have all data points grouped in bins f0,0g and f1,1g, indicating that the CDF models the step functions exactly and thus all predicted and actual values coincide. The fact that these bins have much higher probabilities than all others in the histogram illustrates the model's accuracy.
In addition to cross-validation, we also tested our singledistribution predictive model using separate training and test datasets. We varied the sizes of the training and test sets, starting with 30% and 70% of the samples, respectively, then increasing the training set by 10% in each experiment, until we had 70% of samples for training and 30% for testing. In each of these set ups, we repeated the experiment 10 times, each time reshuffling the samples among each class. The results of these experiments are presented in the next section.
In a slightly more elaborated version of the predictor, we used the same predictive model but with separate inter-tweet delay distributions for each hour of the day. Each inter-tweet delay data point was associated with an hour of the day based on the timestamp of the tweet that occurred before that delay. This resulted in in 24 different probability distributions for the intertweet delay, one for each hour of the day. After computing the 24 distributions, we selected which distribution to use according to the timestamp of the sample user's last tweet. Although they do not follow a standard model for prediction, both our models are based on simple probability and statistics principles [18].

Tweeting Activity Analysis
We now present the statistical analysis of the dataset retrieved through our Twitter crawler application. This dataset contains 100 manually identified and verified Twitter accounts for each of the three account classes, namely ''personal'', ''managed'', and ''botcontrolled'', and was used for analysing and comparing the behaviour of users in each account class. All managed accounts selected are maintained by large, well-known corporations, and the bot-controlled accounts were chosen based on online lists of Twitter bots. Apart from manual selection, the collected data was not filtered in any way. Table 1 shows the average, minimum and maximum number of days that accounts were active for each class. We present an analysis of the periodicity of tweets in Figure 2, which contains the periodogram power spectral density estimation of tweeting activity for each account class. No relevant predominant frequencies were found in this analysis.
The two main properties of the data studied in this paper were the tweet time (hour of the day in the respective timezone and day of the week) and the inter-tweet delay, i.e., the amount of time elapsed between two consecutive tweets by the same user. The timestamps of tweets were adjusted to the timezone of each user and users who did not specify their timezone were hence discarded from this analysis. Consequently, our dataset was reduced to 86 personal accounts, 91 managed accounts and 67 bot-controlled accounts, and we used a total of 51,924 tweets from personal accounts, 67,436 tweets from managed accounts and 45,615 tweets from bot-controlled accounts.
We begin by studying the inter-tweet delay distributions in each class. Figure 3 shows, for each class, a scatter plot of individual inter-tweet delay standard deviations vs. inter-tweet delay means (black line denotes linear proportionality). The linear fits show that the variability of inter-tweet delay is closely proportional to mean inter-tweet delay, i.e. inter-tweet delays exhibit signal-dependent noise characteristics. Figure 4 (a) shows the probability density function (PDF) for the inter-tweet delay in each class, while Figure 4 (b) shows the complementary cumulative distribution function (CCDF) for each class (blue curve -personal; pink curvemanaged; orange curve -bot), as well as the power-laws fitted to the tail of each class distribution.
The power-law decrease in the tail we have observed for this instance of broadcast communication is in accordance with results previously obtained for peer-to-peer communication, such as interletter, inter-email and inter-webpage delay [10,14,15]. To correctly fit the power-laws, we adopted maximum likelihood estimators and a goodness-of-fit approach for estimating the lower cutoff of the power-laws [20]. For the personal accounts intertweet delay distribution we obtained a slope of -2.38, from which we conclude that the tail of this distribution is well approximated by a power-law P(t)!t {a , where a~2:38+0:059. For managed accounts, typically controlled by more than one person, we obtained a~2:88+0:16, and for bot-controlled accounts we obtained a~2:71+0:18. The detailed statistics of the power-law fit for each account class are shown in Table 2. In order to verify that these distributions were not generated by the same model, we performed the two-sided Kolmogorov-Smirnov test between each pair of classes, which rejected the null hypothesis at the 5% significance level in each pair. Thus, the inter-tweet delay distributions were statistically significantly different.
We analysed the time of day tweet statistics for each user in each class using circular statistics and fitted a von Mises distribution to each account. To characterise tweet time variability around the mean, we converted the concentration parameter k of the von Mises distribution into a dispersion measure (s~ffi ffiffi 1 k r ), which is unit equivalent to standard deviations for the Gaussian distribution. Figure 5 shows, for each class, a polar plot of tweet hour of the day means (in the accounts local time zone) against individual tweet time variability. The CDF corresponds to the probability that a tweet will be posted t seconds after the previous tweet (predicted probability), while the step functions correspond to the observed probability for the occurrence of tweets (observed or actual probability). A perfect prediction for a specific tweet would mean that the CDF coincides exactly with the step function for that tweet. (b) In this histogram, the axis on the left of the plane corresponds to the value of the CDF obtained for the inter-tweet delay (predicted value), while the axis on the right corresponds to the value of the step function obtained for the same delay (actual value, which is either 0 or 1). A perfect predictive model would have all data points grouped in bins f0,0g and f1,1g, indicating that the CDF models the step functions exactly and thus all predicted and actual values coincide. The fact that these two bins have much higher probabilities than all others in the histogram illustrates the model's accuracy. doi:10.1371/journal.pone.0065774.g001 Figure 6 shows the pooled empirical PDFs for the hour of the day for all tweets in each class. We can observe that personal accounts increase their tweeting activity level as the day progresses, peaking at 9pm. Managed accounts tend to tweet more during work hours, between 9am and 6pm. The dip in the distribution at 12pm can probably be explained by lunch hour breaks. Finally, the distribution for bot-controlled accounts exhibits a variety of peaks, which is probably because their behaviour is not associated with a structured daily routine.
The weekly tweeting patterns for the 65 most active users from each class are shown in Figure 7, where each tile is associated with a user and a day of the week, and the tile's colour intensity is proportional to the amount of tweets posted by that user on that day. Managed accounts have higher tweeting activities during work days, while personal accounts present a homogenous behaviour throughout the week. The activity for most botcontrolled accounts shows little correlation with the days of the week. We do not distinguish between week days and weekends in our analysis since we are interested in the global timing behaviour of each user, regardless of the day of the week. The hourly tweeting patterns for the same users are shown in Figure 8, where each tile is associated with a user and an hour of the day, and the tile's colour intensity is proportional to the amount of tweets posted by that user at that hour. In this figure we can clearly observe the differences in behaviour between the three classes: personal accounts tend to tweet more in the afternoons and evenings; managed accounts tweet more during work hours; and bot-controlled accounts either have a regular behaviour, tweeting at an approximately constant rate throughout the day, or display a low tweet rate with a very high peak at one or a few specific hours. These behavioural plots show that the tweeting patterns for both personal and managed accounts are intrinsically related to a real life daily routine, whereas bot-controlled accounts exhibit an artificially designed behaviour. The very distinct patterns obtained for the three account classes allowed us to use tweeting behaviour as a classification criterion.

Automatic Recognition of User Account Types
We now analyse the results obtained with our classification algorithms. In the cross-validation phase, four attempts of classification were made with each algorithm: using only intertweet delay distributions (ITD); using only tweet time distributions (TT); using the joint distribution of both features as independent variables (JI); and using the joint distribution of both features as non-independent variables (JNI). Table 3 shows the percentage of correct classification for the 2-Classifier in each of the four trials, with 86 samples from each class. We can see from this table that using the marginal distribution for tweet time yielded better results than using the one for inter-tweet delay (78.5% vs. 71.5%), which is reasonable since the tweet time distributions, presented in Figure 6, exhibit particularly distinct shapes among the three tweeter classes. As one would expect, using both ITD and TT features yielded better results than using only one feature (83.1% vs. 71-79%). Moreover, the classifier using the joint distribution of the two variables under the independence assumption, with 83.1% correctness, generated better results than the one with the nonindependence assumption, with 82.6%. We believe this is due to subsampling of the joint distribution, which causes interpolation to be poor. Table 4 shows the percentage of correct classifications for the 3-Classifier, in which we used 67 samples from each class. The 3-Classifier performed slightly worse than the 2-Classifier due to the larger number of classes. From this table, we can see that again the tweet time marginal distribution led to better classification results than the inter-tweet delay distribution and that in the 3-Classifier this difference was even more pronounced (70.6% vs. 54.2%). Similarly, the variable independence assumption again yielded better results than the non-independence assumption (73.1% vs. 52.7%). The good performance under the independence assumption suggests that the tweet time and inter-tweet delay variables are  rather independent in terms of account class. To check for independence, we performed both Pearson and Kendall's t tests between the values obtained for these variables in each class. As shown in Table 5, the correlation values obtained for the two variables was very low in all cases, which proves that they are indeed independent. We also used separate training and test datasets in order to evaluate the performance of our most successful classification system, which uses the joint distribution of both inter-tweet delay and tweet time features as independent variables. The samples in the training set were used to generate the probability distributions, then each sample in the test set was classified following the procedure described in the Methods. Figure 9 shows the average percentage of correct classification obtained with the 2-Classifier and the 3-Classifier when varying the training dataset size from 5% to 70% of the total number of user accounts. Both classification algorithms are shown to be robust to the decreasing size of the training dataset.
To further verify the validity of ours models, we compared their performance to that of randomised models, created by randomly shuffling the data across the 3 classes, thus generating probability distributions that included data from all classes. The performance of these randomised models was therefore governed by chance, yielding an average 52:2%+18:2% correct classification for the randomised 2-Classifier and an average 32:3%+5:6% correct classification for the randomised 3-Classifier. Despite having an expected decrease in classification correctness as we decreased the size of the training dataset, both our algorithms performed well above the randomised models, even when the training dataset comprised only 30% of the samples (81.2% vs. 52.2% for the 2-Classifier, and 70.8% vs. 32.3% for the 3-Classifier).
Previous research applying content-based classification achieved correctness ratios from 82.8% to 94.9% when distinguishing   Table 2. doi:10.1371/journal.pone.0065774.g004 between human, bot and cyborg users [17]. In contrast, our approached using the timing of tweets alone resulted in an average 75:8%+4:8% correctness when distinguishing between personal, managed and automated accounts (using 70% of samples for training). The classification correctness percentage we have obtained is only slightly worse than those presented in related Power-law fits to the tail of each class inter-tweet delay distribution in terms of power-law exponent (mean a+SD and cut-off value x min above which powerlaw tails are observed). The p-value for the fit statistics was obtained by using the Kolmogorov-Smirnov statistic as a distance measure between the data and the fitted power-laws. doi:10.1371/journal.pone.0065774.t002  work, with two important advantages: we did not decide a priori what features were characteristic of each account class, and our classification is based solely on tweeting behaviour and does use any other account feature or require parsing of tweet contents. The fact that both our classifiers generated good results when operating under the assumption that the inter-tweet delay and tweet time variables are conditionally independent is an unexpected yet interesting result. Intuitively, it is reasonable to assume that these two variables are closely related, since the amount of time a user waits before tweeting must be influenced by the time of the day when their last tweet was posted. However, our results in both the classification algorithms and the correlation tests indicate that these two variables are in fact not so closely related, which could be explained by the existence of external factors which influence them more strongly than they influence each other. For the tweet time, this external factor is probably the daily routine of a user, which has a much bigger impact on the time of a user's tweet than the inter-tweet delay. For the inter-tweet delay, we conjecture that this factor is the universal laws that govern the timing of many human activities, as found in previous research in other modes of communication [10,14,15,21] and observed in our own data analysis.

Prediction of Next Tweet's Time Distribution
We can predict, using a simple method, the time to the next tweet, based on the time that has passed since the last tweet with good accuracy. We tested two ways of predicting tweet times, 1. a simple one using just the inter-tweet-time distribution for each class and 2. a more complex representation that accounts for the inter-tweet-time distribution on a given time-of-the-day. This means that the first prediction method would ignore for an individual's next tweet time prediction whether the person may be currently in the middle of their night or in their working hours. Surprisingly the first method is about as good in prediction performance as the second prediction method. In these algorithms, we used 67 samples from each class, and computed the coefficient of determination R 2 as a goodness of fit measure between our simple model and the data. Table 6 shows the average R 2 obtained for each account class by our predictive models. Differences between the two methods were negligible for individual and managed accounts (R 2~0 :66 vs 0:66 and R 2~0 :72 vs 0:72, respectively), or small for robot accounts (R 2~0 :52 vs 0:57).
The average R 2 results of 0:66+0:07 for personal accounts and 0:72+0:09 for managed accounts are in a good range for human data. In order to evaluate the statistical significance of these results, we applied the same predictive model (CDFs generated for each account class) to predict randomly generated data. We used a pseudo-random number generator, drawing numbers from a uniform distribution over range 1 to 1,000,000, thus creating random test samples. In these tests, the average R 2 obtained when measuring the fit of the CDF to the test samples' step functions was 0:33+0:08, which is much lower than any of the values obtained for the real data. Conversely, we used the randomly generated data to create a null model and used this model to predict our real test samples. In these tests, we again obtained average coefficients of determination well below our prediction results, as shown in Table 6. We can conclude that our results are statistically significant, but could potentially be improved by the use of additional information about the tweeting patterns observed.
To test the robustness of our single-distribution predictive model, we performed experiments using separate training and test datasets. We varied the sizes of the training and test sets, starting with 30% and 70% of the samples, respectively, then increasing the training set by 10% in each experiment, until we had 70% of samples for training and 30% for testing. In each of these set ups, we repeated the experiment 10 times, each time reshuffling the samples among each class. Table 7 shows the average R 2 obtained in these experiments. We can see from these results that the predictive model is robust to the decreasing size of the training dataset. Furthermore, the lower average R 2 and larger standard deviations obtained for the bot-controlled class indicate that the behaviour of these accounts in less predictable by the model than those of human-controlled accounts. This is expected, since botcontrolled accounts have programmed activities and are therefore less uniform in their behaviour.

Discussion
We have investigated the nature of broadcast communication by first developing a system to collect large-scale datasets from Twitter, then studying the behaviour of different types of user accounts: personal, belonging to a single individual; managed, belonging to a corporation; and bot-controlled, which are administered by a computer program. We examined the intertweet delay and tweet time distributions for each account class, and found that they present very distinct tweeting patterns, allowing us to distinguish between them in an automated manner. We also found that the distribution of a user's tweets throughout the day is closely related to their daily routine, and that the distribution of the inter-tweet delay, i.e., the time interval between two consecutive tweets by the same user, displays a power-law decrease in its tail. This last result agrees with and extends the findings of many other studies in Computational Social Science [9,10,14,15,21,22], reinforcing the idea that a bursty, fat-tailed behaviour is characteristic of the time of many human actions.
All three classes of Twitter accounts considered, even the botcontrolled one, did not exhibit a characteristic time scale in their tweet periodicity, but rather a scale-free behaviour. Characterising the fluctuations in tweet activity, we found that inter-tweet delay variability scales in proportion to the mean inter-tweet delay, known as signal-dependent noise in neuroscience [1]. This abstract decision to post a tweet thus shows the same characteristic variability structure of both neuronal and behavioural variability (reviewed in [1]). We found that the power-law distributions of inter-tweet delays, particularly the tails of the distribution, exhibit a pronounced difference across the three classes of Twitter Correct classification percentage for the 2-Classifier in four attempts during the cross-validation phase: using the marginal distribution for inter-tweet delay (ITD), using the marginal distribution for tweet time (TT), using the joint distribution of both properties as independent variables (JI), and using the joint distribution of both properties as non-independent variables (JNI). doi:10.1371/journal.pone.0065774.t003 Correct classification percentage for the 3-Classifier in four attempts during the cross-validation phase: using the marginal distribution for inter-tweet delay (ITD), using the marginal distribution for tweet time (TT), using the joint distribution of both properties as independent variables (JI), and using the joint distribution of both properties as non-independent variables (JNI). doi:10.1371/journal.pone.0065774.t004 accounts considered. Bot accounts describe a more unstructured tweet time behaviour both throughout the day and the week, when compared to human-driven accounts. Personal accounts tweeted more evenly throughout the week and on each day more tweets were recorded during typical awake time hours (7am to midnight). Managed accounts were more active during the 5 working days and during reasonable working hours (8am-8pm). Thus, our Twitter activity analysis showed that there are different patterns of tweeting activity across the Twitter account classes, suggesting that classification of account holders is possible without having to parse the content of tweets. We created two naive Bayes classification algorithms based on the empirical probability tweet time distributions, the first one to distinguish between personal and managed accounts, and the second one to classify all three types of accounts studied (personal, managed and bots). Both classifiers performed well, resulting in 84:6%+2:2% correctness for the 2-Classifier and 75:8%+4:8% for the 3-Classifier, with the best results being generated by the use of joint probability distributions of inter-tweet delays and tweet times, assuming independence of the two variables. The fact that our classification algorithms performed well under the assumption that these two properties are independent indicated that they are not closely related, which we have proved by performing correlation tests between the two variables for all three account classes. Previous research using contextual analysis and tweet content analysis achieved correctness ratios from 82.8% to 94.9% [17]. In contrast, our approached using tweet timing alone resulted in 75:8%+4:8% correctness when distinguishing between the three account classes studied.
Additionally, we implemented two predictive models in order to attempt predicting when the next tweet of a user would be posted. In these probabilistic models, we used the inter-tweet delay distribution of a given class in order to predict the next delay for a user of the same class. In our first attempt at probabilistic prediction, we used only the inter-tweet delay distribution of a given class in order to predict the next delay for a user of the same class. We then tried using separate distributions for each hour of the day, adding to our model the information about the time of the tweets. The use of separate prediction hours based on the time-ofday only marginally improved the prediction results, if at all. Interestingly, we were better able to predict human-driven nexttweet times than for the robot-driven accounts. Thus the fact that robot-driven tweet times are less predictable than human tweet times may be the result of a. bot-controlled accounts having programmed activities which vary considerably across individual bots and are therefore less uniform in their behaviour or b. the result of bot-controlled accounts being more driven in a reflexive mode responding to external events (e.g. news). To the best of our knowledge, there has been no previous research attempting to predict the timing of tweeting activity or related human activities, and we present our model results here as a first benchmark.
The identification and classification of specific types of users on Twitter can be useful for a variety of purposes, from the computational social sciences, focusing advertisement and political campaigns, to filtering spam, identity theft and malicious accounts. The occurrence of spamming and campaigning on Twitter has prompted several studies on methods for identifying certain types of behaviour that are characteristic of 'manipulators'. Chu et al. Figure 9. Classification correctness obtained with varying training dataset size. We evaluated the robustness of our classification algorithms by testing with different sizes for the training and test datasets. The horizontal axis shows the percentage of user accounts used for training, as well as the number of accounts used for training in the 2-Classifier (in blue) and in the 3-Classifier (in red). The remaining accounts were used for testing. Both algorithms perform well above a randomised model in all experiments, even when the training dataset comprised only 30% of the samples (81.2% vs. 52.2% for the 2-Classifier, and 70.8% vs. 32.3% for the 3-Classifier). In these experiments, we used the joint distribution of inter-tweet delay and tweet time as independent variables, and used a total of 86 accounts from each class in the 2-Classifier and 67 accounts from each class in the 3-Classifier. Each experiment was repeated 10 times, and at each time the samples were randomly shuffled among each class. doi:10.1371/journal.pone.0065774.g009 Average + SD coefficient of determination (R 2 ) obtained for each class by the two probabilistic prediction models during cross-validation. We compare the performance of our models to the results of a null model, which was created with random samples generated from a uniform distribution over range 1 to 1,000,000. doi:10.1371/journal.pone.0065774.t006 Table 7. Tweet-time predictive model R 2 for varying training set sizes.

Personal
Managed Bot Average + SD for the coefficient of determination (R 2 ) obtained for each class by the predictive model when varying the size of the training dataset, starting with 30% of samples and increasing up to 70% of samples (the remaining samples were used for testing). doi:10.1371/journal.pone.0065774.t007 [17] investigated the differences between Twitter accounts controlled by humans, bots, and cyborgs by studying the message content of tweeting behaviour, tweet content and account properties. Despite attaining a high classification success rate, the system heavily relies on processing the contents of tweets, which can be expensive and amenable to manipulation as in email-based spam. Similarly, Lumezanu et al. [16] investigated how Twitter is used to spread propaganda by studying the Twitter behaviour of ''hyperadvocates''. In contrast to these studies, which rely on context and content-parsing to operate, we used the timing of tweet actions as only variable, abstracting away from complications of natural language processing and context-factors. We were able to reliably detect the difference between true individuals and public-relations managed accounts -ultimately due to the nature of one being personal activity, the other resulting from an employment type activity. Moreover, we can reliably distinguish these human tweeters from robot-based tweeters on their relative tweet-timing distributions.
In the context of computational approaches to study human behaviour, we have used a free, publicly harvestable resource to study human behaviour patterns. We have measured and shown that Twitter-using individuals have a distinct and characteristic structure in their tweeting behaviour, characterised by the tails of their inter-tweet time distribution and their rather more unstructured hourly tweet probability. Related work in Computational Social Science [10,14,15] has been concerned with the timing of peer-to-peer human communication, such as emails, letters and phone calls, for which the power-law slopes obtained were between 1 and 1.5, while our results show 2:38+0:059 for personal accounts and 2:88+0:16 for managed accounts. In contrast to previous studies, we have obtained results for broadcast communication that extend the general conclusions about the nature of human communication behaviour to this more novel form of personal communication. Our findings may be easily applied and extended to other forms of broadcast communication in public spaces, be it social networks or information sources such as blogs. Our work suggests that inter-communication intervals may show characteristic scaling-law exponents in human broadcast communication and may also be applied to the analysis of animal and plant broadcast communication timings, as in the case of mating calls or chemical signals. The finding that individual communication and broadcast communication are markedly different in human electronic communication may suggest that different (neuronal) mechanisms are at play in decision making about communication initiation, however this would need to be verified and compared to non-electronic forms of interaction. We note that the inter-event statistics of electrical impulses (spikes) of single neurons exhibit the same variability structure and powerlaw tails in their inter-event statistics. Thus, some of the statistical features we observe in our broadcast communication data and others in peer-to-peer communication, may be a more general feature of distributed communication networks, applicable from neural circuits to human society.


extfieldsofstudy: Medicine Physics

text:
Record-Breaking Avalanches in Driven Threshold Systems

Record-breaking avalanches generated by the dynamics of several driven nonlinear threshold models are studied. Such systems are characterized by intermittent behavior, where slow buildup of energy is punctuated by an abrupt release of energy through avalanche events which usually follow scale invariant statistics. From the simulations of these systems it is possible to extract sequences of record-breaking avalanches, where each subsequent record-breaking event is larger in magnitude than all previous events. In the present work, several cellular automata are analyzed among them the sandpile model, Manna model, Olami-Feder-Christensen (OFC) model, and the forest-fire model to investigate the record-breaking statistics of model avalanches which exhibit temporal and spatial correlations. Several statistical measures of record-breaking events are derived analytically and confirmed through numerical simulations. The statistics of record-breaking avalanches for the four models are compared to that of record-breaking events extracted from the sequences of independent identically distributed (\emph{i.i.d.}) random variables. It is found that the statistics of record-breaking avalanches for the above cellular automata exhibit behavior different from that observed for \emph{i.i.d.} random variables which in turn can be used to characterize complex spatio-temporal dynamics. The most pronounced deviations are observed in the case of the OFC model with a strong dependence on the conservation parameter of the model. This indicates that avalanches in the OFC model are not independent and exhibit spatio-temporal correlations.

I. INTRODUCTION
A record-breaking event is an event that is larger (smaller) than all previous events in the observational sequence of some physical measurements. Sequences of record-breaking events can be extracted for various natural phenomena where records of the variability of different measurements exist. The best known example is weather, where daily observations of temperatures are documented and new record-breaking temperatures are often reported. Monitoring such record-breaking temperatures can be used to infer possible trends in the weather variability [1][2][3][4][5]. For example, the excess of maximum record-breaking temperatures over minimum record-breaking ones is taken to be direct evidence of global warming. Other natural phenomena for which record-breaking events play an important role include earthquakes [6][7][8][9][10], floods [11], forest-fires, volcanic eruptions, and solar events [12], to name only a few.
Record-breaking events or records form a subsequence of all recorded measurements and are only distinguished based on their sizes or magnitudes. Since a record-breaking event is defined as the largest event up to its time of occurrence, a subsequent event becomes the next record-breaking event only if it exceeds the size of the previous record and all events in between are smaller than both records [1,4,[13][14][15]. The mathematical framework to analyze record-breaking events is related to extreme value statistics and in most cases assumes events which are independent and identically distributed (i.i.d.) [13][14][15][16]. In this case several statistical properties of records are independent of their underlying distribution and yield to analytical treatment [17]. This standard treatment of records assumes that the underlying events occur equidistant in time. As shown in Ref. [18], one can generalize the analysis of record-breaking events to situations in which the occurrence of events grows stochastically in time. It was found that such a temporal modification preserves several fundamental properties of the original model, i.e. independence from the underlying distribution of the i.i.d. events.
In contrast the situation can be significantly different if one considers more complicated stochastic processes. In recent years some progress has been made in analyzing records and record-breaking events in non-i.i.d. random sequences including sequences with time-varying underlining distributions as well as sequences with strongly correlated events. For example, a theoretical and observational analysis was performed for daily record temperatures in Philadelphia to investigate the effects of temperature trends and correlations [3]. In that work the authors found that the small rate of increase in the mean temperatures does not significantly influence the statistics of record-breaking events over the time span of the observations. The authors also considered the role of temporal correlations between events by generating synthetic time series with correlations decaying as a power-law. It was reported that the frequency and magnitude of the record-breaking events was almost identical to those obtained using i.i.d. random variables and were not sensitive to the power-law exponent of the correlation function. Events drawn from independent random variables but with progressively broadening or sharpening distributions were considered in [5]. It was shown that for a power-law growth of the width of the distribution three regimes can be identified if one considers the growth of the mean number of records. These regimes correspond to the asymptotic behavior of the tail of the distribution according to Gumbel, Fréchet, and Weibull classes often studied in the context of extreme value statistics [19]. To further investigate the effects of correlations in time series, record-breaking events in sequences generated by random walks and Lévy flights were also considered [20]. It was found that the statistics of records are independent of the distribution of jumps. The average number of records and its standard deviation grow as a square root of the number of time steps which implies that fluctuations remain large and of the same order as the mean. For related stationary and non-stationary stochastic processes with long-range memory, the situation is more complicated but memory typically leads to significant deviations from the i.i.d. behavior [12,21].
Despite this progress for some non-i.i.d. processes, we are still far away from a general understanding of recordbreaking events and their properties in complex systems that are often observed in nature. To understand such systems and specifically the emergence of complex patterns and long-range spatio-temporal correlations in physical systems, various cellular automata have been studied in the past. Usually such models possess rather simple and intuitive rules of evolution which mimic real observed properties of dissipative complex systems. Of particular interest are driven nonlinear threshold systems exhibiting avalanche-like behavior with scale invariant statistics where effects of correlations and self-organization play a crucial role [22,23]. A plethora of models has been introduced to illustrate such behavior. Among them the sandpile model [24,25], Manna model [26], the Olami-Feder-Christensen model [27], and the forest-fire model [28,29] are the most prominent ones. They all exhibit scale-invariant behavior in the distribution of avalanches with respect to their size and duration when driven into a steady state. They also possess nontrivial spatio-temporal correlation patterns and serve as conceptual models to explain such diverse phenomena as earthquakes, forest-fires, landslides, etc.
In the present work, we investigate the influence of spatio-temporal correlations on records by studying the recordbreaking statistics of avalanches in the four cellular automata described above. For each of the models we study the record statistics both in magnitude and time domains by constructing various distributions and the corresponding averages. Our numerical simulations show that there are significant differences to what would be expected for i.i.d. random variables. The most pronounced deviations are observed in the case of the OFC model with a strong depen-dence on the conservation parameter of the model. In addition to numerical simulations and analysis of records of the four cellular automata, we also derive an exact analytical expression, Eq. (18), for the probability that the kth record is broken after m time steps in the case of i.i.d. random variables. This allows us a more rigorous comparison between the models and the i.i.d. case.
The paper is organized as follows. In the next section we introduce basic known facts concerning the statistics of record-breaking events extracted from sequences of i.i.d. random variables. We also present the analysis of recordbreaking events generated from the truncated power-law distribution which is relevant for our analysis of threshold cellular automata. In Section III we recall the definition of the four cellular automata which we have studied to infer the structure of their record-breaking avalanches. In Section IV we report on the results of computer simulations of those models. Section V concludes our analysis.

II. RECORD-BREAKING STATISTICS
Consider a sequence of measurements, {m(t i )}, of a particular outcome of a physical experiment or computer simulation of a model over a period of time, where t i are the times of measurements. In general, this represents a stochastic variable. Examples include temperature measurements, earthquake magnitudes, flood areas, forest-fire areas, etc. A record-breaking event x(t n ) up to time t n , is defined as the event with the largest value among all previous measurements, x(t n ) = max{m(t 1 ), m(t 2 ), ..., m(t n−1 )} [13]. Subsequent measurements can produce recordbreaking events if they exceed in value the current record-breaking event. Therefore, for any historical data or for ongoing observations the sequence of record-breaking events can be extracted and analyzed. In the analysis which follows, we use discrete time steps n = 1, 2, 3, ... which define the moments when measurements are done or time steps when random variables are generated. In this case we use a simplified notation for records: x(t n ) ≡ x(n). To enumerate the record-breaking events in a sequence we use a subscript k = 1, 2, 3, ... to indicate their record order. For example, x k indicates the value of the kth record-breaking event, therefore, for a given realization of observations one has the following sequence of record-breaking events: {x k } = x 1 , x 2 , x 3 , .... To specify in addition the time of the occurrence of the kth record we write: x k (n).
In the most simple case it is assumed that the records are extracted from the sequence of i.i.d. random variables drawn from a particular probability density function f (x). This density can be specified a priori from known distributions or it can be the outcome of the measurements of natural systems or the computer simulations of model systems. It can be bounded or unbounded depending on the system under consideration. For the bounded distributions the corresponding record-breaking events will be bounded as well. The probability that the records will not exceed x is given by where x min is the lower bound of the support over which the distribution is defined.
To quantify different aspects of record-breaking events several fundamental statistical measures can be estimated. The existing theory of record-breaking events is mostly developed for the i.i.d. random variables drawn from a particular distribution [14,15]. In the rest of the section we assume such kind of random variables.

A. The average number of record-breaking events
The average number of record-breaking events, N n , characterizes the occurrence of records in time. It can be estimated by averaging over ensemble realizations of sequences of random variables drawn from a particular distribution with the probability density function f (x). Consider a random variable N n which is the number of records that occurred up to a time step n. This number can be expressed as where I j is the indicator function given by , if a record occurs in the jth time step, 0, otherwise.
Assuming that the events are i.i.d. random variables the probability of the occurrence of the jth record-breaking event is simply P j = 1 j [14]. This represents the harmonic decrease in time of the probability of the occurrence of record-breaking events. Therefore, the mean value of I j can be computed using this probability: Using Eqs. (2) and (4) the average number of records, N n , can be evaluated where H n is a harmonic number, H n = n j=1 1 j , and γ ≈ 0.577215665... is the Euler-Mascheroni constant [30]. This result indicates that the average number of record-breaking events grows as N n ∼ ln n for large values of n and is independent of the underlying probability density function f (x). For record-breaking events generated from processes with memories or long range correlations the average number can deviate from Eq. (5).
The variance of the average number of record-breaking events N n is given by [5] Var For i.i.d. random variables the ratio of the variance, Eq. (6), to the mean, Eq. (5), tends to unity as n → ∞ and the distribution of the N n becomes Poisson distribution with the mean ln n. As a consequence the record times can be described by a log-Poisson process [5]. In this case the probability that k records have occurred up to a time step n is For large values of n the distribution for the number of records N n tends to the Gaussian distribution.

B. Distribution of magnitudes
The probability density function that the kth record-breaking event is equal to x has the form [3]: where F (x) is the distribution function for the random variable x given by Eq. (1). Eq. (8) represents a recursive formula to compute the distribution of magnitudes for the kth record-breaking event by knowing the distribution for the (k −1)st record. The p 1 is simply equal to the underlying distribution from which the random variables are drawn, . It is also possible to evaluate Eq. (8) for k = 2, 3, .... For example, by noticing that p 1 (x) = f (x) = dF dx , one has for k = 2 One can prove by induction that the general form of Eq. (8) for any arbitrarily k is Eq. (10) is valid for i.i.d. random variables drawn from any underlying distribution function F (x) [14].
The average magnitude, x k , of the kth record-breaking event can be estimated using Eq. (10) where S = [x min , x max ] is the support of the distribution function F (x). It is also possible to define the average magnitude x(n) of a record-breaking event at a time step n: where q(x, n) is the probability density function for the records at a time step n. So that the probability to find a record-breaking event to be between x and x + dx at time n is q(x, n)dx. This probability density function q(x, n) can be expressed in terms of F (x) as [4,14] q Noticing that f (x) = dF dx , Eq. (12) can be rewritten as follows It is important to note that the results, Eqs. (10)- (14), are valid only for i.i.d. random variables. In Section IV we consider simulations of the four cellular automata and compare the obtained results for record-breaking avalanches with the behavior given by Eqs. (10)- (14) in order to quantify the deviations from the i.i.d. case.

C. Distribution of interevent times
The interoccurrence or interevent times between two subsequent record-breaking events is an important measure which gives an estimate as for how long one has to wait for the next record to occur. It was shown that for i.i.d. random variables the distribution is independent of the underlying distribution of the record-breaking events and the non-normalized histogram follows a simple power law: G(m) = 1/m, where m is the time difference between kth and (k − 1)st record-breaking event, m = t k − t k−1 [4]. This distribution is constructed by incorporating all interevent times including those early in the sequence as well as those later in the sequence of records.
A more detailed structure of time intervals between record-breaking events is given by the distribution, w k (m), which is the probability that the kth record is broken after m time steps. The formal expression for this quantity in the case of i.i.d. random variables is [3,15] where F (x) is the underlying cumulative distribution function for the random variables from which records are drawn.
gives the probability that a new record, (k + 1)st, is set after m ≥ 1 time steps after the kth record of value x k . To obtain Eq. (15) one has to average this probability over all possible values of x. It is possible to evaluate Eq. (15) for the 1st record (k = 1) to be broken after m time steps, w 1 (m). This can be done explicitly by noticing that p 1 (x) = f (x) = dF dx . Substituting this into Eq. (15) and performing integration by parts one obtains Due to power-law nature of the distribution, Eq. (16), the average waiting time between the first record and the second record is infinite, m 1 = ∞ m=1 m w 1 (m) = ∞. Using Eq. (10) the Eq. (15) can be rewritten in terms of the cumulative distribution function F (x) as follows The integration can be carried out explicitly (Appendix A) with the result: For k = 2 and 3 we have where H m is a harmonic number and Ψ(m) is a digamma function [30]. These results show that the probability distribution w k (m) is independent of the underlying distribution F (x) as long as the random variables are i.i.d..
The derived equation for w k (m) can be used to compute the non-normalized histogram G(m) as follows where to evaluate the integral we exchanged the operation of summation and integration and used the known result for the expansion of the exponent into the infinite sum Using Eq. (22) the integral in Eq. (21) can be evaluated easily. The probability distribution, u k (n), for the time of occurrence of the kth record-breaking event at a time step n can also be analyzed. Using u k (n) one can compute the average time n k = ∞ n=1 n u k (n) of the occurrence of the kth record-breaking event. The first event is always a record-breaking event so n 1 = 1. The distribution for the second record-breaking event u 2 (n) coincides with the distribution w 1 (n) which defines the distribution of interevent times between the first record (k = 1) and the second record (k = 2).
The time of the occurrence of the kth record, n k , is a random variable and can be represented as a sum of two random variables, n k = n k−1 + m k , where n k−1 is the time of the (k − 1)st record and m k is the interevent time between records k and k − 1. Therefore, the density u k (n) can be obtained recursively from the discrete convolution of the two distributions u k−1 (n) and w k (m): The above defined distributions and averages can be used to describe statistical properties of record-breaking events generated by different models and processes and documented in observations of various natural systems. Based on their behavior it is possible to infer the structure of record-breaking events and underlying processes which generate them.

D. Power-law distributed random variables
Many natural phenomena with nonlinear threshold dynamics exhibit scale-invariant statistics in their distribution of output events (avalanches). These statistics can often be approximated by a power-law distribution. One prominent example is earthquakes where the frequency-magnitude statistics follows Gutenberg-Richter power-law scaling [31]. The probability density function for a power-law distribution can be written as where γ is an exponent. Usually observed physical processes exhibit power-law scaling only for several orders of magnitude and are truncated on both ends. This introduces finite-size effects into distributions and can lead to non-trivial scaling behavior. For this analysis we consider the power-law distribution with a bounded support [x min , x max ], where x min ≥ 0 and γ = 1. In this case the distribution is well defined and can be normalized. The probability density function f (x) and the distribution function F (x) are: For i.i.d. random variables the distribution function of the kth record-breaking event, p k (x), can be evaluated using Eqs. (8) and (26). Using Eq. (11) one can compute the mean value x k of the magnitude of the kth record-breaking event. The integration can be performed analytically for specific values of k but the expressions are rather cumbersome so one can evaluate them numerically. Similarly the average magnitude x(n) of the record-breaking events at a time step n can be evaluated numerically using Eqs. (14) and (26).

III. DRIVEN NONLINEAR THRESHOLD MODELS
Observations of many natural phenomena reveal that the behavior of these systems is governed by the threshold dynamics where a slow build up of energy is punctuated by the fast release of the accumulated energy through avalanches. The analysis of the frequency-magnitude statistics of these events shows that in many instances the distributions follow power-law statistics with anomalous exponents. These systems also exhibit non-trivial correlations in time and space. To replicate such behavior a plethora of cellular automata has been proposed. Among them the sandpile model [24,25], Manna model [26], the Olami-Feder-Christensen model [27], and the forest-fire model [28,29] play a prominent role. These models were also introduced as examples to illustrate the hypothesis of self-organized criticality (SOC) where systems evolve into a steady state which lack any characteristic length and time scales [22,24,25]. The initial assumption was that this state can be reached without fine tuning of any model parameters, i.e. the model dynamics organizes the system into a scale-free state. Subsequent analyses of the above and other models revealed that the tuning of parameters is present even in the case of "parameterless" models such as the sandpile model. In many cases the tuning parameter is the separation of slow and fast time scales present in the models though it is still debatable to which extent such a separation is necessary [32]. In the limit of large separation (slow driving and fast release of energy) some systems exhibit power-law tails in their avalanche size distributions.
Another feature which plays an important role in the dynamics of these models is a local conservation of energy. During avalanche events in some models the energy is conserved at each toppling site and only allowed to dissipate at the boundaries (the sandpile model, Manna model). Other models introduce a possibility of dissipation at interior sites (OFC model, fores-fire model). In some cases presence of local dissipation can alter the critical state.
A. The sandpile model The sandpile model was introduced by [24,25] as a paradigm of SOC to illustrate the emergence of scale invariant behavior in driven nonlinear threshold systems. The model is defined on a graph where each vertex has a height variable z ij ∈ {1, 2, ...} representing the amount of sand grains. The system is driven by slowly adding new sand grains at random into the system. Each vertex is also associated with a critical height z c . When a sand grain is redistributed to the vertex with z c amount of sand present the vertex becomes unstable and topples redistributing sand grains to its nearest-neighbors. This process can initiate an avalanche of toppling events. Starting from any initial distribution of heights after a transient period of time the model self-organizes into a steady state where the distribution of avalanches with respect to their sizes and durations follow a power-law behavior and long-range spatial correlations emerge [33]. This state has been identified as a SOC state in analogy to the second-order phase transitions encountered in statistical mechanics systems.

B. The Manna model
The Manna model is another example of a model exhibiting a SOC type behavior [26]. It was introduced as a model for particle collision and repulsion. The dynamics is pretty similar to the sandpile model except the height variable describes only two states: empty or occupied. If a new particle arrives at an empty site the site becomes occupied. Otherwise if the site is already occupied by a particle the new particle makes it unstable and two particles are transferred randomly to its nearest neighbors. This dynamics mimics the collision and repulsion of particles without conservation of momentum. After passing the transient regime the model enters into a steady state where distribution of avalanches display power-law behavior with well defined critical exponents [26].

C. The Olami-Feder-Christensen (OFC) model
The Olami-Feder-Christensen (OFC) model was introduced by [27] to analyze the behavior of massless blocks interconnected by springs and pulled over a surface with friction. In turn, the spring-block model [34] can be considered as a simple realization of a fault surface with threshold avalanche events representing the analog of earthquakes. The OFC model is locally nonconservative which results in different regimes of behavior depending on the degree of non-conservation specified by the parameter α.
In the model, each site (i, j) of a square lattice is assigned a continuous variable h ij ∈ [0, h c ] which defines a stress level. The model is driven uniformly with a constant rate until one of the sites reaches the threshold value h c = 1.0. This site becomes unstable and redistributes stress to its nearest neighbors and is reset to zero. The amount of stress transferred to each nearest neighbor site is αh ij , where 0 < α ≤ 0.25 defines a degree of non-conservation. This transfer may cause some of the neighbors to topple initiating an avalanche. The model is locally non-conservative for α < 0.25. The parameter α plays an important role in the dynamics of the model with smaller values of α leading to longer transient regimes [35] and also affecting its critical properties. However, the question of criticality of the OFC model still remains open [35].
The model produces rich behavior that can be explored in order to understand its dynamics. After passing the transient regime the model forms highly correlated spatial structures with power-law distributed patches of different heights. The existence of non-conservation in the model introduces several scales over which the distributions of avalanches are different at different scales. Large patches are mostly responsible for large avalanches while small avalanches are generated by an interplay between small patches [35].
Recent results concerning the spatio-temporal organization of avalanches in the OFC model show that the model generates, in addition to regular avalanches, near periodic avalanches associated with asperity structures formed as a highly concentrated stress state [36,37]. These asperity events form spatio-temporal clusters with pronounced characteristic features. The processes of synchronization and desynchronization in the OFC model with short-term correlations among events were also observed [38]. It was suggested that synchronization in the model drives it into a state with scale-invariant power-law distributed events while desynchronization results in the occurrence of aftershocks and foreshocks.
It was also suggested that the OFC model can exhibit the maximum in the entropy production rate for the range of the parameter α ∈ [0.15, 0.2]. This was associated with finite spatial memory of past avalanches associated with asperity events [39][40][41]. There is a tendency of similar large avalanches to reoccur as the model is driven forward. These avalanches rupture approximately the same area and occur quasi-periodically in time. The range of the parameter α where the maximum of the entropy production rate is achieved corresponds to the subcritical regime of the OFC model. A critical point was suggested to exist at α = 0.25 with a diverging correlation length and average avalanche sizes [41].

D. The forest-fire model
A model which resembles the occurrence and propagation of natural forest fires was proposed in [28,29]. The model is considered on a two-dimensional square lattice where each site can be occupied by a green tree or be empty. The dynamics of the model defines how trees are planted and how fires destroy the trees. Two parameters define the behavior of the model. A rate of planting of new trees, p, and a firing rate f p with which lightning can strike a site and start a fire. At each time step a random site is chosen. With the probability p a tree is planted and with the probability f lightning strikes a green tree. The fast time scale defines how model fires propagate in the system. A green tree adjacent to a nearest-neighbor burning tree will catch fire. This processes of spreading fires can initiate an avalanche which propagates until all connected green trees are burnt. It is assumed in the model that the process of the burning of fires is instantaneous. The control parameter in the model is the ratio θ = p/f [42].

IV. MODEL SIMULATIONS AND RESULTS
One of the main goals of this paper is to investigate the statistics of record-breaking avalanches in the above mentioned four cellular automata. Since these cellular automata are simple conceptual models mimicking real threshold systems with avalanche type dynamics, this will help us to answer questions related to the occurrence of recordbreaking events in models and real systems which display spatio-temporal correlations and non-trivial power-law scaling.
In order to perform numerical simulations of the above mentioned four cellular automata we used a square lattice of size L × L with varying L = 64, 128, 256 and performed in most cases 10 4 ensemble realizations of sequences of 10 5 − 10 6 time steps, where time steps were defined as the time of slow driving of the models not taking into account the time of propagation of avalanches. The statistics was collected after each model entered into a steady state to exclude transient effects. For each of the four models we constructed the various statistical measures of record-breaking events defined in Section II. The avalanche size for the sandpile, Manna, and forest-fire models was defined as the number of sites that toppled at least once during an avalanche. For these models the value of the avalanches were bounded by the lattice sizes. As a result, the record-breaking avalanches were bounded from above.
When simulating the OFC model for smaller values of the conservation parameter α the model generates a significant number of avalanches where only a single site relaxes. This directly affects the statistics of record-breaking events for early records since it leads to a non-negligible probability that ties among record values occur. Such ties do also occur for later records. To eliminate ties as well as the pathological limit states we introduced a low amplitude noise, η, into our simulations: during each toppling or relaxation event, a random number was drawn from a uniform distribution in the interval [−0.00005; 0.00005] and the toppled site was reduced by the value of its height plus/minus this noise η. This allowed us to eliminate ties by defining the size of an avalanche in the OFC model as the total amount of stress redistributed and dissipated during an avalanche for α < 0.25.

A. Results: The average number of record-breaking events
The average numbers of record-breaking events N n versus time step n obtained from the simulations of the four cellular automata are given in Fig. 1a. For the sandpile, Manna, and forest-fire models the behavior is close to the record-breaking events drawn from the i.i.d. random variables, Eq. To model the asymptotic behavior of the average number of record-breaking avalanches, N n , for large time steps we propose the following asymptotic equation where A, B and δ are constants. We fitted this equation to the average number of records obtained from the model simulations with the varying conservation parameter α and lattice sizes L = 64, 128, and 256. One such fit is given in Fig. 1a as a solid curve. The analysis indicates that the exponent δ increases with decreasing values of α. Then it reaches a plateau and starts to decrease moderately. To illustrate this the exponent δ is plotted versus the conservation parameter α in Fig. 3a. This figure also illustrates the dependence of the exponent δ on the lattice size L. We also plot the dependence of the parameter A on δ in Fig. 3b, that shows a similar trend. When the parameter α tends to the limiting value of 0.25, which corresponds to the full local conservation of energy in the OFC model, the average number of record-breaking events approaches the behavior corresponding to the records drawn from i.i.d. random variables and given by Eq. (5).
To confirm that the observed behavior for the OFC model was due to correlations among avalanche magnitudes we reshuffled the original sequences of avalanches which effectively destroyed any correlations among magnitudes. For these reshuffled sequences the average number of record-breaking avalanches followed exactly Eq. (5) for all considered values of α (not shown).

B. Results: Statistics of magnitudes of record-breaking avalanches
The avalanche frequency-magnitude statistics in the steady state for the most of the models follow power-law trend with finite-size effects. Both the sandpile and Manna models produce well defined power-law distributions for the sizes of avalanches with known power-law exponents [33,43]. The forest-fire model exhibits some variation in the distributions depending on the firing rate θ with presence of system wide fires for high values of θ and characteristic fires for small values of θ [42]. The behavior of the OFC model varies for different values of the conservation parameter α [35].
In this subsection we consider the distribution of record-breaking avalanche magnitudes. For the sandpile, Manna, and forest-fire models by magnitude we assume the area of toppled sites during an avalanche whereas for the OFC model it is the amount of stress transferred during toppling events in one avalanche as defined earlier. The distributions of the magnitudes of the kth record-breaking event, p k (x), are given for the OFC model with α = 0.18 in Fig. 4a. In the same figure we also plot as solid curves the theoretical distributions obtained from the record-breaking events drawn from the truncated power-law distribution, Eq. (25), and given by Eq. (10) with γ = 1.99, x min = 1.0, and x max = 256 2 . The model distributions deviate from the theoretical curves for several reasons. It is assumed that the magnitudes are not i.i.d. random variables. It is also evident from the distributions for avalanche magnitudes that they are not precisely power-law distributed. The average magnitudes x k of the kth record-breaking event for the four models are shown in Fig. 4b. These are compared for reference with the theoretical results applicable for i.i.d. random variables drawn from the truncated power-law distribution, Eq. (25), and given by Eq. (11). For the theoretical curves we used x max = 128 2 and the corresponding values of the exponent γ estimated from model simulations. These comparisons with the i.i.d. random variables drawn from the truncated power-law distribution are for illustrative purposes to show how model simulations deviate from them.
It is appropriate to note that when the record order k tends to infinity one expects that the distribution of magnitudes of the records has to converge to one of the three known classes (Gumbel, Fréchet, and Weibull) for the statistics of extreme events. Particularly, for the records drawn from the bounded distributions the distributions have to converge to a Weibull distribution. Another measure which characterizes the records is the average magnitude of record-breaking avalanches x(n) at a time step n. These have been computed for the four models and are plotted in Fig. 4c. It is observed that the average magnitude x(n) increases monotonically and approaches the upper bound set by the system size L 2 . For comparison we also plot as solid curves the average magnitudes obtained from Eq. (14) with the power-law distribution function Eq. (26). It is evident from the plots that the distributions and averages deviate from the theoretical curves. This can be attributed to the fact that the frequency magnitude distributions of the models on finite latices are not exactly given by Eq. (26) and the avalanches are not strictly i.i.d. random events. We also plot in Fig. 4d the average magnitude x(n) of record-breaking avalanches for the OFC model with the varying conservation parameter α and a fixed lattice size L = 128. The curves show a prominent trend where they experience a change in slope when going from small avalanches to large ones. This break in the behavior approximately corresponds to the transition observed in the average number of record-breaking events shown in Fig. 2a.
To confirm the presence of correlations between avalanche sizes in the OFC model we considered the size differences between subsequent avalanches, ∆a i = log 10 (a i+1 ) − log 10 (a i ), where i is the time step. A similar approach was used to study correlations between earthquake magnitudes [44][45][46]. In the absence of correlations the distribution of ∆a's should not deviate significantly from the distribution of the size differences constructed from the reshuffled catalogs, ∆a * i = log 10 (a i * ) − log 10 (a i ), where i * corresponds to the randomly picked index among N avalanches considered. In order to quantify the existence of correlations between avalanche sizes we considered the difference between cumulative distributions: ∆P (d) = P (∆a < d) − P (∆a * < d). This is illustrated in Fig. 5a for the four models. The dashed curve that fluctuates around 0 corresponds to the randomly reshuffled catalog where size correlations were destroyed. The error bars are given for 3σ confidence intervals. The OFC model displays strong deviation from the uncorrelated case indicating that avalanches in the OFC model are correlated. The deviations imply in particular that avalanches in the OFC model have the tendency to be followed by ones of similar size. On the other hand the avalanches in the forest-fire model show the opposite trend. This is expected for the forest-fire model as large fire avalanches are usually followed by small fires. This anti-correlation of avalanches for the forest-fire model does not change the statistics of records significantly [12]. The size differences between successive avalanches do not show strong correlations for the sandpile and Manna models. In Fig. 5b we also plot ∆P (d) for the OFC model for different values of the conservation parameter α. The strongest correlations are observed for the values of α below ≈ 0.16 which also corresponds to the higher values of δ observed in Fig. 3a.

C. Results: Interevent time distributions
Here we summarize the results obtained for the temporal structure of record-breaking avalanches of the four cellular automata. We computed the distribution that the kth record was broken after m time steps. This is shown in Fig. 6a for the OFC model with α = 0.17 for the first several orders of k = 1, 2, ..., 10. For comparison we plot the densities computed for the record-breaking events drawn from i.i.d. random variables given by Eqs. (19) which are independent of the underlying distribution F (x) of records. The comparison shows that the distributions for the OFC model deviate from the i.i.d. case for higher orders of k and large interevent time steps m. The other three models do not show such a departure and the corresponding distributions resemble the i.i.d. case.
The results of simulations for the interevent times between all subsequent record-breaking events for all four models are given in Fig. 6b. The normalized distributions, g(m), for the sandpile, Manna, and forest-fire models follow a simple power law: ∼ 1/m which is similar to the case of record-breaking events drawn from i.i.d. random variables [4]. These distributions were constructed by counting all interevent times for all record orders k. The OFC model shows departure from the i.i.d. case where one can observe a change in the behavior for small and large interevent times and strong dependence on the conservation parameter α (Fig. 6c). The probability density functions, u k (n), for the time of the occurrence of the kth record-breaking event versus a time step n are shown in Fig. 7a for the OFC model with α = 0.17 for the first several orders of k = 1, 2, ..., 10. In addition, the average time n k of the occurrence of the kth record-breaking avalanche is given in Fig. 7b. For comparison we also plot the average times obtained from the simulations of i.i.d. random variables for sequences with 10 6 time steps. This also illustrates the effect of the finiteness of the sequences on the average numbers. Fig. 7b demonstrates that the results for the four cellular automate are pretty close to the i.i.d. case with some deviations for the OFC model. On the other hand the finiteness of the sequence plays an important role as the average times of the occurrence are limited by this upper time limit that affects the average times for larger record orders k. In Fig. 7c we also plot the average time n k for the OFC model with the varying parameter α.

V. CONCLUDING REMARKS
Record-breaking events can be extracted from various natural and model systems. Their statistics can be used to analyze underlying trends and effects of correlations in the dynamics of those systems. Previous analyses of recordbreaking events drawn from i.i.d. random variables have revealed robust and universal behavior of temporal statistics of records independent of the underlying distribution functions from which the random variables were generated. The main goal of this paper is to find deviations from that universal behavior by studying record-breaking avalanches generated by driven nonlinear threshold systems. Particularly, we have considered four prominent cellular automata: the sandpile, Manna, OFC, and forest-fire models. In addition, we have also derived an exact formula, Eq. (18), for the probability that the kth record is broken after m time steps. The performed analysis reported in this paper revealed that among four models studied the OFC model displayed the most significant deviations in the statistics of record-breaking events compared to ones drawn from i.i.d. random variables. It was found for the OFC model that the average number of records for large avalanches can be described by Eq. (27) with the exponent δ being dependent on the conservation parameter α. This result signifies the departure from the case of i.i.d. random variables and is the direct effect of non-trivial correlations between avalanches in the OFC model. The other three cellular automata (the sandpile, Manna, and forest-fire models) do not show such a drastic departure indicating that their record-breaking avalanches are weakly correlated.
The distributions of record magnitudes for the four models also displayed deviations from the i.i.d. case. These distributions are controlled by the statistics of avalanche sizes which are not pure truncated power-laws with strong finite-size effects and dependence on model parameters. The correlations between avalanches also contribute to the observed behavior. This results in a more complex form of the distributions of record magnitudes and the corresponding averages.
The temporal structure of the occurrence of records also displays deviations compared to the i.i.d. case. Specifically, the distributions of the times of the occurrence of the kth record, u k (n), and the distribution for the kth record being broken after m time steps, w k (m), computed for the OFC model display departure from the corresponding distribution constructed from the i.i.d. random variables. On the other hand, the distribution of the interevent times computed for all record orders k do not show strong departure as well as the average time of the occurrence of the kth record, n k .
The finite-size effects both in time and magnitude domains play a prominent role in model simulations. They strongly affect statistics of records and have to be taken into account. The finiteness of the lattice introduces an upper bound for the magnitude of record-breaking avalanches. More thorough understanding of these effects requires development of more sophisticated simulation techniques in order to generate enough statistics of large records.
The OFC model can be considered as a simplification of a single fault system to simulate earthquake occurrences.
The avalanches in the model reproduce some important aspects of statistics of earthquakes such as a power-law frequency-magnitude distribution. Therefore, the analysis of correlations in the OFC model can shed light on the dynamics and triggering mechanisms of real earthquakes. The obtained results concerning the record-breaking avalanches in the OFC model can be used to constrain the effects of correlation and self-organization in seismogenic zones. Concerning real seismicity, the sequences of record-breaking earthquakes were extracted and analyzed for world-wide earthquakes [9,10]. The analysis of annual maximum earthquake magnitudes for the global earthquake catalog was also performed [47]. The results from those studies indicate that global earthquakes are independent and their magnitudes follow approximately the exponential distribution. The analysis of record-breaking earthquakes occurring on a single fault is problematic due to the limited records of existing catalogs. This makes it difficult to directly compare the correlations between avalanches we observe in the OFC model to real earthquake data.