Temporal and spatiotemporal perturbations in paced finger tapping point to a common mechanism for the processing of time errors

Paced finger tapping is a sensorimotor synchronization task where a subject is instructed to keep pace with an external metronome, as when following along with the beat of music, and the time differences (asynchronies) between each stimulus and its response are recorded. The usual way to study the underlying error correction mechanism is to make unexpected temporal perturbations to the stimuli sequence and then let the subject recover average synchronization. A critical but overlooked issue in traditional temporal perturbations, however, is that at the moment of perturbation two things change: both the stimuli period (a parameter) and the asynchrony (a variable). In terms of experimental manipulation, it would be desirable to have separate, independent control of parameter and variable values. In this work we perform paced finger tapping experiments combining simple temporal perturbations (tempo step change) and spatial perturbations with temporal effect (raised or lowered point of contact). In this way we decouple the parameter-and-variable confounding of traditional temporal perturbations and perform novel perturbations where either the parameter only changes or the variable only changes. Our results show nonlinear features like asymmetry and are compatible with the idea of a common mechanism for the correction of all types of asynchronies. We suggest taking this confounding into account when analyzing perturbations of any kind in finger tapping tasks but also in other areas of sensorimotor synchronization, like music performance experiments and paced walking in gait coordination studies.


Introduction
Sensorimotor synchronization (SMS), the mainly specifically human ability to keep synchrony with an external periodic metronome, is the basis of all music and dance [ReppSu2013, Patel2009, Schachner2009, Zatorre2007]. It is a spontaneous behavior found in humans and, despite its simplicity, the mechanisms underlying it remain mostly unknown; the processing of temporal information is an open area of research in neuroscience and how time is represented and manipulated in the brain is still one of the most elusive concepts, particularly in the millisecond timing range (hundreds of milliseconds) [IvrySpencer2004, IvrySchlerf2008, BuonomanoLaje2010, Grondin2010, Paton2018].
The most used experimental paradigm in SMS is auditorily paced finger tapping, a task where a subject is instructed to move a finger (tap) in synchrony with a periodic stimulus sequence of short tones (beep) as in keeping pace with music ( Figure fig.tapping a). The main observable, both for behavior description [Chen1997] and for mathematical modeling [ReppSu2013], is the asynchrony e n = R n -U n , that is the difference between the occurrence time of each response R n and the occurrence time of its corresponding stimulus U n . It is simple to demonstrate the existence of an underlying error correction mechanism---if the stimuli sequence is suddenly silenced while the subject keeps tapping, the "virtual asynchronies" between responses and extrapolated beeps become very large very soon drifting away from average synchrony. In order to understand the inner workings of such error correction mechanism, temporal perturbations to the stimuli sequence are usually performed and then study the resynchronization [Bavassi2013]. Traditional temporal perturbations consist in for example an abrupt change in the period of the sequence from a given beep on, which is known as a tempo step change ( Figure fig.tapping b). There is, however, an overlooked issue with traditional perturbations. Let's consider a simple perturbation like the mentioned tempo step change. At step n when the perturbation occurs, both the parameter value (stimuli period) and the variable value (asynchrony) change, because the stimulus occurrence time U n changes and thus e n changes instantaneously and arbitrarily---importantly, without involving any dynamics of the underlying system ( Figure  fig.tapping b). The critical difference between period and asynchrony is that the former has no dynamics by itself and its value is up to the experimenter, while the latter is subject to the dynamics of the system under study; they might have, for instance, different neural correlates. Note that this issue appears in any traditional perturbation that modifies any stimulus occurrence in the sequence, like for instance the above-named tempo step changes [ReppKeller2004, Thaut1998], but also sigmoidal changes [SchulzeCordesVorberg2005, PecenkaKeller2009, vanderSteen2015], linear changes [Loehr2011], quadratic changes [Pecenka2013], sinusoidal changes [Michon1967, PecenkaKeller2011], random changes [Michon1967, HaryMoore1987, MadisonMerker2004, ReppPenel2004], phase shifts [Repp2001phaseshifts], event onset shifts [Repp2002phasecorrection, Repp2002automaticity, Flach2005], and adaptively timed sequences [ReppKeller2008, Repp2012, Mills2015]. Up to our knowledge, this is the first work to note this issue.
Beyond the finger tapping literature, research in other areas within SMS are quite likely impacted by this parameter-variable confounding. Two areas are easily spotted. First, research on music performance and perception also relies on perturbation experiments and ecological experimental conditions with natural variability in the stimuli, like accelerando, ritardando, and natural expression [Rankin2009], groove and microtiming [Madison2011], musicians following a conductor [Luck2006], phase shift and tempo perturbations [Levitin2018, Zanto2005, LargeKelso2002], auditory delayed feedback [Pfordresher2006], and even on theoretical works about oscillator models of synchronization to a musical beat [LargePalmer2002]. Second, research on gait coordination also makes use of experimental designs that might be confounding parameter and variable perturbations, like phase shift perturbations [Roerdink2009, Pelton2010], interpersonal coordination during side-by-side walking [vanUlzen2010], and varying belt speeds and directions in treadmill walking [Choi2007].
An analogy may help convey our point and our proposal. Let's consider the circadian clock, a well-characterized genetic oscillator that in mammals is in the suprachiasmatic nucleus of the hypothalamus. The circadian clock displays an autonomous oscillatory activity with a period of around 24 hs capable of synchronizing to the daily cycle mainly through the light/dark stimulus [Golombek2010, Laje2018]. In a minimalistic description, the system works as follows: a gene produces a protein (BMAL/CLOCK) that activates a second gene, which in turn produces a second protein (PER/CRY) that inactivates the first gene. In this system the period of the external light/dark stimulus is a parameter, while the concentrations of CLOCK/BMAL and PER/CRY proteins are the variables whose time evolution is set by the dynamics of the system. Both the parameter and the variables may be experimentally manipulated independently: the parameter value can be modified by changing the period of the external light/dark stimulus (called T-cycle [Golombek 2013), and the value of the variables can be modified for instance by applying an acute dose of inhibitor PF-670462 that modifies the concentration of PER [Kim2013].
It would be desirable in paced finger tapping to manipulate parameter and variable values independently. In this work we perform paced finger tapping experiments with simple temporal perturbations (traditional tempo step changes), and also novel spatial perturbations with temporal effect (raised or lowered point of contact) and combined perturbations with the aim of probing the error correction mechanism. This allows us to decouple the effect of traditional temporal perturbations and perform novel perturbations where the parameter only changes (a change in stimuli period without a change in asynchrony) and where the variable only changes (a change in asynchrony without a change in period). Our results show nonlinear effects even when the perturbations are 10% of the stimulus period and are compatible with the idea that the origin of the asynchrony doesn't influence the subsequent resynchronization. The issue of whether motor timing and sensory timing depend on the same neural circuits is still open [BuonomanoLaje2010], but our results suggest that asynchronies produced by purely temporal perturbations and those produced by spatiotemporal perturbations are processed by a common error correction mechanism.

