Response to perturbations as a built-in feature in a mathematical model for paced finger tapping

TL;DR

This paper introduces a new mathematical model for paced finger tapping that inherently accounts for tempo perturbations, providing better insights into sensorimotor synchronization and error correction mechanisms.

Contribution

The authors propose an alternative variable and a phase space-based model that naturally incorporates perturbations, improving upon traditional asynchrony-based models.

Findings

Model predicts responses to novel perturbations

Trajectories in phase space reveal synchronization dynamics

Applicable to interpersonal and non-periodic synchronization

Abstract

Paced finger tapping is one of the simplest tasks to study sensorimotor synchronization. The subject is instructed to tap in synchrony with a periodic sequence of brief tones, and the time difference (called asynchrony) between each response and the corresponding stimulus is recorded. Despite its simplicity, this task helps to unveil interesting features of the underlying neural system and the error correction mechanism responsible for synchronization. Perturbation experiments are usually performed to probe the subject's response, for example in the form of a "step change", i.e. an unexpected change in tempo. The asynchrony is the usual observable in such experiments and it is chosen as the main variable in many mathematical models that attempt to describe the phenomenon. In this work we show that although asynchrony can be perfectly described in operational terms, it is not well defined as a model variable when tempo perturbations are considered. We introduce an alternative variable and a mathematical model that intrinsically takes into account the perturbation, and make theoretical predictions about the response to novel perturbations based on the geometrical organization of the trajectories in phase space. Our proposal is relevant to understand interpersonal synchronization and the synchronization to non-periodic stimuli.

Tables

Table 1. Table 1: Fitted parameter values. All simulations in this work use this unique set of parameter values unless stated otherwise. Linear coefficients a 𝑎 a , b 𝑏 b , c 𝑐 c , d 𝑑 d are nondimensional; quadratic coefficient δ 𝛿 \delta has units of ms -1 ; cubic coefficients α 𝛼 \alpha , β 𝛽 \beta , γ 𝛾 \gamma have units of ms -2 . The pre-perturbation interstimulus interval is T 0 = 500 subscript 𝑇 0 500 T_{0}=500 ms.

$a$	=	$-0.0485$			$\alpha$	=	$5.67\times {10}^{-5}$
$b$	=	$0.467$			$\beta$	=	$7.71\times {10}^{-5}$
$c$	=	$-0.491$			$\gamma$	=	$9.74\times {10}^{-5}$
$d$	=	$0.987$			$\delta$	=	$4.61\times {10}^{-3}$

Table 2. Table 2: Fitted parameter values for a representative solution of Type-II phase space (Fig. 4 (b)). Units as in Table 1 .

$a$	=	$0.943$			$\alpha$	=	$-1.15\times {10}^{-5}$
$b$	=	$0.176$			$\beta$	=	$-9.29\times {10}^{-5}$
$c$	=	$-0.842$			$\gamma$	=	$7.54\times {10}^{-5}$
$d$	=	$0.0690$			$\delta$	=	$3.50\times {10}^{-3}$

Table 3. Table 3: Fitted parameter values for a representative solution of Type-III phase space (Fig. 4 (c)). Units as in Table 1 .

$a$	=	$0.751$			$\alpha$	=	$4.34\times {10}^{-5}$
$b$	=	$0.0167$			$\beta$	=	$1.72\times {10}^{-6}$
$c$	=	$0.975$			$\gamma$	=	$-5.02\times {10}^{-5}$
$d$	=	$0.371$			$\delta$	=	$-3.80\times {10}^{-3}$

Equations

$$e_n=R_n-S_n$$

$$e_{n+1}=f(e_n,T_n)+\text{mbox}noise$$

$$T_n=S_n-S_{n-1}.$$

$${\Delta }_n=T_n-T_{n-1}$$

$$p_n=e_n+{\Delta }_n$$

$$p_{n+1}=f(e_n,T_n)$$

$$p_{n+1}=f(p_n-{\Delta }_n,T_n)$$

$$e_{n+1}$$

$$x_{n+1}$$

$$F(e_n,x_n,T_n)$$

$$G(e_n,x_n,T_n)$$

$$p_{n+1}$$

$$+F(p_n-{\Delta }_n,x_n,T_n)$$

$$x_{n+1}$$

$$+G(p_n-{\Delta }_n,x_n,T_n)$$

