Temporal complexity measure of reaction time series: Operational versus event time

Abstract Introduction Detrended fluctuation analysis (DFA) is a well‐established method to evaluate scaling indices of time series, which categorize the dynamics of complex systems. In the literature, DFA has been used to study the fluctuations of reaction time Y(n) time series, where n is the trial number. Methods Herein we propose treating each reaction time as a duration time that changes the representation from operational (trial number) time n to event (temporal) time t, or X(t). The DFA algorithm was then applied to the X(t) time series to evaluate scaling indices. The dataset analyzed is based on a Go–NoGo shooting task that was performed by 30 participants under low and high time‐stress conditions in each of six repeated sessions over a 3‐week period. Results This new perspective leads to quantitatively better results in (1) differentiating scaling indices between low versus high time‐stress conditions and (2) predicting task performance outcomes. Conclusion We show that by changing from operational time to event time, the DFA allows discrimination of time‐stress conditions and predicts performance outcomes.

TA B L E 1 The relation among the scaling indices , , and , which can be measured, respectively, from detrended fluctuation analysis (DFA), power spectral densities (PSD) S(f), and the waiting-time probability density function (PDF) ψ(t) of the time series

Scaled functions Parameter relations Parameter range
Waiting-time PDF ψ(t) ∝ t −μ μ = 3 − β 1 < μ < 3 Power spectrum Note: The value μ = 2 is the boundary between the underlying process having a finite (μ > 2) or an infinite (μ < 2) average waiting time and is also the point at which β = 1 where the process is that of 1/f-noise.
detrended fluctuation analysis (DFA) (2α−1 = ) and determines the well-known 1/f-noise when = 1. Simola et al. (2017) argued that the LRTC arises from critical dynamics and use this perspective to resolve a long-standing controversy concerning the source of LRTC in response time fluctuations. They concluded: "Our findings thus favor the hypothesis that LRTCs are caused [by] the system being in a critical state over the idea that these LRTCs would reflect long-memory dynamics" (p. 2).

Introducing crucial events (CEs)
Researchers have shown that individuals performing more difficult tasks requiring greater cognitive effort exhibit relatively smaller IPL indices of reaction time series, Y(n), than those observed for simpler tasks. For example, in a series of experiments, Kello et al. (2007) found smaller IPL indices for variable versus fixed trial intervals, random versus patterned cues, and unpreviewed versus previewed trials. Wijnants et al. (2009) showed that IPL scaling indices become greater (more clearly patterned 1/f-variability) as participants become more skilled in a precision-aiming motor task over blocks of practice. Correll (2008) found smaller IPL indices in participants who reported a high level of effort to avoid racial bias while responding to black versus white target stimuli under threat/no-threat conditions (see also Grigolini et al., 2009 for a theoretical explanation). Overall, the evidence suggests that greater IPL indices are indicative of greater levels of cognitive and behavioral performance and/or performance of simpler tasks.
In a similar Go-NoGo paradigm, Simola et al. (2017) showed that IPL scaling indices were inversely related to cognitive flexibility, as measured by errors of commission (i.e., by greater IPL indices associated with a lower number of errors). Although Simola et al. (2017) found support for criticality being the basis of LRTCs, they do not consider the alternative explanation that CEs might offer in a self-organizing critical system. CEs are defined to be a sequence of events separated by time intervals that are statistically independent of one another. The time intervals between independent events are generated by an IPL PDF with the IPL index that lies in the interval 1 < < 3. The relationship between complexity indices α, β, and μ is summarized in Table 1.
Finally, the statistical distribution of the CEs is renewal and consequently possesses a new kind of memory, one associated with renewal processes. This new kind of memory was discovered by Allegrini et al. (2002) and was given the apt name "memory beyond memory." It is this new kind of memory that is often confused with LRTC (West & Grigolini, 2021), the latter reflecting long-memory dynamics, whereas the former is a statistical property of the renewal CEs. The leading question therefore arises as to whether the true source of the IPL PDF of the Y(n) can be quantitatively determined. If the Y(n) time series is renewal, which is to say it is a CE time series, the above question could be answered because the time series would necessarily have the following properties: "The time series generated by complex processes are characterized by three regimes: the short-time regime, where the true complexity of the process is not yet perceived; an intermediate-time regime driven by the CEs; and a long-time regime where the process can be mistaken for an ordinary statistical process" (West & Grigolini, 2021, p. 16).