Results
We performed an auditorily paced finger tapping experiment with perturbations, where the subject is instructed to keep in synchrony at his/her best and keep tapping to resynchronize in case a perturbation appears at a random beep. Perturbations can be any of the following: a) simple temporal perturbations "T" (traditional tempo step change perturbations by an amount +/-DeltaT); b) simple spatial perturbations "S" that have temporal effect, achieved by raising (+) or lowering (-) the point of contact where the subject has to tap and thus advancing or delaying the time of contact; c) combined simultaneous perturbations "ST". Perturbations were classified according to two criteria: by size (small and large perturbations), and by type (simple temporal +/-T; simple spatial +/-S; combined analogous +S+T and -S-T; and combined opposite +S-T and -S+T). Note that in the simple spatial +/-S perturbations the stimuli period doesn't change and thus they are variable-only perturbations; in the combined opposite +S-T and -S+T the stimuli period changes but the asynchrony doesn't because both components compensate each other on average at the perturbation beep, and thus they are parameter-only perturbations. See Methods for a detailed description and rationale.  One of the effects of a tempo step change perturbation in paced finger tapping is a shift in baseline relative to its pre-perturbation value [Bavassi2013]. Such baseline, which represents the state of the system in stationary conditions, is the extensively studied "negative mean asynchrony" or NMA [Repp2005]. In this work, on the contrary, we are interested in analyzing the system dynamics right after being perturbed and how it converges to its post-perturbation value [Bavassi2013]. In order to do that, from here on we choose to plot every time series after subtracting its post-perturbation baseline such that all time series converge to zero after perturbation (see Methods).