$$p_{n+1}$$

$$+F(p_n-(T_n-s_n),x_n,T_n)$$

$$x_{n+1}$$

$$+G(p_n-(T_n-s_n),x_n,T_n)$$

$$s_{n+1}$$

$$e_{n+1}$$

$$x_{n+1}$$

$$F=-\sqrt{\frac{1}{2\times 26}\sum _{j=1}^2\sum _{n=-5}^{20}{\left(e_n^j-E_n^j\right)}^2+P}$$

Full text

I Introduction

Sensorimotor synchronization (SMS), that is the ability to synchronize our movement to a periodic external stimulus, underlies specifically human behaviors like music and dance Repp (2005) and involves time processing in the millisecond timing range (i.e., hundreds of milliseconds Ivry and Spencer (2004)). SMS is a rare ability among animals and it apparently correlates with the also rare ability of vocal learning Schachner et al. (2009); Patel et al. (2009); Hasegawa et al. (2011), with potential evolutionary implications. The simplest task to study this behavior is paced finger tapping, where a subject is instructed to tap in synchrony with a periodic sequence of brief stimuli (for instance, tones or flashes) like keeping pace with music and while registering the occurrence time of every response. The natural observable and one of the variables most used to quantify the behavior Chen et al. (1997); Repp and Su (2013) is the difference between the occurrence times of every response ($R_{n}$) and the corresponding stimulus ($S_{n}$), called synchrony error or simply asynchrony:

[TABLE]

The asynchrony of a single trial is a relatively noisy time series with a standard deviation of up to a few tens of milliseconds (Fig. 1). Despite that none of the $e_{n}$ is zero, most people can achieve average synchronization with a mean value that is typically negative (called Negative Mean Asynchrony or NMA, hypotetically representing the point of subjective synchrony) Repp (2003). The main goal is to understand how the brain can mantain average synchrony or recover it after a perturbation.

The number of theoretical and experimental works dedicated to understand this behavior and its neural bases is rapidly growing, especially imaging and electrophysiology studies like EEG, MEG, and fMRI Bavassi et al. (2017); Praamstra et al. (2003); Pollok et al. (2008); Bijsterbosch et al. (2011); Nozaradan et al. (2018); Jantzen et al. (2018); Iversen and Balasubramaniam (2016); Merchant et al. (2015). It is very simple to show that there is an error correction mechanism in the brain in charge of keeping average synchrony that operates based on past performance Bavassi et al. (2013); Repp (2005). On the theoretical side Bavassi et al. (2013); van der Steen and Keller (2013), such a mechanism is easily conceptualized as a map or difference equation for the observable $e_{n}$:

[TABLE]

where the asynchrony at the next step $e_{n+1}$ depends on its previous value $e_{n}$ (or several previous values in some models Pressing and Jolley-Rogers (1997)) and probably on some parameter like the sequence period or interstimulus interval $T_{n}$:

[TABLE]

Studies that aim at finding a linear correction function $f$ make use of mean values, standard deviations and auto-correlation functions and thus they analyze synchronization to periodic stimuli sequences without any perturbation Pressing and Jolley-Rogers (1997); periodic sequences are also used by works that study the structure of the noise term Wagenmakers et al. (2004). Alternatively, in order to find the best correction function $f$ that is the deterministic part of the equation one can perform perturbations to an otherwise isochronous sequence and analyze the resynchronization Bavassi et al. (2013). That is the approach we chose.

The main variable for quantifying the behavior, the asynchrony $e_{n}$, is always operationally well defined according to Eq. 1 and Fig. 1. In this work, however, we show that $e_{n}$ is an ill-defined variable in terms of a map or difference equation when the stimulus sequence has perturbations. This issue is relevant not only during a controlled perturbation experiment in a laboratory setting, but also in a more natural, ecological setting like music performance where the stimulus sequence is not strictly periodic (e.g., a choir or orchestra conductor performing a ralentando) or when the stimulus sequence is intrinsically variable (e.g., interpersonal synchronization where the stimulus sequence is the other person’s production). We propose an alternative variable and a mathematical model that reproduces the observed data including the effect of perturbation and make theoretical predictions.

II Results

II.1 Predicted versus observed asynchrony