Empirical data
In the present work, the existence of IPL scaling indices is examined using data previously recorded during a Go- The shooting task required the participants to scan the environment for targets appearing over spatially distributed locations, detect and orient the weapon toward the target upon its appearance, identify whether the target is friend or foe, decide whether to shoot or not, and if the decision is to shoot, the participant must accurately aim and fire the weapon within the constrained time interval in which the target is exposed. As a task designed to assess inhibitory executive control mechanisms by evaluating errors of commission (shots fired at friendly targets), this shooting task variant represents a more ecologically valid, complex, and challenging Go-NoGo task than has been implemented in previous research. Typically, Go-NoGo tasks consist of stimuli appearing on a computer screen in a fixed location, and Go responses are simple button presses on a keyboard (Correll, 2008;Simola et al., 2017). We also have a relatively small number of trials in each condition in each session (n = 360) because the previous study was structured to limit the overall length of the experiment, which could introduce additional effects of boredom or fatigue. Thus, these data present a challenge for rigorous testing of theoretical models and analytical approaches for identifying power laws, IPLs, and scaling from shorter time series generated under more complicated task conditions than in previous applications.
Because of the comparatively small dataset resulting from our relatively short reaction time series, the previously established methods do not provide reasonable estimates of the IPL scaling index . The scaling index in previous studies is typically estimated by finding the slope of a linear approximation to the PSD function in log-log coordinates.
The PSD estimate is usually achieved by applying the fast Fourier transform to the reaction time trial series (1:N), which requires a relatively large number of stationary data points to obtain low variance and stable PSD estimates. Typically, researchers implement tasks consisting of 1000 or more trials to enable stable spectral estimates, and the evolution of longer periods to observe slow fluctuations, although as few as 200 trials have been employed (Correll, 2008). Because the IPL selection of interest in a PSD occurs at low frequency, we need to have long time series to obtain meaningful IPL behavior for at least one decade of frequency. As such, short time series create noisy PSD estimates, and consequently IPLs of uncertain slope, with a corresponding IPL index β.

New way to process reaction time series
DFA has been widely used to measure the scaling indices α of reaction time series (Delignieres et al., 2006;Simola et al., 2017;Wijnants et al., 2009). In previous studies, the signals analyzed by DFA have been the Y(n) time series measured at each trial n of the experiments. This represents a time series expressed in trial intervals referred to as operational time (Turalska & West, 2018). For simple tasks, IPL relations have been found using the IPL index for the PSD , or power law relations using DFA α (Correll, 2008;Gilden, 2001;Kello et al., 2007;Simola et al., 2017). However, we show herein that for reaction time series from a more realistic task having short trial sequences, neither the PSD method nor the DFA technique is able to extract a reliable measure of the scaling indices of the process.
Herein, we propose a new perspective by considering each reaction time as a time interval (t in ms) rather than a trial number (n) and filling each reaction time latency interval (i.e., event time) with the noise of fixed amplitude (of magnitude 1 and a random chosen sign ± for the interval; see Figure 2). This represents a latency or duration time series expressed as the new variable X(t), which is a function of temporal intervals (not trial numbers) referred to as chronological time, see Section 2.5 for details (Turalska & West, 2018). 1 This secondary time series, X(t) (Figure 2, panel c) represents the rigidity (constant noise) and flexibility (changing sign) of the process. By measuring the complexity of the X(t) time series, using DFA, we find a clear classification of scaling indices between low and high time-stress conditions. We also find clear trends between scaling indices and errors of commission.
Neither of these findings was observed for the traditional DFA analysis of Y(n) time series.
To show the connection of the scaling index α of the X(t) time series, measured using DFA, to the temporal complexity index μ, we introduce a simulation where the relation between the two scaling indices can be established and tested. The simulated Y(n) time series, with an IPL PDF having an IPL index μ, is generated using the well-known Manneville map (Manneville, 1980). Then, we created the corresponding simulated X(t) time series from the simulated Y(n) and measured the scaling index α via DFA. Consequently, we find the empirical relation between α and , = 4 − 2 , making it possible to connect α to the IPL index of the PSD (β = 3 − μ) and therefore adapt the technique to empirical reaction time datasets.
1 Hereafter we replace the term chronological time with event time because the time intervals between reaction time trials during the continuous task are ignored (i.e., time after a response has been made until the time of the next target onset), whereas the durations of the reaction time latencies are taken as events with temporal durations. As mentioned, see Section 2.5 for details of how X(t) is obtained from Y(n) by means of an analytic transformation.