The response to simple perturbations is asymmetric
In Figure fig.simples a and b we show the averaged time series for perturbations +T and -T. Large perturbations produce asymmetric responses--the response to the positive perturbation converges to zero more rapidly than the negative one, in accordance with previous reports for this type of perturbation [Thaut1998, Bavassi2013]. Small perturbations, on the other hand, do not show asymmetry. Both responses are consistent with a nonlinear underlying system, where nonlinear behavior like asymmetry is evident only when the perturbation size is large enough [Loehr2011].
Spatial perturbations +S and -S, on the contrary, produce responses that are mostly symmetric (Figure fig.simples c and d). It is worth noting that the asynchrony values produced by +/-S perturbations at the perturbation beep are relatively small and similar in size to those produced by the small +/-T perturbations, which also have symmetric responses ( Figure  fig.simples a). Spatial perturbations +/-S are novel in that they produce a change in the value of the variable (asynchrony e n ) without changing the stimulus period. Large temporal perturbations produce asymmetric responses, and the grey area indicates steps where the asymmetry is significative (p=0.0025 each). All traces are mean +/standard error across subjects.

The response to combined perturbations is asymmetric
In this section we analyze the degree of asymmetry of the combined perturbations, and for that we group them in "analogous" (+S+T and -S-T, because their components S and T individually would produce asynchronies of the same sign) and "opposite" (+S-T and -S+T, whose components would individually produce asynchronies of opposite signs; see Methods). We note that the analogous perturbations produce asynchronies with values at the perturbation beep (n=0) greater than the temporal perturbations, while the opposite perturbations produce values lesser than the temporal perturbations (Wilcoxon sign-rank test, significant after Bonferroni correction; simple temporal vs. analogous p=0.0134; simple temporal vs. opposite p=0.00073; analogous vs. opposite p=0.00048). The result in this section may be understood within the same framework that the result from section 3.1 : it is the behavior of a nonlinear underlying system whose asymmetric response is only evident when the asynchrony values are large enough, and whose degree of asymmetry is larger for greater asynchronies. Response to analogous perturbations -S-T and +S+T (whose components S and T would individually produce asynchronies of the same sign) and degree of asymmetry (ASYM). (b) and (c) Response to opposite -S+T and +S-T (whose components S and T would individually produce asynchronies of opposite signs) and degree of asymmetry (ASYM). Left: small perturbations; Right: large perturbations. When compared to the large temporal perturbations of the previous section, analogous perturbations produce more asymmetric responses while opposite perturbations produce less asymmetric responses. Mean +/-standard error across subjects. The grey area indicates steps when the asymmetry is statistically significant (p=0.022, 0.0025, 0.0025, 0.0067, 0.022, respectively).

Simple perturbations are additive
One of the aims for proposing the combined perturbations (+/-)S(+/-)T was to compare every combination with the simple perturbations they consist of. For example, and to show the notation, we would like to compare between the actual combined perturbation +S-T and the sum of the simple perturbations (+S)+(-T).