From a behavioral point of view, one of the approaches to unveil the form of the error correction mechanism in a paced finger tapping task is to find the best correction function $f(e_{n})$ (Eq. 2) that, based on the observed asynchrony $e_{n}$ (Ec. 1) and the interstimulus interval $T_{n}$ (Eq. 3), predicts the asynchrony at the next step. Figure 1 shows a graphical definition of all variables and parameters.

According to Eq. 1, the asynchrony $e_{n}$ takes as a reference point the occurrence time of the stimulus $S_{n}$. The traditional way of perturbing the system is to unexpectedly modify the stimulus period, which is done by shifting in time one or several consecutive stimuli, i.e. modifying $S_{n}$ (and perhaps the following stimuli too; see Fig. 1). An example is the “step change” perturbation where the stimulus period changes unexpectedly by a fixed amount $\Delta$ from a given stimulus $S_{n}$ on (in musical terms it is a sudden change in tempo). A critical issue that has been overlooked in the literature, both experimental and theoretical, is that the change of the occurrence time of the perturbed stimulus $S_{n}$ is arbitrary and thus is not a well defined time reference López and Laje (2019). The consequence of this is that the variable $e_{n}$ becomes ill-defined because its value changes at the moment of perturbation but not because of its own dynamics. In other words, if an unexpected perturbation occurs at step $n$, the actual asynchrony will be different from the value predicted by the subject (or the model) because the corresponding stimulus was shifted by an arbitrary amount equal to the size of the perturbation.

In this work we propose the following way to resolve this issue: we distinguish between the predicted asynchrony value $p_{n}$ and the actually observed asynchrony value $e_{n}$ (Fig. 1). If a change in period occurs at step $n$:

[TABLE]

then the predicted asynchrony $p_{n}$ and the actually observed asynchrony $e_{n}$ are related by the following expression:

[TABLE]

It is important to note that Eq. 5 is not a theoretical assumption but the actual relationship between the variables. Note also that if no perturbation occurs then $\Delta_{n}=0$ and thus $p_{n}=e_{n}$.

Taking this distinction into account, a model like Eq. 2 must be written as (discarding the noise term for simplicity):

[TABLE]

and by using Eq. 5 we can get a closed model for the variable $p_{n}$:

[TABLE]

We can always go back to the observed asynchrony $e_{n}$ by means of Eq. 5.

II.2 Model implementation

In this section we proposed an improved model based on previous work Bavassi et al. (2013) and apply the distinction proposed above. In Bavassi et al. (2013) we showed that nonlinear effects are important even for small perturbations of 10% of the period. There are two main effects when step change perturbations are performed. If the tempo is made faster (negative perturbation $\Delta_{n}<0$; Fig. 2(a)) the resynchronization is monotonic until a new baseline is reached; however, if the tempo is made slower (positive perturbation $\Delta_{n}>0$; Fig. 2(b)) the asynchrony overshoots before reaching the new baseline. The existence of the overshoot makes necessary to introduce a second variable $x_{n}$—otherwise the deterministic nature of the model would be violated Bavassi et al. (2013); Schoner (2002). In addition, the asymmetry of the response in front of symmetric perturbations requires the introduction of nonlinear terms of even order (e.g. quadratic). In order to correctly reproduce these two observations we proposed a two-variable, nonlinear model Bavassi et al. (2013):

[TABLE]

where $F$ and $G$ are nonlinear functions of their arguments.

Before introducing the distinction between predicted and observed asynchrony, we propose improved functions $F$ and $G$ based on new experimental evidence. We take the experimental data of step change perturbations from Bavassi et al. (2013) and reanalyze them as follows. First we perform an embedding of the time series shown in Figs. 2(a) and 2(b) to reconstruct the qualitative geometrical arrangement of the underlying phase space Gilmore (1998). Figure 2(c) shows the result, where we added the detection of return points for every trajectory (see methods in Appendix A.1). The geometrical organization of the return points shows a certain degree of asymmetry and saturation (especially for the vertical axis), which tells us that quadratic and cubic terms are needed. There are six possible quadratic terms and eight possible cubic terms, and we want to use the smallest number of terms that correctly represent the behavior. However, we cannot use normal form theory to choose among the nonlinear terms because within the regime analyzed in this work (synchronization to a periodic sequence or resynchronization following a tempo step change of fixed size) the behavior does not show any bifurcations but a robust convergence to a single fixed point representing average synchrony Bavassi et al. (2013); Thaut et al. (1998); Repp (2001); Repp and Keller (2004). After testing many combinations of nonlinear terms with qualitatively similar results, we choose the following:

[TABLE]

In Section II.3 we display a summary of the obtained phase spaces with our selection of nonlinear terms. We emphasize that our choice of nonlinear terms is not unique—it is not possible to solve for unique $F$ and $G$ from the shape of the return points in the embedding.

Now we incorporate the distinction between predicted and observed asynchrony. According to the previous subsection we must substitute $e_{n+1}\rightarrow p_{n+1}$ and $e_{n}\rightarrow p_{n}-\Delta_{n}$:

[TABLE]

and to get a closed system we define the auxiliary variable $s_{n+1}=T_{n}$. This leads us to

[TABLE]

where $F$ and $G$ are defined according to Eq. 9. It is worth noting that $T_{n}$ is a parameter whose value is up to the experimenter.

II.3 Model fitting, model simulations, and fitting analysis

As we did before Bavassi et al. (2013), we use a genetic algorithm to fit the model (Eq. 11) to the experimental time series of $\Delta=\pm 50$ ms only (Fig. 3(a); see methods in Appendix A.2). The obtained parameter values are shown in Table 1 and the corresponding phase space is shown in Fig. 3(b).

The fitting goodness is noteworthy, especially when compared to previous attempts in the literature with a similar or even greater number of parameters (for instance our own previous work Bavassi et al. (2013) and many others’ work, e.g. Thaut et al. (1998); Repp (2001); Schulze et al. (2005); Large et al. (2002); Repp and Keller (2004); Repp et al. (2012)). It must be taken into account that, contrary to ours, modeling and fitting efforts in the finger tapping literature usually perform separate fits for different perturbation sizes, effectively multiplying the number of parameters by a factor of 2 at least.

In order to analyze the robustness of the fitting results we proceeded as follows. Every run of a genetic algorithm provides a solution that best fits the data according to the fitness function. As we describe in the Appendix A.2 (also in our previous work Bavassi et al. (2013)) we decided to run the algorithm 200 times to choose the absolute best among those 200 converged solutions. This, in addition, allows us to perform a statistical and dynamical analysis of all solutions. The distributions of obtained values for every parameter can be seen in Appendix A.3. All distributions are mostly unimodal, which speaks in favor of the robustness of the fitting procedure (with a slight bimodality in some of them). On the other hand, there is some interdependency between some of the linear parameters (i.e., a correlation between their converged values, for instance between $a$ and $d$, but not between nonlinear parameters, but a small correlation among coefficient $\alpha$ and coefficients $a$ and $d$ (see Appendix A.3). This means that the model might be overparameterized and might fit the experimental time series with a smaller number of linear terms, but our choice of nonlinear terms is appropriate.

We plotted the phase space of every obtained solution and made an exhaustive visual search in order to qualitatively classify them. We found three types of phase spaces (Fig. 4). The phase space we chose for this work (Figs. 3(b) and 4(a)) is the fittest of all 200 solutions and in addition belongs to the most frequent type of phase space (type I); the parameter values corresponding to phase space types II and III are in Appendix A.4).

II.4 Benefits of the proposed approach

The most important feature of our model is that it incorporates the perturbations to the stimulus period in the modeling approach, leading to autonomous dynamics once the stimulus period sequence $T_{n}$ is chosen as input. We achieved this by developing a closed equation for variable $p_{n}$ (Eq. 11) based on a model-free relationship between $p_{n}$ and the observable $e_{n}$ (Eq. 5). An advantage of having a closed model with built-in perturbations is the ability to perform a bifurcation analysis on it, which is the subject of future work.

In the absence of bifurcations, as is the case in this work, our approach still offers advantages in the form of autonomous dynamics without any need to modify the value of the variable “by hand” whenever the stimulus period changes. This is illustrated in Fig. 5, where we show three common experimental manipulations of the stimulus period $T_{n}$ (panel (a), step change; panel (b) sinusoidal variation; panel (c) random variation). Once the stimulus period sequence is set (top row), our model evolves in time without any intervention from the experimenter (second row). The traditional way of doing this (third row) consists in considering the model without distinguishing between $p_{n}$ and $e_{n}$, that is Eq. 8, but this of course needs adjusting the value of the variable $e_{n}\rightarrow e_{n}-\Delta T$ by hand (shown in the figure as orange circles) every time there is a change in the parameter $T_{n}\rightarrow T_{n}+\Delta T$ to produce the correct time evolution (in the sinusoidal and random variations the experimenter needs to adjust $e_{n}$ at every step). Finally, a naive version is shown in the bottom row where Eq. 8 is used without adjusting the value of $e_{n}$ by hand.