F I G U R E 2
The schematics of changing data from operational time , each duration is filled with +1 or −1, assigned randomly (in this example, each latency interval is filled with +1, −1, +1, +1, and −1, respectively); (Panel d) the S(t) trajectory is the cumulative sum of X(t), notice that because of the random assignment of signs, there are many possible trajectories in event time X(t)'s. Thus, we use an ensemble average for our analysis to consider this variety, especially for short time series. In panels (b) and (c), t represents accumulated durations of reaction times. Note that the trajectory X(t) hosts both CEs and non-CEs from the empirical trials. The trajectory is then analyzed using detrended fluctuation analysis (DFA) to determine the scaling index α, which is not the same index as that obtained by applying DFA to Y(n). This difference is discussed in the following section.

Task design
The Go-NoGo task was implemented in virtual reality using HTC Vive (https://www.vive.com/us/). In each session, the participants com- and right) and exposed at variable onset intervals (1000 ± 500 ms over a Gaussian-distributed range of 100 ms increments) for various target exposure durations (see Figure 1). The probability of targets (enemy; red) to nontargets (friendly; green) was 0.90/0.10, respectively (i.e., 324 enemy targets and 36 friendly targets were presented) to induce a prepotent response bias (see Kerick et al., 2023 for more details).
The participants were instructed to "shoot enemy targets as quickly and accurately as possible, while refraining from shooting at friendly targets." Time-stress conditions were individualized based on a pretesting performance thresholding procedure to account for individual differences in participants' ability to perform the task.
This was done by empirically determining target exposure durations corresponding to the 50th (High time-stress) and 90th (Low timestress) percentile hit-rates in response to 100 enemy targets using psychophysical methods (method of limits; Ehrenstein & Ehrenstein, 1999). The mean (SD) target exposure time in the low time-stress condition was 842.78 ms (216.65 ms) and 547.96 ms (223.20 ms) in the high. Figure 1 provides example reaction time series from a representative subject in low and high time-stress conditions in one session.
Herein we analyze the scaling index of reaction times series and examine their relation to the performance measure errors of commission, defined as shots fired at friendly targets (failure to inhibit incorrect responses).

Data analysis
Across

DFA analysis of Y(n) time series
For each subject under each condition and time scale, we computed the percentage of errors of commission and applied DFA analyses to Y(n) time series. DFA is applied to reveal long-range correlations in time series (Peng et al., 1995). Here we explain the steps for DFA on Y(n): (1) integrating Y(n) over n trails to obtain S(n); (2) dividing S(n) into windows of length w; (3) deriving least squares line to fit the dataset in each window, S w (n); (4) calculating the root mean square amplitude in each window w, F(w) of Equation (1); (5) average the F(w)s; (6) repeat 1 through 5 for different ws, and (7) evaluate the scaling α as the slope of F(w) versus w in a log-log plot: (1)

DFA on simulated reaction times with different temporal complexity μ
To introduce the method used herein, and the relating temporal complexity to scaling measured by DFA, we used simulated data, temporally complex, using turbulence intermittency as prescribed by Mannevil-lleManneville (1980). To make the numerical generation of the events faster, we adopted an idealized version of this map in the following way (Buiatti et al., 1999): Each reaction time can be calculated by transforming the sequence y i of random numbers uniformly distributed in the interval (0, 1) into the temporal sequence τ i : where μ (μ > 1) is the temporal scaling index of the generated time series. The τ s are generated by a hyperbolic PDF ψ(τ) defined by which is properly normalized. The time constant T defines the time necessary to turn microscopic dynamics into a process with an evident IPL index μ. For 1 < μ < 3, this PDF represents a complex system, which is a renewal process of CEs. For μ > 3, the dynamics falls into the region of normal statistics (Annunziato & Grigolini, 2000) and is no longer complex. We take the set of τ i generated by the Manneville map, with a given scaling index μ, as the simulated reaction times Y(n) and create the corresponding X(t) time series. The DFA was then used to process the resulting X(t) time series generated using the prescription in

Relation between α and 1/f-variability
Here we connect the scaling index α measured via DFA processing of the X(t) time series to the power-law index of f − variability. It is known that the series of events created by the Manneville map with scaling index μ has an IPL PDF (Lukovic & Grigolini, 2008): Substituting the relation between index given by Equation (4) into this expression, we obtain So, only α = 1 corresponds to a perfect 1/f-noise.

Statistical analysis
Wilcoxon matched-pairs signed rank tests were applied to determine   Table 2 provides median and quantile (0.25 0.75) values for each dependent variable under each time-stress condition. Figure 5 shows differences of scaling indices between low and high time-stress conditions for the DFA analysis of Y(n) and X(t) time series. The DFA analysis of the X(t) time series was able to categorize the scaling indices of the two conditions, whereas Y(n) was not.

Differences between time-stress conditions for errors of commission and DFA scaling indices
Next, we evaluated whether shuffling the reaction time series, Y(n), prior to DFA analysis destroys the long-term correlations inherent in the intact reaction time series and lowers scaling indices. Figure 6 F I G U R E 6 The left panel shows the results of applying detrended fluctuation analysis (DFA) to the Y(n) (black squares) and to shuffled Y(n) (red dots). The right panel depicts the results of applying DFA to X(t) created from Y(n) (black squares) and to X(t) created form the shuffled Y(n) (red dots). compares the DFA applied to randomly shuffled Y(n) and their corresponding X(t) time series. After shuffling, the DFA analysis of Y(n) time series yielded scaling indices close to α = .5 (randomness), whereas in the case of DFA analysis of the X(t) time series, created from the shuffled Y(n), does not significantly change the scaling indices.

Regression analyses of DFA scaling indices for predicting errors of commission
Regression analyses revealed that scaling indices derived from DFA analysis of Y(n) time series was not a significant predictor of errors of commission in either low or high time-stress conditions for either individual sessions or all sessions appended (all p values >.05; see Table 3).
However, scaling indices derived from the DFA analysis of the X(t) time series were all highly significant (all p values <.001; see Table 4). Higher scaling indices were associated with lower errors of commission. Further, the associations were stronger for data in the high versus low time-stress condition.

DISCUSSION
We developed and applied a new approach to measuring the scaling index of reaction time series data using DFA, from a short, complex cognitive-motor decision-making task. In this method, we consider  Although Simola et al. (2017) observed IPL scaling of Y(n) time series from a Go-NoGo task using autocorrelation functions, PSD, and DFA, their task consisted of 1000 equally spaced trials and was presented on a computer monitor requiring keyboard responses. These two methodological differences (lower trial numbers and greater task complexity) may account for why we did not observe the hidden IPL scaling of Y(n) time series in the present study. Further, the new method of DFA processing of the transformed time series X(t) was able to discriminate between low and high time-stress conditions, whereas the DFA analysis of Y(n) time series was not. Our finding by applying DFA to X(t) time series is also consistent with previous research, which has observed higher scaling indices in Y(n) time series associated with lower task demand or decreased task complexity with longer trial series (>1000; Kello et al., 2007) (see also Delignieres et al., 2006).
The stronger negative relation between scaling indices and errors of commission in the high time-stress condition suggests that under higher time demand conditions (less time to complete tasks), it is beneficial for the system to shift from a more ordered, predictable state to a more disordered, unpredictable state. However, a shift too far in the direction toward disorder results in the deterioration of inhibitory control. Lower scaling indices in the high versus low time-stress condition suggest that behavior is more disordered or unpredictable (more random) in the high time-stress condition, an effect that appears to be environmentally coupled to the time constraints imposed by the task.
Together, these two findings support previous research suggesting that system complexity is associated with greater degrees of freedom, thus facilitating greater flexibility and adaptability of the system to internal and external perturbations or task demand conditions (Kello et al., 2007;2010;Van Orden et al., 2003, 2005Wijnants et al., 2009).
Traditional DFA provides a method to quantify long-range correlations in time series and to index repeating patterns over different time scales (Peng et al., 1995). As our results reveal, randomly shuffling The results of shuffling analysis show that the X(t) time series, unlike the traditionally used Y(n) time series, is insensitive to the order of the reaction times. In our experiment, the time intervals within sessions were randomly distributed and there were longer time intervals between sessions. Due to this property of the method, we could concatenate trials even when temporal distance between trials differed significantly (e.g., between sessions). This made the method suitable for even more realistic experimental and real-world conditions where intervals between responses to stimuli may vary significantly. In other words, the statistics of X(t) are renewal, and the time interval between events are statistically independent of one another. Shuffling the Y(n) will not change the statistics of its corresponding X(t), so they will again be an IPL with index μ.
Now we consider a second time series to which we apply DFA and get a scaling index λ. But we know when we shuffle the time series, we are going to obtain a different scaling index, say α, because the index λ we obtained had to do with long-range correlations that are destroyed when we shuffle. If the second time series (before shuffling) is fractional Gaussian noise then after shuffling, we would obtain α = .5. However, if the time series is mixed, we could obtain α = μ so that the shuffling acts like a filter, and the deviation of α from 0.5 is that part of the mixed time series that is renewal.
We note that the method of subordination to ordinary diffusion proposed by Sokolov (2000) to illustrate the popular theory of continuous time random walk (Scher & Montroll, 1975) is widely adopted to generate the events that are the source of temporal complexity revealed by the statistical analysis of experimental data. For instance, Turalska and West (2018) used this theoretical interpretation to describe the temporal complexity of individuals of a decision-making social system.
Subordination makes it possible to convert the fluctuations of data that would not generate a departure from ordinary diffusion into rare events yielding temporal complexity. Our new method is based on converting the intensity (latency) of reaction times into the time intervals between consecutive events. The time interval between consecutive events is filled with either +1 or −1, tossing a fair coin (see panel c of Figure 2). In Bologna (2020), this approach was used to illustrate the anomalous diffusion properties of events, with a waiting-time PDF known to be proportional to 1/τμ with τ denoting the time interval between consecutive events. However, in this paper, the number of reaction times available to us is not large enough to allow us to determine with precision the power law index μ from the distribution of the reaction times. We have directly converted the reaction times in the chronological order into a diffusion process without finding a deviation from ordinary scaling. Quite surprisingly we found that the adoption of the current model, used for the first time to the best of our knowledge, to detect the complexity of the experimental fluctuations allows us to evaluate with precision the parameter μ.

CONCLUSIONS
In this study, we introduce a novel conceptual and analytical framework for estimating the complexity of behavioral time series data based on CEs. Analyses of the X(t) time series, using DFA, we found a clear classification of IPL scaling indices between low and high time-stress conditions. We also found clear relations between IPL scaling indices and errors of commission. Neither of these findings was observed for the traditional DFA analysis of Y(n) time series. This new approach was borne out of necessity due to the relatively short reaction time series (360 trials over 8-10 min task durations), which has strong implications for future research in cognitive science and neuroscience experiments, where task durations may be relatively short and where data are collected intermittently within and across recording sessions.

AUTHOR CONTRIBUTIONS
Scott E. Kerick conceived and conducted the original experiment; Korosh Mahmoodi, Paolo Grigolini, and Bruce J. West conceived and applied the novel DFA analysis. All authors contributed to analyses of results and writing and reviewing the manuscript.