A coarse-grained conceptual model for the error correction mechanism
The results shown in the previous sections allow us to distinguish between various ways of how temporal information is processed in this task and so we will consider a very simple, coarse-grained, conceptual model for the error correction mechanism. The input is the temporal information coming from performance monitoring (time perception) and that could enter the system through different paths depending on the condition--for instance, temporal information from a simple temporal perturbation T and from a simple spatial perturbation S could travel different paths. On the other end, the output is the result of processing such temporal information and drives the motor output. The central part is the time processing itself where hypothetical processes like time comparison, decision making and time prediction and production take place [Grondin2010, Barne2018, Paton2018, IvrySpencer2004, IvrySchlerf2008, Bueti2011, Coull2011], and where the estimation of the appropriate asynchrony correction is made. It should be noted that the parts of this model are not necessarily associated with individual brain regions. The central part of time processing might include sensory components (at its beginning) and motor components (at its end), besides the hypothesized processes for time comparison, decision making, etc. That is, part of the processing itself might take place in sensory and/or motor cortices. We propose the following hypothesis: 1. The processing of temporal information is nonlinear. In particular: a. the response to symmetric perturbations is not symmetric; b. the response to a combined perturbation is not equal to the sum of responses of the individual perturbations. 2. The processing of temporal information coming from temporal perturbations and from spatial perturbations is performed by a common mechanism. Our hypotheses are represented in four versions of the conceptual model as combinations of 1) linear vs. nonlinear processing, and 2) common vs. separate processing. The four combinations are displayed in Figure fig.models . A straightforward result is the observed asymmetric responses, sections 3.1 and 3.2, that validate hypothesis 1.a and thus the linear models (c) and (d) must be discarded. It is worth noting that the symmetric responses to the spatial perturbations S ( Figure fig.simples b) don't necessarily imply that the processing of asynchronies coming from spatial perturbations is linear and consequently discard hypothesis 1.a--as we noted above, spatial perturbations S in our experiments produced small asynchronies in comparison to the ones produced by the large temporal perturbations T and thus they might be insufficient to elicit typical nonlinear behaviors. Now consider hypotheses 1.b and 2. An unwarranted conclusion from the additivity results presented in Figure fig.add would be that T perturbations and S perturbations are processed either by a common linear mechanism or by separate nonlinear mechanisms and then added up. However, in front of the discussion of the previous paragraph, if the S perturbations are indeed too small to elicit nonlinear behavior then we would rather expect they might show as additive. That is, the most parsimonious conclusion is that we cannot discard a single nonlinear mechanism and thus we cannot decide whether hypotheses 1.b and 2 are true or false. The two models backed up by our results are then (a) and (b), and to distinguish between them we would have to perform perturbations where the asynchronies produced by S perturbations were comparable in size to those produced by the large T perturbations. However, the remarkably consistent behavior across perturbation sizes, perturbation signs, and perturbation types is noteworthy, suggesting there might be a single underlying mechanism in charge of correcting asynchronies regardless of their origin.

Experimental phase space
In a previous work, we showed experimental and theoretical evidence suggesting that the correction mechanism underlying resynchronization after a temporal perturbation in paced finger tapping is bidimensional and nonlinear [Bavassi et al 2013]. If one of the variables of such a system is the observable e n , how could we obtain information about the other variable? The embedding technique [Gilmore1998] allows us to get information about the geometric organization of the underlying system's trajectories by analyzing the time series of a single observable, in this case e n .
To illustrate the procedure, in Figure

Trajectories of the combined perturbations
We've shown that the simple perturbations S produce a novel effect: a change in asynchrony without a change in the stimulus period. On the other hand, the combined opposite perturbations (+S-T and -S+T) are also novel: they produce a change in the period without changing on average the expected asynchrony value (see Methods and Supplementary Figure fig.schematic ). For instance, let's consider a traditional simple temporal perturbation -T representing a tempo step decrease (a change in stimulus period by an amount -DeltaT). It normally produces an increase in asynchrony by an amount -DeltaT at the perturbation step. By combining -T with +S to get a +S-T perturbation, a decrease in stimulus period is produced due to the -T component but with a lesser asynchrony value because of the temporal compensation produced by the +S component. Conversely, the combined analog perturbations (+S+T and -S-T) make the asynchrony value greater (in absolute value) than the one produced by a T perturbation alone.
In Figure fig.embedST a we compare the time series of the traditional perturbations -T and +T (size 50 ms only) and the corresponding combined opposite and analogous perturbations. All three types of perturbations (traditional, opposite, analogous) produce similar time evolutions except for their initial asynchrony value (n=0). In the phase space, this can be seen as the three trajectories on the same side making similar paths after a quick convergence from differing initial points. That is, the responses to perturbations that have a change in period are similar despite different initial asynchrony values.

All perturbations: new system states
In this section we continue analyzing the experimental phase space shown above by including the information from the novel perturbations presented in this work. The combined perturbations, as well as the already discussed simple spatial perturbations S, allow us to decouple the effect of traditional temporal perturbations and access as yet unexplored system states. We will normalize the trajectories from all perturbations so we can make evident their positions in phase space relative to one another. . Shaded areas represent the regions in phase space (i.e. system states) explored by the traditional tempo step change perturbations. In Figure fig.areas b we include the rest of the perturbations (simple spatial +/-S, combined opposite +/-S-/+T, and combined analogous +/-S+/-T, after normalization and non-dimensionalization). It is evident that the regions in phase space explored by the novel perturbations overlap with the traditional +/-T but only partially, and thus they allow us to probe the system in novel ways and access system states as yet unexplored as a result of decoupling the experimental manipulation of period and asynchrony.