Note that we decided to plot $e_{n}$ instead of $p_{n}$ throughout this work mainly for historical reasons so it is easier to compare to previous models’ results and experimental results. Keep in mind, however, that all numerical simulations of our proposed model in this work were made by solving the closed equation for $p_{n}$, Eq. 11, and then translated $p_{n}\rightarrow e_{n}$ by means of Eq. 5.

II.5 Theoretical predictions

II.5.1 Geometrical organization of phase space

The geometrical arrangement of the trajectories in the experimental phase space is remarkable (Fig. 2(c)). The trajectories corresponding to the largest positive perturbations share a region of phase space during the overshoot with some of the trajectories of the negative perturbations (that do not have overshoot). This suggests that the error correction mechanism in the brain might not distinguish between both types of trajectories while they are roughly in the same region of phase space. The approximate location of such region is clear and is labeled “B” in the bottom panel of Fig. 6(a), but corresponds to very different time instants in the time series (top panel).

The comparison between the reconstructed phase space (Fig. 2(c)) and the model phase space (Fig. 3(b)) is to be understood as a qualitative similarity in the geometrical arrangement of the trajectories: both have trajectories with asymmetric overshoot and both have trajectories that share the same region in phase space despite coming from opposite perturbations. Our model allows us to make the following prediction: a perturbation to the variable in the region labeled “B” in phase space (see Fig. 6(a)) should show the same post-perturbation time evolution no matter what original trajectory belongs to.

For the very same reasons exposed in Section II.1, in the history of paced finger tapping experiments it has been intrinsically difficult to perturb the value of the stimulus period $T_{n}$ without perturbing the variable $e_{n}$, and conversely to perturb the variable without perturbing the stimulus period. In a recent work, however, we showed the experimental feasibility of such manipulations López and Laje (2019). We then propose to use these novel perturbations to study the response of the model in front of two consecutive perturbations:

first a traditional step change perturbation (in two conditions, positive and negative); 2. 2.

second, and while the resynchronization from the first perturbation takes place, a perturbation to the variable without changing the stimulus period.

Our prediction is that the time evolution following the second perturbation will be the same for both conditions (positive or negative first perturbation), provided the second perturbation is performed when the system is approximately in the same region of phase space. The rationale behind this proposal is that, if we on the contrary performed two consecutive traditional step change perturbations, we would not be able to resolve the following confounder: in case an overshoot appears in both conditions after the second perturbation it may be due either because our hypothesis is valid or because it is the known response to the (second) step change perturbation.

Numerical simulations supporting our prediction are in Fig. 6, where we show the results of (a) a single step change (negative control), and the definition of perturbation points; (b) the proposed manipulation; (c) two consecutive perturbations as in (b) but performed in different points of phase space (positive control). In the consecutive perturbations (comparison (b)-(c)) the time evolution after the second perturbation is either very similar between conditions (b) or different (c).

II.5.2 Perturbations to the variable only

Our second prediction is that the response to large enough symmetric perturbations to the variable might be asymmetric. This can be seen in the phase space of Fig. 7, after noting that large negative perturbations (i.e., jump to the left) have a larger overshoot than large positive perturbations (i.e., jump to the right). This would prove that a step change perturbation is not needed to display nonlinear behavior—perturbations to the variable without changing the stimulus period might elicit it too. On the other hand, if the perturbation size is small enough then the response to symmetric perturbations might be symmetric (see Fig. 7).

II.5.3 Large perturbations