Conclusions
Experimental work on sensorimotor synchronization and, in particular, both paced and unpaced finger tapping date back to the end of the XIX century [Stevens1886], while the first mathematical models were published 50 years ago, e.g. [Michon1967]. Much work has been done in this area of research to advance our knowledge about the underlying mechanism responsible for achieving average synchrony [Repp2005, Repp2013]. In particular, perturbations to the sequence period are a usual way to probe the system's inner workings.
As we described in this work, these perturbations are in fact confounded parameter and variable manipulations and, up to our knowledge, no published work has addressed this issue.
Our observation applies to any perturbation to the sequence period: local step changes, phase shifts, and event onset shifts, and global changes like accelerando or ritardando, etc. A different, less usual way of perturbing the system is also worth discussing: temporally delayed or advanced auditory feedback from the taps [AscherPrinz1997, Wing1977, PfordDalla2011, MatesAscher2000]. Note that the period of the sequence is not changed in this kind of perturbation. However, it introduces a dissociation between auditory feedback and proprioceptive and tactile feedback from the taps and thus its relationship to the (perturbed) asynchrony value remains to be elucidated.
We showed that an experimental phase space can be reconstructed from the averaged time series via embedding and that all trajectories (corresponding to either simple or combined perturbations of any size and sign) organize in a geometrically remarkable way. This supports the notion that the underlying system in charge of the error correction can be considered as a single mechanism [Bavassi2013], as opposed to many different mechanisms that are turned on or off depending on the magnitude, sign, and type of perturbation, a usual---and frequently implicit---assumption when perturbation data are considered [Thaut1998]. Note that this doesn't imply that the details of such a mechanism are simple---in fact, we expect that the hypothetical processes within such a mechanism like time perception, comparison, and production interact in complex ways, particularly in this task where processing stimuli and producing timed responses is ongoing and the processes related to step n are likely overlapped in time with those from the previous step n -1 [Bavassi2017]. Yet, according to our results, it appears that the response of the system is consistent regardless of whether the asynchrony has a purely temporal or spatiotemporal origin, and of the size and sign of the perturbation.
Research about how temporal information is processed in the brain has grown very rapidly in the last decade, particularly research related to the neural underlying mechanisms and involved brain regions [Paton2018]. Many important issues remain open, like whether timing is centralized and dedicated or distributed or intrinsic; whether it reflects properties of individual neurons or is the emergent behavior of neural networks; whether sensory and motor timing rely on the same circuitry or not; and so on. Our results show that traditional perturbations are confounded parameter-variable manipulations, and that a broader set of perturbations (to either the parameter, the variable, or both) leave the system in different points of the phase space. This suggests revisiting the usual assumptions about perturbations and interpreting results with this distinction in mind. In the music domain, for instance, our results open the possibility that different brain regions and processes might be recruited in front of a perturbation to either the asynchrony or the tempo [Zatorre2007], and that period and asynchrony might have different neural representations or correlates. Whether the known associations between functions and regions would split after solving the confounding is yet to be elucidated.

Task and subjects
The task was an auditorily paced finger tapping that may, at some random point in the series, present a perturbation. The subject was instructed to keep in synchrony at his/her best by using his/her index finger (holding the wrist in place during the experiment) and keep tapping to resynchronize in case a perturbation appears. Each trial was either isochronous (constant stimulus period and fixed point of contact) or perturbed. Perturbations were any of the following three types: a) simple temporal perturbations "T" (traditional tempo step change perturbations); b) simple spatial perturbations "S" that had a temporal effect (see below); c) combined simultaneous perturbations "ST". Subjects were volunteers; one subject was not able to complete the task, and three subjects were excluded according to the criteria below. The final number of subjects was N=30 (ages 19-40, mean 29.1, 13 female, 28 right-handed, all subjects used their dominant hand).

Perturbations
We performed three types of perturbations: a) Simple temporal perturbations +/-T. These perturbations were the traditional tempo step changes where the stimuli period changes once by an amount +/-DeltaT. b) Simple spatial perturbations +/-S. These novel perturbations consisted in raising (+) or lowering (-) the contact point between finger and sensor, which had the temporal effect of advancing or delaying the time of contact (when the sensor was in its higher or lower position, respectively) and thus made the asynchrony at the perturbed step more negative or more positive, respectively. c) Combined perturbations ST, where simultaneous (i.e. at the same step n) S and T perturbations take place: +S+T, -S+T, +S-T, -S-T. In any case the perturbation was unexpected (it occurred at a randomly chosen step).

Setup
The setup consisted of two interconnected microcontrollers (Arduino Mega) communicating via the Wire library. The "master" microcontroller received instructions from the control script (MATLAB, The MathWorks, Inc) and started a trial (it received trial parameters, generated the sequence of auditory stimuli, sent instructions to the slave microcontroller, registered taps, sent trial data back to the control script). The master microcontroller had a custom-designed shield [Bavassi2017] for interfacing and signal conditioning (sum and amplification of audio signals, virtual ground, voltage divider for the force sensor, etc.). Every time a tap was detected, auditory feedback was sent as a brief tone. Stimuli and feedback tones were 50-ms duration square waves of 440 Hz (A4) and 587 Hz (around D5) respectively. Visual feedback from the subject's hand was avoided by means of a blocking screen. The "slave" microcontroller was in charge of producing the spatial perturbations S whenever the master sent the instruction. The slave microcontroller drove a small platform by means of a servo motor (Savox SC-1258TG Super Speed Titanium Gear Standard Digital) and a rotating arm, displacing the platform upwards or downwards along a linear ball bearing. The platform had a force-sensitive resistor (FSR 406, Interlink Electronics) attached on top of it to detect the subject's response. Sound was played diotically through Sennheiser HD419 headphones, and subjects adjusted sound volume to a comfortable level.

Experimental design
The combination of temporal perturbations T, spatial perturbations S, and perturbation direction (positive or negative) resulted in 9 experimental conditions: isochronous (no perturbation), 4 simple perturbations (temporal +T and -T and spatial +S and -S), and 4 combined perturbations (+S+T, -S+T, +S-T, -S-T). Each subject participated in a single session with three phases: Introduction, Calibration, Test. Each phase ended when all trials were successfully completed. Each trial consisted of a sequence of 30 auditory stimuli with a period (interstimulus interval) of 500 ms until a perturbation occurred randomly between the 15th and the 20th stimuli. In the temporal perturbations the period of the sequence was either increased (positive, +T) or decreased (negative, -T). In the spatial perturbations the vertical position of the platform was either moved upwards (positive, +S) or downwards (negative, -S) from its resting position by a distance that depended on the subject (see calibration below) and the new position was kept until the end of the trial. In the combined perturbations both T and S manipulations were performed simultaneously. The introduction phase was an exposure to the task and the 9 experimental conditions (one trial per condition, 9 trials total) to allow the participant to familiarize with the experimental task. The calibration phase consisted of 8 trials with simple spatial perturbations +/-S with the most extreme positions of the platform (+/-1.5 cm, 16 trials total), so we can estimate the longest temporal effect of the spatial perturbations for each subject. The test phase consisted of 8 trials from each experimental condition (72 trials total), in three blocks separated by short breaks to prevent the participant from becoming tired. Trials were presented in random order within each phase.

Calibration
The temporal effect of raising or lowering the platform is subject-dependent since every subject moves the finger at its own speed and, thus, a given spatial displacement of the platform is translated to different times for the finger to reach the new contact point. In the calibration phase we estimated for each subject the average time error produced by the unexpected movement of the platform from its resting position to its two most extreme positions (highest 1.5 cm, lowest -1.5 cm).

Experiments
We performed two experiments: Experiment 1 Matched temporal and spatial perturbations (i.e., the temporal effect of the spatial perturbation is equal to the assigned temporal perturbation size). Subjects were assigned to either one of the following groups, depending on the recorded average time error during the calibration phase:

Experiment 2
Unmatched temporal and spatial perturbations. Subjects with a large temporal effect in spatial perturbations were very infrequent, so in order to test larger temporal perturbation magnitudes we created a fourth category with +/-50-ms temporal perturbations and +/-15-ms spatial perturbations. Number of subjects: N 4 =8.

Exclusion criteria
A trial was correct if the subject began tapping before the 6th stimulus and if the subject didn't miss any response from there on. Each incorrect trial was repeated until completed successfully. Asynchrony outliers: after finishing the experiment, any trial with very large asynchronies (greater than 145 ms in absolute value, with respect to the pre-perturbation baseline) was discarded. Standard deviation outliers: the distribution of the standard deviation of pre-perturbation asynchronies for each subject was computed, and any trial with a standard deviation greater than 1.5 times the IQR was discarded. After discarding asynchrony and deviation outliers, subjects with less than 4 correct trials for every condition were removed from further analysis (3 subjects, 3 different conditions).