We note the saddle node appearing in the upper right quadrant of the phase space, Fig. 3(b). Its stable manifold is a separatrix between bound trajectories (i.e., converging to the origin) and unbound trajectories (i.e., diverging). The existence of the saddle node was not explicitly considered in the modeling effort; a post-hoc interesting interpretation is that it separates successful from unsuccessful resynchronization responses (bound and unbound trajectories, respectively). The unbound trajectories develop only when the perturbation is large enough, meaning the subject is not able to successfully resynchronize after large perturbations. Unbound trajectories display asynchronies with eventually diverging values; within the scope of our model, a very large asynchrony means that the response is to be asociated to the following (or previous) stimulus. Our model does not take into account this reassociation process and thus the validity of our description ends when the asynchrony crosses a threshold value, for instance half period.

Lastly, the predictions described above are valid if the perturbations are performed on time series displaying asymmetric overshoot, otherwise predictions and conclusions would be flawed.

III Conclusions

We showed that the asynchrony $e_{n}$, the most important variable in the literature of paced finger-tapping experiments and theoretical models Chen et al. (1997); Repp and Su (2013), is actually an ill-defined variable for a map or difference equation model when perturbations to the stimulus sequence are present. We proposed a distinction between predicted and actually observed asynchrony and developed the first mathematical model to solve the inherent ill-definition. This is also the first mathematical model in sensorimotor synchronization that takes into account the response to a temporal perturbation as a built-in feature. This is an important issue when considering interpersonal or group synchronization and leader-follower relationships (like in choirs and orchestras), as any naturally occurring variability in the timed actions of any participant or tempo change by the leader will act as a perturbation to the rest. Our own previous attempt Bavassi et al. (2013), though successful at unifying, fell short of completely including the perturbation in the model dynamics.

Our model is able to fit the step change perturbation data remarkably well (Fig. 3(a)), reproducing the response time series very accurately at all steps (including the perturbation step) with no modification of the variable “by hand” and with basically the same number of parameters than comparable models. It is only surprising that no published works so far in the paced finger tapping literature, up to our knowledge, deal with the issue of ill-definition of the main model variable when the sequence period is perturbed. On the other hand, models based on forced or coupled nonlinear oscillators Haken et al. (1985); Large et al. (2002); Loehr et al. (2011); Egger et al. (2019), traditionally classified as belonging to a “dynamical systems” approach, are naturally well defined even in the presence of perturbations to the period.

Temporally displaced auditory feedback (either delayed or advanced) is also a usual way of probing the system Aschersleben and Prinz (1997); Wing (1977); Pfordresher and Dalla Bella (2011); Mates and Aschersleben (2000) and it does not suffer from the issue of ill-definition of the variable we addressed in this work. It remains to be shown, however, how it relates to changing the asynchrony in the models as it produces a modification in the asynchrony value but it also introduces a dissociation between auditory feedback and proprioceptive and tactile feedback Pfordresher and Dalla Bella (2011).

We also showed that nonlinear behavior (asymmetry of responses) might be observed when the variable only is perturbed, i.e. even in the absence of a perturbation to the stimulus sequence, if the perturbation is large enough. Experimental perturbations like the ones proposed in López and Laje (2019) but with larger magnitudes are needed. We acknowledge that similar results can be obtained by using a different set of nonlinear terms, and this calls for more experimental data showing any kind of bifurcation in the behavior so as to choose the correct set of parameters via normal form theory.

Our model assumes, as many others, that the origin of the asynchrony is not important for keeping average synchrony or for achieving resynchronization. Qualitative features of resynchronization after a perturbation are thus similar independently of whether the asynchrony was produced by a perturbation to the parameter or to the variable or both. This is a common theoretical assumption in the literature only recently supported by experimental results López and Laje (2019).

Our theoretical results show that past (observed) and future (predicted) asynchronies play different roles in the model, and the remarkable fitting to the experimental data thus offers indirect evidence for a separate cerebral account of predicted versus actually observed asynchrony. Further experimental work is needed to decide whether this holds true.

Acknowledgements.

This work was supported by Universidad Nacional de Quilmes (Argentina), CONICET (Argentina), and The Pew Charitable Trusts (USA). Author contributions: CRG and RL analized data, wrote code and performed numerical simulations; RL conceptualized the work and developed the model; CRG, MLB and RL wrote the manuscript.

Appendix A Methods and Parameter distributions

A.1 Estimation of return points

The return points of a map $e_{n+1}=f(e_{n})$ are the values $e_{n}$ such that $e_{n+1}=e_{n}$ (the analogous concept in a continuous-time flow is the nullcline, that is the points in phase space where the rate of change associated to a given variable is zero). In the time series $e_{n}$ the return points appear as local maxima or minima.