Data pre-processing
In order to average across trials and subjects, we aligned all trials at the perturbation step (renamed as n=0) and defined the analysis range from n=-10 to n=14 such that within that range all subjects have all responses. We define the trial pre-perturbation baseline as the average of all asynchronies between n=-7 and n=-1, so we can avoid adaptation effects at the beginning of the trial. The trial post-perturbation baseline is the average asynchrony between n=9 and n=14.

According to their size
In order to increase the statistical power, a posteriori we grouped the data in two categories according to the perturbation size: small and large perturbations. Small perturbations include 15-ms and 30-ms perturbation sizes (N=N 1 +N 2 =17 subjects), while large perturbations include 45-ms and 50-ms sizes (N=N 3 +N 4 =13 subjects).
According to their type According to their origin and expected effect, perturbations can be classified as: Simple: temporal +/-T and spatial +/-S. Combined: simultaneously changing the period and raising/lowering the point of contact. If the temporal effect of both components is in the same direction (i.e., the asynchrony at the perturbation step either increases or decreases) the perturbation is Analogous (-S-T and +S+T); if the effect goes in opposite directions (for instance one component increases the asynchrony while the other decreases it) the perturbation is Opposite (-S+T and +S-T). See Supplementary Figure fig.schematic .

Averaging
Time series in Figures fig.todas through fig.embedST  In order to compute the p-values we generated null distributions of ASYM time series and DIFF time series. We illustrate the procedure with an example. To generate the null distribution of ASYM for the small +T/-T comparison ( Figure fig.simples a) we first pooled all trials from all subjects from both conditions +T and -T; then randomly pick 8 trials, averaged them and assigned them to the "surrogate +T" time series of surrogate subject 1; similarly with other 8 trials to the "surrogate -T" of same subject; these surrogate +T and surrogate -T average series were summed to get the surrogate ASYM for this surrogate subject; we then repeated the steps to generate 17 and 13 surrogate subjects for small and large perturbations respectively in order to match the actual number of subjects; the average across surrogate subjects makes one surrogate grand average ASYM time series, and finally repeated the whole procedure 1000 times to get the null distribution of ASYM time series, i.e. the null distribution of ASYM for each step n. Since we are interested in the transient part of the resynchronization, we restricted the FDR analysis to the range n=1 to n=5, both included (n=0 is excluded because it is the perturbation step and the asynchrony there is a forced error; steps n>5 are too close to the post-perturbation baseline). We then compared the true ASYM value at each step between n=1 and n=5 to the null distribution of ASYM at the same step and compute a p-value as the proportion of null ASYM values above the true value or below its opposite (two-tailed). We applied the same procedure to all comparisons between times series in this work. fig.embedT and  fig.embedST ) . Given a univariate time series e n , there are several possible alternatives to implement time series embedding for phase space reconstruction [Gilmore1998], like the usual time-delay embedding ( e n ; e n -1 ). We chose the following: the mean between consecutive values ( e n + e n -1 )/2 for the first variable and the semi-difference ( e n -e n -1 )/2 for the second variable. Our choice was based on visual clarity and greater trajectory separation (we also tested other embeddings like the above-mentioned time-delay embedding with qualitatively similar results). All reconstructed trajectories in this work were computed directly from the corresponding grand average time series. Normalized trajectories (Figure fig.areas ) . Normalized trajectories were computed by first taking the regular trajectories from 45 and 50 ms perturbation sizes and averaging between them (separately for positive and negative perturbations), then taking each of these and dividing it by its corresponding absolute maximum value (separately for the first and second embedding components). Then we kept the normalized +/-T trajectories as standards and multiplied every trajectory (+/-S, +/-S+/-T, +/-S-/+T) by a factor equal to the ratio between its absolute maximum value and the absolute maximum value of the T trajectories. Areas ( Figure  fig.areas ) . Areas in Figure fig.areas are a rough measure of the region in phase space where the system is left when the perturbation arrives (step n=0), in units of the normalized trajectories. Areas are ellipses with a center and two radii. The center (which roughly corresponds to the first point of the corresponding normalized trajectory) is the average across subjects of the point in trajectory corresponding to n=0, after pooling 45 and 50 ms perturbation sizes (separately for positive and negative perturbations). The radii are the standard deviation in each component.

Color blind-friendly plots
We used the ametrine colormap [GeissbuehlerLasser2013].