In a 2D map:

[TABLE]

the return points for the variable $e_{n}$ are the solutions of the implicit function $e_{n}=f(e_{n},x_{n})$ and for the variable $x_{n}$ are the solutions of the implicit function $x_{n}=g(e_{n},x_{n})$. Here it is clear that the return points give us information about the shape of $f$ and $g$.

In our case we want to find the return points in the chosen bidimensional embedding $(e_{n};z_{n})$, where $z_{n}=e_{n}-e_{n-1}$. That is, on the one hand we want the points $(e_{n};z_{n})$ such that $e_{n}$ is a local maximum or minimum (return points for the variable $e_{n}$), and on the other hand the points $(e_{n};z_{n})$ such that $z_{n}$ is a local maximum or minimum (return points for the variable $z_{n}$).

We exemplify the procedure with the calculation for the variable $e_{n}$ (see Figure 8):

Define a 5-point time window from $n=0$ through $n=4$; 2. 2.

Fit a 4th-order polynomial to the $e_{n}$ time series in such window; 3. 3.

Find the local maximum or minimum $e_{\mbox{ret}}$ of the fitted function and the corresponding value $n_{\mbox{ret}}$; 4. 4.

Interpolate the $z_{n}$ time series with a polynomial of the same order and compute the value $z(e_{\mbox{ret}})$; 5. 5.

The return points are the set of values $(e_{\mbox{ret}};z(e_{\mbox{ret}}))$.

Same procedure applies for the variable $z_{n}$ after switching $e_{n}\leftrightarrow z_{n}$.

A.2 Genetic algorithm

We fitted the model to the data by using a genetic algorithm in C with both custom-written code and the GAUL libraries (http://gaul.sourceforge.net).

The eight model parameters were arranged into a single chromosome with eight genes, and were initialized randomly from a uniform distribution in the ranges $-1.0<a,b,c,d<1.0$, $-0.01<\delta<0.01$, and $-0.0001<\alpha,\beta,\gamma<0.0001$. The number of generations was 200, the population size 2000, the crossover rate $0.9$, and the mutation rate $0.1$.

The fitness function was defined as minus the square root of the average squared deviation between model series and experimental series:

[TABLE]

where $e_{n}^{j}$ is the experimental time series and $E_{n}^{j}$ is the model time series; $n=-5,\dots,20$ is the step number as in Fig. 3(a), and $j=1,2$ represents the two conditions $\pm 50$ ms used to fit the model. The fitness function $F$ decreases as the differences $e_{n}^{j}-E_{n}^{j}$ get larger in absolute value. Note that fitting and simulations in this work encompass all steps $n=-5,\dots,20$ in a sequence including the perturbation step $n=0$ as shown in Fig. 3(a), making no distinction whatsoever among pre-perturbation, perturbation, and post-perturbation.

In order to prevent survival of unrealistic solutions (for instance damped oscillations or alternating series), penalties were included as a positive term $P$ inside the square root that depends on the linear coefficients only and takes a large value in any of the following cases:

the eigenvalues are complex (in order to avoid oscillatory approach to the equilibrium); 2. 2.

the eigenvalues are real but any of them is either greater than 1 or negative (in order to avoid solutions with unstable manifolds, and convergent solutions that alternate sides);

otherwise $P=0$.

In order to prevent the selection of a surviving local optimum and to perform a statistical and dynamical study of the obtained solutions, the whole procedure described so far was repeated 200 times; the chosen solution was the one with the highest fitness of all.

To improve fitting, a post-perturbation constant baseline was added to the model variable $p_{n}$ with a fixed value equal to the experimental post-perturbation baseline of the corresponding perturbation size.

A.3 Fitted parameter distributions

See Figures 9, 10, 11, and 12.

A.4 Parameter values of Type-II and Type-III phase spaces

A.5 Data and code

See Supplemental Material Sup at [URL will be inserted by publisher] or at the Sensorimotor Dynamics Lab’s webpage: www.ldsm.web.unq.edu.ar/perturbations2019 for C and MATLAB code to reproduce all figures and data in this work.

We use the morgenstemning colormap Geissbuehler and Lasser (2013) for color blind-friendly and grayscale-friendly plots.