## Abstract

Using a rhythmic task where human subjects bounced a ball with a handheld racket, fine-grained analyses of stability and variability extricated contributions from open-loop control, noise strength, and active error compensation. Based on stability analyses of a stochastic-deterministic model of the task—a surface contacting the ball by periodic movements—open-loop or dynamic stability was assessed by the acceleration of the racket at contact. Autocovariance analyses of model and data were further used to gauge the contributions of open-loop stability and noise strength. Variability and regression analyses estimated active error compensation. Empirical results demonstrated that experienced actors exploited open-loop stability more than novices, had lower noise strength, and applied more active error compensations. By manipulating the model parameter coefficient of restitution, task stability was varied and showed that actors graded these three components as a function of task stability. It is concluded that actors tune into task stability when stability is high but use more active compensation when stability is reduced. Implications for the neural underpinnings for passive stability and active control are discussed. Further, results showed that stability and variability are not simply the inverse of each other but contain more quantitative information when combined with model analyses.

## INTRODUCTION

It is widely recognized that stability is a core requirement for many motor tasks. Despite external perturbations and variability generated by the sensorimotor system itself, humans and animals can maintain static postures and perform movements consistently and reliably. The neurophysiological system has developed a variety of mechanisms at different levels to achieve stable postures and movements, ranging from mechanical stiffness of the muscular and tendonous tissues (Loeb 1995; Wagner and Blickhan 1999), to stiffness generated by muscular reflexes with short and long loops (Nichols and Houk 1976; Popescu et al. 2003), to voluntary intervention on the basis of perceptual information (Hasan 2005; Loram and Lakie 2002; Morasso and Sanguineti 2002). Besides these mechanisms, there is another contribution to stability originating from the interaction between the actor and the environment. Since humans interact with objects in the environment that pose physical constraints but also afford inherent stability, further understanding of the challenges for the neurophysiological control system can be derived from analysis of the dynamics of this interaction in the context of a task (Saltzman and Kelso 1987; Warren 2006).

Some tasks are inherently stable, others inherently unstable. Upright posture when viewed as an inverted pendulum is unstable but it is stabilized by active control processes in the central and peripheral nervous system, as evidenced by the continuous sway during quiet standing. Other behaviors have inherent stability and require no additional control. For example, your arm hanging at the side of your body does not require active muscle activations and perturbations die out without error correction. Although this trivial example shows static stability, more complex behaviors can have dynamic stability. One illustrative example is the passive dynamic walker. McGeer (1990) first demonstrated that a mechanical biped could walk down a gentle slope without any control. Small fluctuations naturally arising in the stepping motion were dissipated, solely by relying on the stability properties of this passive mechanical system. Based on the understanding of the passive dynamic stability, actuated robots have been designed that exploit the self-stabilizing properties and thereby can function with minimal control (Collins et al. 2005). The implication for the study of biological systems is that the dynamic interactions of the actor with the physical environment can similarly afford passively stable coordination strategies. We and others have hypothesized that the CNS seeks solutions to a task that exploit such passive stability to reduce the burden on control (Warren 2006). Hence, the research strategy in the present work is first to examine the stability that is afforded by the task without control and, on this basis, to identify additional features of active control that the nervous system may apply. Specifically, the present study pursues to identify active and passive contributions to the dynamic stability of steady-state rhythmic bouncing of a ball with a racket.

This complex interactive task and its dynamics have already been investigated in a series of studies on human sensorimotor control but also in robotic and theoretical studies (Bühler et al. 1994; de Rugy et al. 2003; Dijkstra et al. 2004; Guckenheimer and Holmes 1983; Rizzi et al. 1992; Ronsse et al. 2006; Schaal et al. 1996; Sternad 2006; Tufillaro et al. 1992; Wei et al. 2007). The task requires the actor (or actuator) to use a handheld racket to rhythmically bounce a ball vertically into the air to a consistent height. Using a dynamical model of the task, mathematical stability analyses revealed that the system has a stable solution if the racket is decelerating in its upward movement before contacting the ball (Fig. 1*A*; for a review of the model see appendix A). This criterion can be intuited with the help of Fig. 1*B*, which shows a time series of racket and ball: if the ball arrives early after a low-amplitude bounce, as for example the impact of the “noisy” (dashed) trajectory at 2.3 s, then the ball is hit with a relatively higher velocity, leading to a higher ball amplitude in the subsequent cycle. Thus the condition for stability is that the ball is hit with a higher velocity when the impact is relatively early and with a lower velocity when the impact is late. This condition is equivalent to stating that the impact should happen during the decelerating portion of the racket trajectory. For an actor moving a racket in this way, this means that small fluctuations or perturbations from this stable racket trajectory do not have to be explicitly corrected; stability can be achieved “passively” (i.e., without error correction or closed-loop control). A series of previous empirical studies of our group demonstrated that actors indeed hit the ball with negative acceleration, consistent with the interpretation that humans exploit dynamical—passive—stability (e.g., Schaal et al. 1996; Sternad 2006; Wei et al. 2007).

Note that the comparison of the mechanical model with the ball–racket system of the actor necessitates some terminological clarification: To obtain passive stability as defined by model analysis the actor still has to (actively) generate a rhythmic racket movement; however, it is sufficient that the actor controls the racket movement in an open-loop manner. In contrast, when the actor applies explicit active error-compensating mechanisms this implies closed-loop control. Hence, in the discussion of human performance we equate passive stability with open-loop control, whereas active control is equated with closed-loop control.

These two control strategies, however, are not simple alternatives. The dynamical stability of the ball-bouncing system can vary in degree, as shown in more detailed analyses in appendix A. This degree of stability, as quantified by the absolute eigenvalue of the Jacobian, is dependent on its parameters, in particular the coefficient of restitution α. The coefficient of restitution expresses the amount of energy lost at the contact.1 Although it is straightforward to expect that control requirements change when more energy has to be imparted on the ball to achieve the same ball amplitude, the mathematical result is not as easily intuited and the reader is referred to appendix A. Given these new model predictions the present study will manipulate the coefficient of restitution and assess whether different degrees of stability lead to different performances.

The finding that actors can exploit passive stability does not necessarily imply that they always use only open-loop control. For example, a learning study revealed that novice actors initially performed the bouncing task with positive-impact acceleration and gradually tuned to the passive regime with negative-impact acceleration only after 30 min of practice (Dijkstra et al. 2004). Importantly, novices could maintain a rhythmic bouncing pattern similar to that of more experienced actors, although their performance variability was higher. It could be concluded that they used an error-correcting closed-loop strategy. A similar difference in strategy was observed across selected individuals with more or less experience in the task (Sternad et al. 2000). These findings suggest that using negative impact acceleration is not an intuitive solution for the actor but rather one that it is learned with practice.

Further, perturbation studies suggested that actors used a mixed strategy involving both open-loop control exploiting passive stability and active control mechanisms based on perceived error information. When experimentally controlled perturbations were applied to the ball trajectory, actors adjusted their racket movements according to target errors, even though racket impact acceleration values returned to the negative range immediately after the perturbation (de Rugy et al. 2003; Wei et al. 2007). As a result, the errors introduced by the perturbation were corrected faster than the relaxation processes of the model would have predicted, suggesting that actors used a blend of passive stability and active control.

Given this evidence in the case of large perturbations it remains an open question whether there are also active control components intertwined with the observed passive stability when the task is performed at steady state. At steady state there are always small fluctuations or perturbations induced by internal and external “noise.” Do actors completely rely on the slower self-correcting properties of the task? Whereas the two previous experiments examined the response following explicit relatively large perturbations, identification of error-compensating strategies during unperturbed steady-state performance requires a detailed time-series analysis of the sequence of bounces. However, as pointed out earlier, error-correcting relations naturally exist and are characteristic of dynamic stability. Therefore a simple autocorrelation analysis of bounce-to-bounce fluctuations will not be able to distinguish between the inherent stabilizing relations and error-compensating control by the actor.

To tease apart the self-stabilizing from the active corrections we extended the current deterministic ball-bouncing model with a stochastic component representing the inherent fluctuations during dynamically stable performance. Covariance analysis of the model across bounces then quantifies the self-correcting properties of the open-loop model. Any deviations from this prediction in an actor's performance can then be attributed to additional active corrections. Figure 1*B* illustrates the model behavior: the ball trajectory with a solid line is a simulation of stable behavior from the deterministic model with no noise; the dashed trajectory shows a simulation where each racket cycle is slightly perturbed by adding random noise to the amplitude and period of the racket. Assuming that these random perturbations are small, they propagate across bounces with a specific relaxation behavior. The ball-bouncing map at dynamic stability will produce a structure in the fluctuations that can be quantified by the covariance of the state variables over successive contacts. Based on the analyses of the autocovariance functions of this stochastic-deterministic model, additional more fine-grained quantitative predictions were derived (appendix B). In short, the auto- and cross-covariances between the state variables of the model, velocity after impact, and time between impacts, are predicted to approach zero over lags. The time of return to zero is inversely related to the degree of stability of the system. As stated earlier, the degree of stability depends on the coefficient of restitution α. Testing these predictions in human performance will provide a fine-grained assessment of how sensitive the physiological system is to stability properties of the task.

To this end the present empirical study examined steady-state bouncing under seven conditions of α. In addition, because different skill levels of actors were previously shown to lead to different performance strategies, we compared novice and expert participants. Specifically, we examined whether skilled actors showed higher contributions of open-loop stability compared with error-based corrections. The predictions can be summarized as follows:

#### Prediction 1a.

Consistent with previous analyses, performance with dynamic stability is indicated by negative accelerations of the racket at contact. This passive stability is expected for all α conditions, yet with a lower degree of stability for higher α.

#### Prediction 1b.

Consistent with previous empirical findings, it is expected that experienced actors rely more on dynamic stability. This is evidenced in more negative accelerations at contact in the performance of experts, in contrast to more positive values in nonexperienced actors.

#### Prediction 2.

Stability, as quantified by the autocovariance function, decreases with increasing coefficients of restitution α. This prediction is expressed in an autocorrelation structure in the data with significant auto- and cross-correlations over longer lags for higher α conditions.

Alternatively, it may be expected that with decreasing passive stability at higher α values, the CNS utilizes more active control to stabilize performance.

#### Prediction 3a.

Given lower stability for higher α conditions, such active corrections are expected to be more prominent at higher α conditions.

#### Prediction 3b.

Experienced actors require less active corrections to stabilize performance and decrease performance variability.

The results with respect to these predictions will give insight into how sensitive human actors at different degrees of expertise are to dynamical properties of the task. The data will provide a rich basis to assess how different neurophysiological processes, both corrective and anticipatory, lead to stable task performance. The results will be discussed in light of new perspectives on the control principles that the CNS uses.

## METHODS

### Participants

Eight volunteers, ranging in age from 21 to 47 yr, participated in this experiment. Four participants had previously taken part in several ball-bouncing experiments and were considered experienced in this task. Four other volunteers had never performed the task prior to this experiment. All participants were right-handed and used their preferred right hand to bounce the ball with the racket. Before the experiment, all participants were informed about the procedure and signed the consent form approved by the Regulatory Committee of the Pennsylvania State University.

### Experimental apparatus

Participants manipulated a real table tennis racket to bounce a virtual ball that was projected onto a backprojection screen (2.5 × 1.8 m) in front of them (Fig. 2). Using custom-made software a PC controlled the data collection and generated the images that were projected onto a backprojection screen (60-Hz refresh rate). Accelerations were measured directly using an accelerometer mounted on top of the racket. The mechanical brake on the rod attached to the racket was controlled by a solenoid. A light rigid rod with three hinge joints was attached to the racket surface and ran through a wheel whose rotation was registered by an optical encoder. The racket could move and tilt with minimal friction in three dimensions but only the vertical displacement was measured. Racket movements in other than vertical directions were minimal because only vertical movements were displayed on the screen and interacted with the virtual ball. The delay between real and virtual racket movements was 22 ms on average with SD of 0.5 ms (for more details on the experimental setup see Wei et al. 2007).

The virtual racket was displayed at the same height as the real racket and its displacement was the same as that of the real racket (the gain between real and virtual movements was 1). The ball displayed on the screen was a virtual ball and its movement was governed by ballistic flight and an instantaneous impact event when the virtual racket impacted the virtual ball. Just before the virtual ball hit the virtual racket a trigger signal was sent out to the mechanical brake that was attached to the rod. The trigger signal was sent out 15 ms before the ball–racket contact to overcome the mechanical and electronic delay of the brake. The brake applied a brief force pulse to the rod to create the feeling of a real ball hitting the racket. Thus the impact caused by the brake and the impact observed on the screen coincided. The duration of the force pulse (30 ms) was consistent with the impact duration observed in a real ball–racket experiment (Katsumata et al. 2003).

### Procedure and experimental conditions

Each trial began with a ball appearing at the left side of the screen and rolling on a horizontal line extending to the center of the screen. On reaching the center, the ball dropped from the horizontal line (0.8 m high) to the racket. The task for the participant was to bounce the ball periodically to the target line (the line on which the ball started). The experiment was conducted in 14 blocks with seven different α values (0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9). Each α condition was presented in two blocks with three trials in each block. For four participants (two from each group), the blocks were presented first in ascending and then descending order of α values, i.e., from 0.3 to 0.9 then from 0.9 to 0.3. For the other four participants (two from each group), the blocks were presented first in descending order and then in ascending order. Each trial lasted 40 s.

### Data reduction and analysis

The first 5 s of data were eliminated from the analysis because subjects stabilized their performance in the initial part of the trial. The raw data of the racket displacement and acceleration were resampled at a fixed frequency of 400 Hz and filtered by a fourth-order Savitzky–Golay filter with a window size of 0.01 s on both sides (Gander and Hrebicek 2004). Compared with conventional filters like Butterworth filters, the Savitzky–Golay filter is superior for smoothing data that have abrupt changes. These abrupt changes occurred when the racket exhibited a sudden drop in acceleration at impact, caused by the brake. The filter order and window size were chosen empirically to remove the measurement noise while not excessively smoothing the signals. The ball-displacement data were generated by the computer and contained no measurement noise. Therefore no filtering was performed.

As a verification of our setup, the racket displacement was double-differentiated using a Savitzky–Golay filter and compared with the acceleration data from the accelerometer (see Wei et al. 2007). On average the accelerometer signal rendered estimates for the impact acceleration that were 0.52 m/s^{2} more positive than the position signal (compare with mean values in Fig. 4). To be conservative and given that the variability of the impact accelerations from the accelerometer was smaller we opted to report the estimates from the accelerometer.

### Dependent measures

Performance was evaluated by the ball height error *HE*, defined as the signed difference between the maximum ball height and the target height (Fig. 1*B*). The amplitude of the racket movement *A* was calculated as half the distance from the minimum to the maximum of the racket trajectory during one cycle. The racket period *T* was calculated from the intervals between the times of peak velocities of successive bounces. The acceleration of the racket at impact *AC* was determined from the raw accelerometer signal one sample before the start of the impact. The mean and SDs for *HE*, *A*, *T*, and *AC* were calculated over occurrences (bounces) within trials, then averaged across trials per condition, and finally averaged over the four participants of each group.

### Correlation function

The vector autocorrelation function *ACF* was calculated from the vector covariance function of the state variables ball velocity at release *v*, and the time between impacts *t* (Fig. 1*B*). The *ACF* captures the correlation of state variables across lags, where lags correspond to bounces in the data and model. The *ACF*, normalized by the lag-0 covariance, encompasses both the autocorrelations *ACF*(v-v) and *ACF*(t-t) and the cross-correlations *ACF*(t-v) and *ACF*(v-t) for ball release velocity and impact period. For a given lag number the vector autocorrelation function equals the Pearson product-moment correlation of state variable (e.g., release velocity) and a lagged copy of a state variable (e.g., time between impacts). Before calculating the *ACF* for each trial, the series of *v* and *t* were detrended because some participants showed slow drifts in their impact position or ball amplitudes. To calculate the mean *ACF*, the *ACFs* of all trials were averaged for a given α and for all participants within a group. To quantify the variability of the *ACF*, the SDs over six trials per α condition were calculated for each participant.

To calculate the *ACF* from the model, the four deterministic model parameters had to be specified. The two noise strengths did not matter because they were divided out by the normalization inherent in the calculation of the *ACF*. Two parameters were specified by the experimental design: *g* was fixed at −9.81 m/s^{2} and α varied with condition. The other two parameters were set from the average values from each group: the average bouncing frequency ω_{r} and amplitude *a*_{r} (*Eq*. A*4*). We chose to fit the parameters at the level of the group, not individual actors or individual trials, because impact acceleration was positive for several trials in the nonexperienced group. The impact phase, which was used in the calculation of equivalent amplitude (*Eq*. A*4*), was calculated from the time difference between time of maximal velocity and time of impact divided by the period *T*. Several choices are possible in calculating an average equivalent frequency and phase from experimental data: from an individual bounce, from a trial average, from a condition average, or from a group average. Since the comparisons are between α conditions and groups, we chose to average over bounces in a trial and over trials of the same condition but not over participants in a group nor over conditions nor groups. From these four parameters, the Jacobian and the *ACF* could be calculated (*Eqs*. B*3*, B*4*, and B*5*).

A problem arises when the impact phase is negative because the *ACF* is not defined in this case. Therefore we eliminated trials for which impact phase was negative. This led to the elimination of 46 trials (14% of all trials). These trials occurred mainly in three nonexperienced participants in the lower α conditions. Note that the trials were eliminated for calculating only the theoretical predictions. The corresponding data were not eliminated when calculating performance measures because we did not regard it justified to eliminate data that did not fit the theory.

### Noise level

The noise strengths *q*_{v} and *q*_{t} are additional parameters that were introduced by the stochastic extension of the deterministic model (Appendix B). The noise strengths were estimated from the experimental data by the unexplained variance in *Eq*. B*1*. In detail, we first obtained an estimate of the Jacobian, as described in the appendix and then calculated the residual ε_{k} from the Jacobian *J* and the observed state variables *y*_{k} and *y*_{k+1}. The noise strengths were then estimated from the root mean square of the residuals.

### Regression analysis

To directly assess compensatory adjustments between a perceived variable with respect to a controlled variable, linear regression analyses were conducted. The idea is that if a perceptual variable is used to control an action variable, then we expect a significant nonzero slope in a regression. This slope can be interpreted as a feedback gain in a proportional controller. The most insightful regression proved to be the one of the racket contact velocity at bounce *k* + 1 on the height error in bounce *k*. For comparison, the same regression analyses were conducted for model simulations that were performed by iterating the stochastic ball-bouncing map (*Eqs*. A*5*, A*6*, and A*7*) with the same set of parameters as for the autocorrelation analysis. For the noise strengths we took small values (*q*_{v} = 10^{−3} m/s, *q*_{t} = 10^{−4} s) because the slopes do not depend on the noise strengths.

### Statistical analyses

All dependent measures were subjected to a two-way mixed-effect ANOVA with group (between-subjects factor, 2 levels) and α (within subjects factor, 7 levels). For this ANOVA averages of the six trials for each α condition per participant were used (to satisfy the assumption of independent samples). Although the ANOVA gave significant results, the relatively small number of subjects may have violated or given insufficient support to the assumptions underlying the ANOVA. Further, due to taking the averages across six trials, the variance from these six trials was lost. Hence, to provide additional statistical tests of the group differences in lag-1 autocorrelation and regression slope, which are a core result of the study, permutation tests were conducted that did not rely on the assumption of independence of samples (Wassermann 2003). For this test, all 336 data points were pooled [6 (trials) × 7 (α) × 4 (subjects) × 2 (group)] and a bootstrapping method was applied: first, from the pooled data two data sets were sampled randomly (without replacement) to represent a new sample for each subject group. Second, for each sampled group a mean was calculated (over six repetitions, seven alpha conditions, and four actors). Third, the difference between the two groups served as the test statistic. Fourth, the three steps generating bootstrap data sets were repeated 1,000,000 times to yield a distribution of difference values. Fifth, the observed group difference was compared against the distribution of differences of the bootstrap data sets. The observed value was significant if it was not very probable. Probability was estimated by calculating the *P* value as the ratio between the number of values from the random permutations that exceeded the observed value and the total number of bootstrap data sets.

To compare the experimentally observed and theoretically predicted lag-1 autocorrelations and regression slopes, we used a *t*-test for each α condition. We could not use ANOVAs for this comparison because the model predicted only one value per group. The significance level for the *t*-test was 0.05, but since we ran seven tests, one for each α condition, we applied a Bonferroni correction to this significance level (0.05/7 = 0.007). We also used the permutation test to investigate the influence of experience on how much the observed values of the lag-1 autocorrelation differ from their values predicted from the ball-bouncing map. The procedure was identical to the one described earlier except that the difference between the observed lag-1 autocorrelations and the theoretical prediction was entered into the calculations.

In most figures, we plotted mean values averaged over trials and actors in a group (exception is the inherent noise level of Fig. 8 where the median was used). The error bars in all data plots denote the SEs over actors within a group (except Fig. 8; see figure caption).

## RESULTS

### Height error

To evaluate task performance both the means and SDs of the height error *HE* were compared across the seven α conditions for all participants. To assess whether the different experience of the eight participants was evident in performance their individual scores were assessed. The mean *HE* values for the two participant groups are listed in Table 1. The average values pooled over trials and participants were positive and showed a tendency to overshoot the target, which was larger for higher α values [*F*(6,36) = 24.54, *P* < 0.0001] (Table 1). The different experience with the task between participants did not lead to a difference in the mean *HE* itself. The apparent tendency that nonexperienced subjects have smaller mean *HE* values does not indicate higher accuracy; because *HE* is a signed error, and large alternations between positive and negative errors caused these seemingly smaller mean errors. Note this difference was also not significant due to the large variance of *HE*. Hence, a better reflection of the different experience levels was given by the SDs of *HE* shown in Fig. 3.

Consistent with subjects’ previous experience, their performance variability shows significant differences between the experts and novices. The ANOVA verified that experienced participants showed significantly higher variability, *F*(1,6) = 112.25, *P* < 0.0001. On the basis of this clear difference in variability, the participants were separated into two groups of four participants each for the following analyses. The interaction was also significant signaling larger separation for the smaller α values, *F*(6,36) = 2.94, *P* < 0.05. The main effect for α was also significant and showed that for higher α, the variability decreased, *F*(6,36) = 11.70, *P* < 0.0001. The decrease in variability at larger α conditions ran counter to Prediction 2 and will be examined in more detail below.

### Racket acceleration at impact

The average impact accelerations for the two participant groups per α condition are shown in Fig. 4*A* . Consistent with Prediction 1a and 1b experienced actors consistently showed negative values indicating that they exploited passive stability of the task for all α conditions, *F*(1,6) = 25.47, *P* < 0.0001. The range of impact accelerations fell within the boundaries predicted by the model (*Eq. 3*). The mean *AC* values of nonexperienced actors were close to zero, with the error bars indicating many trials with positive *AC*, especially for the lower α conditions. More precisely, 28% of the trial means of the nonexperienced actors were positive; in contrast, none of the trial means of the experienced actors were in the positive range. Additionally, an interaction was observed, *F*(6,36) = 2.41, *P* < 0.05.

Consistent with these results, the SDs of *AC* were significantly lower in experienced actors [*F*(1,6) = 65.46, *P* < 0.0001]. However, inconsistent with Prediction 2, the variability decreased with increasing α for both groups [*F*(6,36) = 80.57, *P* < 0.0001] (Fig. 4*B*).

### Racket amplitude and period

The means of racket amplitude and period within trials were calculated and then averaged across participants. Tables 2 and 3 show clear changes in the mean values of both variables across the α conditions: the racket amplitudes *A* decreased and the periods *T* increased with increasing α. Results for the 2 × 7 ANOVA for *A* rendered the two main effects significant [group: *F*(1,6) = 19.82, *P* < 0.0001; α: *F*(6,36) = 350.54, *P* < 0.0001]. The decrease in *A* with α can be understood from the physics of the task: the more “bouncy” the racket, the lower the racket velocity at impact necessary to bounce the ball to the fixed target height. The mean *T* values showed the opposite trend: both groups increased their periods for higher α values [*F*(6,36) = 53.06, *P* < 0.0001]. This increase in *T* is partly due to the increase in *HE* because larger ball amplitudes increase the time intervals between bounces. Further, participants tended to lower the racket impact position for the higher α values. This was possible given that participants were free to choose their preferred impact position while the target height was fixed.

The variability estimates of *A* and *T* also showed clear trends that resembled the pattern seen in both *HE* and *AC* variability results. In Fig. 5 the SDs of *A* showed a highly significant decrease with higher α [*F*(6,36) = 87.04, *P* < 0.0001] and experienced actors had lower variability values [*F*(1,6) = 33.36, *P* < 0.0001]. The interaction revealed that this difference was observed mostly for the lower α conditions [*F*(6,36) = 3.31, *P* < 0.05]. The periods showed the same pattern of results, only the interaction was no longer significant [group: *F*(1,6) = 83.27, *P* < 0.001; α: *F*(6,36) = 8.45, *P* < 0.0001]. This decrease in variability in performance measures is in contradiction to Prediction 2.

### Covariance functions

Figure 6 illustrates the autocorrelations between release velocities *ACF*(v-v) for 10 lags (bounces) both from the model and from a typical participant for the seven α conditions. As expected from the theoretical stability analyses shown by crosses, the correlations for lower α conditions were smaller and quickly converged to zero due to the higher stability of smaller α values. In contrast, higher correlation values for more lags indicated persistence of the decay of the deviations, indicative of a slower relaxation time. Higher α conditions showed a distinct pattern across the 10 lags with a change from positive to negative. Although all α conditions showed an exponential return of the correlation values to zero, higher α values showed an oscillatory pattern that is the result of the imaginary part of the eigenvalues of the Jacobian.

Comparing the empirical with the theoretical correlations a difference can be observed that becomes more prominent for larger α conditions. Particularly for α > 0.6 the theoretical predictions are outside of a band of 1SD of the observed mean, indicated by the bars. In contrast to the theoretical results, the experimentally observed autocorrelations decayed to zero faster than the model and beyond lag-2 the observed autocorrelations were rarely significantly different from zero. The oscillatory behavior displayed by the theoretical *ACF*s was only minimally reflected in the experimentally observed *ACF*s with the exception of α = 0.9. The autocorrelations of the times between impact, *ACF*(t-t), and the cross-correlations between release velocity and period, *ACF*(v-t) and *ACF*(t-v), revealed patterns similar to those of *ACF*(v-v) and were thus not shown. This similarity was caused by the strong coupling between the two state variables, ball release velocity, and time between contacts. These results are only partially consistent with Prediction 2.

Because the main differences between model and data were predominant at short lags, further analyses focused on lag-1 only. Figure 7 depicts *ACF*(v-v) for lag-1 only but compares the correlation values against α for experienced and nonexperienced participants separately. The lag-1 *ACF*(v-v) from the data increased with increasing α, corresponding to the model's prediction for both groups. Thus overall, Prediction 2 was supported by the data. However, although the *ACF* values of model and data matched closely for α ≤0.5, for higher α conditions the data exhibited smaller *ACF* values and the model values were no longer within the error band of the data. This indicates that in real bouncing fluctuations were compensated faster than in the model, probably due to active error compensation. Separate *t*-tests compared model and data for each α condition for both groups and supported this impression: for α >0.5 there were significant differences between the model and the data (all *P* < 0.001); for α ≤0.5, the differences were not significant (all *P* > 0.05).

Although the plots at first sight suggest similar behavior for both participant groups, a direct comparison of the performance data shown in Fig. 7*C* revealed that there were significant differences between groups. The 2 × 7 ANOVA rendered the comparisons between groups and α significant [group: *F*(1,6) = 11.32, *P* < 0.01; α: *F*(6,36) = 12.62, *P* < 0.0001]. To corroborate the significant difference between the small number of participants a permutation test was applied (Wassermann 2003). Results revealed that the observed difference of 0.0751 between the two groups was highly significant (*P* = 0.00003).

The results suggested that experienced actors deviated more from the model predictions than did nonexperienced actors. To test this observation directly, we calculated the difference of the observed lag-1 autocorrelations of all 336 data and the respective predicted values and applied the same permutation test as described earlier. The observed difference between the two groups was 0.0328, rendering *P* = 0.04. Hence, it can be concluded that experienced subjects deviated more from the predictions.

### Noise level

The noise strength in the velocity dynamics (matrix element 1,1 of matrix *Q* in *Eq*. B*2*) is plotted in Fig. 8 as a function of α, separately for the experienced and nonexperienced groups. As expected from the analysis of the residuals, the estimates of pure noise strength were more variable between trials than the performance measures. Thus medians and interquartile ranges rather than means and SDs were calculated. As Fig. 8 shows, noise strength was generally lower for the experienced group and decreased with increasing α. The 2 × 7 ANOVA revealed a significant interaction between group and α [*F*(6,36) = 2.61, *P* < 0.05]. Both main effects of group and α were significant [group: *F*(1,6) = 36.04, *P* < 0.0001; α: *F*(6,36) = 6.61, *P* < 0.0001]. The decrease in noise strength is in agreement with the decrease in height error. The small difference in noise strength for α = 0.9 is an observation that resembled the *HE* results where three of the four actors in the nonexperienced group obtained height error variabilities comparable to those of actors in the experienced group.

### Regression analyses

To uncover potential contributions of active control based on perceived errors more directly, regression analyses were performed of the racket velocity at contact *v* on the previous height error *HE* (see Fig. 1*B*). Regression slopes were negative for each actor in each α condition. *R*^{2} values varied between 0.2 and 0.3 with an increase of *R*^{2} with increasing α values. However, these findings cannot be taken as direct evidence for feedback control because open-loop stability of the model implies that similar compensatory relations exist due to the inherent dynamics that are the basis for the observed stability. Hence, the negative slopes had to be compared against the negative slopes predicted from the stochastic model.

As before, these model simulations were performed separately for each group and α value because the specific parameter values of individual participants have an influence on the relaxation behavior of the model. The resulting regression slopes for the participant groups are plotted in Fig. 9, *A* and *B* (as crosses) and were compared against the regression slopes obtained from the data (points with error bars). Across all α values participants showed more negative slopes than would have been expected from open-loop stability alone. For the nonexperienced group, significant differences were supported by *t*-tests performed separately for each α condition (all *P* < 0.0001). Similar results were found for the experienced group: all *t*-tests were significant with *P* < 0.0001 (α = 0.4, *P* < 0.05), except for α = 0.3, where *P* = 0.43. For the experienced group the differences were more pronounced, especially for the higher α conditions. These results mean that racket impact velocity was influenced and modulated by the height error of the preceding cycle.

Figure 9*C* directly compares the slopes of two groups. Experienced actors showed more negative slopes than those of novices, especially for the higher α conditions. The 2 × 7 ANOVA revealed both main effects were significant [group: *F*(1,6) = 6.70, *P* < 0.05; α: *F*(6,36) = 5.18, *P* < 0.0005]. A permutation test corroborated the significant result between the two groups. The mean observed difference between the groups was 0.0883, which rendered *P* = 0.00012 against the distribution obtained from the permutations. These results mean that experienced actors applied more active corrections, counter to Prediction 3b.

## DISCUSSION

The present study aimed to reveal the contributions of both open- and closed-loop components in the performance of a rhythmic ball-bouncing task. The foundation for this investigation is a mathematical model of the task defined by a racket interacting rhythmically with a ball in a gravitational field. It is assumed that the racket movements are generated by actors viewing and intercepting the ball trajectories with the goal to produce stable rhythmic performance. Hence, fine-grained analyses of the kinematics of racket and ball movements are expected to provide insight into the actors’ control strategies and the required underlying neurophysiological mechanisms. As such, the approach takes the actor–environment system as the window of analysis paying tribute to the fact that the task is defined over the actor in interaction with the ball that sets constraints on the generation of action (Warren 2006). An essential characteristic of the task as revealed by stability analyses of the model is that it can be performed with open-loop or dynamic stability if one condition is satisfied: the racket's trajectory should decelerate directly before and around the moment when contacting the ball assuming approximately sinusoidal movements of the racket. Dynamic stability implies that small deviations from periodic performance die out without requiring any error compensation in the actor's/racket's movements—which is the reason for the equivalent label of “passive stability” (Dijkstra et al. 2004; Schaal et al. 1996).

Empirical evidence that human actors execute the task with this passive strategy has been provided in a series of experiments where participants performed the task in two- and three dimensions and with different experimental equipment (Katsumata et al. 2003; Schaal et al. 1996; Sternad 2006; Sternad et al. 2000, 2001). However, in addition, two previous studies have demonstrated that actors also actively compensate for errors when they are caused by relatively large experimentally applied perturbations (de Rugy et al. 2003; Wei et al. 2007). More specifically, Wei et al. (2007) probed the relevance of the theoretical notion of the domain of attraction of the ball-bouncing map by different perturbation magnitudes. The relaxation times following different perturbations were a reflection of the domain of attraction, although generally shorter than predicted. Together with the observed negative impact accelerations, an indicator of passive stability, these results were interpreted as evidence for a blend of active control with passive stability. These results raised the question whether such intertwining of active control and passive stability is also present at steady state.

To pursue this question, the model of the bouncing ball was extended by a stochastic component. During steady state, the sequence of bounces constitutes a time series where the sequential structure of the values of the two state variables (ball release velocity and time between impacts) can be characterized by auto- and cross-correlations (Dijkstra et al. 2004). As analyses of the autocovariance function of the stochastic dynamics revealed, there is indeed a characteristic signature that depends on the coefficient of restitution of the ball–racket contact α. In particular, the model predicted that dynamic stability decreases with increasing α. Hence, the present study examined steady-state performance and manipulated the key parameter α. Given some previous evidence that the strategy of passive stability was not a trivial solution immediately adopted by all actors, we included participants with different levels of experience (Dijkstra et al. 2004; Morice et al. 2007).

### Open-loop stability

Consistent with the previous studies human actors, and in particular those with experience, displayed negative accelerations at impact. It should be emphasized that from a mechanical perspective this solution is not the most energy efficient one: hitting the ball at the moment of maximum velocity during the upward swing would yield the highest ball amplitude for a given racket movement. Decelerating toward contact is therefore not an intuitive or optimal solution from this perspective. However, based on our analyses it can be argued that it is the most efficient solution with respect to required computations. Although the actor needs to generate rhythmic racket actions with the required contact in an anticipatory fashion, he/she does not need to correct for every small deviation from periodicity or, equivalently, from the target height of the ball amplitude. The present data showed that less experienced participants did not adopt this strategy with the same consistency as did experts because the error bars include positive-acceleration values, indicating that some actors performed with positive values.

The key role of dynamic stability and its indicator acceleration at impact has also been highlighted in an extension of the bouncing-ball model into two dimensions. Ronsse and colleagues derived a model for two surfaces configured in a V-shape in the plane that bounce the ball between them by rhythmic antiphase movements. One key prediction from model analyses was that robustness of the task performance was determined by specific values of the acceleration at ball contact. This prediction was confirmed by human actors who performed the rhythmic task in an apparatus simulating this V-shaped paddle arrangement (Ronsse et al. 2006, 2008). These lines of work on uni- and bimanual ball manipulation are related to the work on passive dynamic walking where theoretical studies and implementation in actual robots have explored to what extent passive stability as provided by the mechanical interactions between the legs and the floor are exploited in biped walking (Adamczyk et al. 2006; McGeer 1990; Ruina et al. 2005). The philosophy behind this vibrant line of research is that understanding of passive mechanics is the foundation for the design of efficient control strategies to produce stable and adaptive walking patterns. This notion of stability and control is also the motivation for many other lines of research.

### Variability and autocovariance structure

Although these results again provide strong support for the presence of a strategy of open-loop stability, the variability results seem to contradict the model predictions. The variability of impact acceleration showed a marked decrease with α, that is, for conditions that were predicted to have less stability. This decrease in variability with higher α was echoed in other less specific variables such as the variability of height errors, racket amplitudes, and periods. That these results are not completely unintuitive is seen in the fact that in all variables novices show higher levels of variability, as is always the case in such comparisons of skill levels. If no further theoretical analyses were available, this decrease would be interpreted as higher stability in the “bouncier” racket–ball conditions. Nonetheless, this mismatch with the model predictions motivated more fine-grained analyses of variability and its temporal structure.

Analyses of the sequential structure of the state variables via auto- and cross-correlations in both model and data provided another window into assessing the relation between open-loop stability and variability. Performance measures like the SDs of height errors within a trial are relatively unspecific and constitute both the self-correcting properties, due to dynamic stability quantified by the eigenvalues of the Jacobian, and the internal noise. In contrast, the autocovariance function (*ACF*) is determined only by the dynamic stability measure. *ACF* analyses of the model yielded that the lag-1 autocorrelation of the state variables increased with increasing α, indicating a decrease of dynamic stability with α. The data on lag-1 autocorrelations were in qualitative agreement with this model prediction, although only for α conditions <0.6. In conditions with higher α subjects showed smaller lag-1 correlations. This mismatch with the predictions was seen in both groups, although it was more pronounced in the experienced actors. One viable explanation for this mismatch is that other than open-loop strategies are present. In particular, lower lag-1 values indicate faster relaxation processes, suggesting that active error corrections were present. It makes sense that such active control processes are more prominent for less stable conditions. It also is not unintuitive that experienced actors showed more of these subtle corrective processes.

Yet, before drawing further conclusions, the *ACF* analyses also provided a rationale for how to assess the amplitude of the internal noise *q* from the data. Given its independence from the autocorrelative structure, the noise *q* was quantified by the residual or unexplained variance of the *ACF* fits. Interestingly, the data showed that this noise amplitude decreased with increasing α and that skilled actors showed a generally lower noise level. Returning to the results on variability in performance measures these two results from the *ACF* analysis suggest that this decrease in performance variability at the higher and theoretically less stable α conditions is brought about by a decrease in the level of internal noise. This decrease in internal noise can be interpreted as a signature of active control. Because the amplitude of movement decreases with α and presumably the control signal also decreases we can interpret the observed decrease in internal noise as also consistent with the idea of signal-dependent noise (Harris and Wolpert 1998).

Another line of research that has used autocovariance analysis as a tool to understand human motor behavior is the work on rhythmic timing by Wing and Vorberg (1996). Although the task of bouncing a ball to a target height shares some similarities with finger tapping—in that both are rhythmic movements with a collision—there are critical differences in the task and the analysis approach. First, in rhythmic tapping the task is only about the time between contacts; the contact force and spatiotemporal aspects of the trajectory are typically not prescribed nor within the purview of the analyses (Balasubramaniam et al. 2004; Sternad, et al. 2000). Hence, corrections for an early or late tap only require that the actor catches up or slows down on the next tap. In contrast, in ball bouncing the effect of an early or late contact depends on the impact acceleration and has a consequence on the subsequent ball amplitude. Moreover, how hard an actor hits the ball is crucial to success in the task. Second, in our approach to ball bouncing the modeled mechanics of the task is hypothesized to constrain the solutions chosen by the actor. The trajectory is analyzed to examine whether the actor makes use of the dynamics of the task. The autocovariance structure is derived from a model of the task dynamics and the parameters of this task dynamics are the window into the CNS. In contrast, in tapping the observed autocovariance structure is interpreted directly in terms of the underlying neural mechanisms because the dynamics of the finger-surface contact is not analyzed.

### Active error corrections

To strengthen the inference of active control we scouted the data for direct evidence for active error compensation. Although many variables can be considered as candidates for error correction the most plausible perceived variable was the height error of the ball amplitude. Because subjects were instructed to bounce the ball to a visible target height with every contact, the height error provided direct information to adjust the racket trajectory at the subsequent ball contact: A higher than desired ball amplitude would necessitate a lower racket velocity at the next contact to compensate for the increased ball velocity at contact. Hence, the perceptual variable height error and the action variable racket velocity at contact were hypothesized to be negatively related to achieve active error correction. This choice of variable was also motivated by results from research on juggling and other ball skills, where the apex of the ball amplitude was shown to carry the most salient information for controlling the next ball catch (van Santvoord and Beek 1994).

The regressions of the action variable on the perceptual variable indeed revealed significant negative slopes that can be interpreted as evidence for servo-control. Again, the theoretical analyses provided information to qualify this result: because the dynamically stable stochastic system also displays such negative slopes, only the comparison between observed and predicted slopes allowed us to draw inferences about active and passive corrections. The slopes in the data were more negative for all α conditions and thereby confirmed the presence of active compensation. Consistent with the better (i.e., less variable) performance in experienced actors, these individuals showed more pronounced compensatory relations. Additionally, these deviations were slightly higher for the higher α conditions. Together with the deviations of lag-1 results for higher α conditions, these results clearly indicate that subjects apply more active control for less stable conditions.

### Neurophysiological mechanisms

The concept of stability has been helpful in much neurophysiological research that aimed to shed light on the mechanisms of the peripheral nervous system. In many elegant studies it was shown how short- and long-latency reflex mechanisms operate to stabilize position and trajectories (Nichols and Houk 1976; Popescu et al. 2003). Such mapping of stability to explicit neural mechanisms is not so easy in more complex tasks such as ball bouncing. Dynamic stability in ball bouncing not only is maintained through servo-control mediated by reflex mechanism, but also is established in a corrective as well as anticipatory fashion, by both active and passive components, as we try to tease apart.

The active control components map more easily onto the traditional error-compensation mechanisms: by correlating perceived height error to subsequent racket velocity at ball contact, we demonstrate that perceived information leads to a modification of the efferent signal controlling the racket probably in a transcortical sensorimotor loop (the time between the height error at the apex and the following racket–ball contact is 350 ms, which allows for visually controlled corrections). However, as we try to show, this is only one part of the complex control system that is intertwined with other mechanisms that ensure passive stability, where error corrections are inherent in the dynamics. These other mechanisms operate in a more anticipatory or task-oriented fashion and are not (yet) easily mapped onto neural substrate.

That corrective mechanisms are not necessarily a function of a fixed set of neural mechanisms has been demonstrated in a recent study on bimanual reaching by Diedrichsen (2007). When both hands perform parallel reaching movements to two separate targets and one hand is perturbed, short-latency responses are observed in the perturbed hand. In contrast, if the task is changed such that the two hands operate a single cursor toward a single target, both hands showed compensatory reactions to the same single-hand perturbation (within 190 ms). These results give evidence that responses to perturbation depend on the task and the goal.

The topic of stabilizing a task goal rather than a fixed position is also a core assumption of optimal feedback control (Todorov 2004; Todorov and Jordan 2002). In this theoretical framework the controller/CNS ensures that a task goal is stabilized; deviations from the optimal solution are corrected for only when they affect the task outcome, whereas deviations that do not affect the end result are not compensated for (minimum-intervention principle). That this proposition opens a new route for identifying neural underpinnings, in particular for the role of the motor cortex, has been recently argued by Scott (2008). He points out that the research over the last 20 years that has been dedicated to map neuronal activations in primary motor cortex (M1) (both single neurons or populations) to selected variables in execution (e.g., direction, load, etc.) has led to a plethora of interesting but also seemingly contradictory results. To give more cohesion to these investigations on M1 he proposes that the search should be for principles of control, rather than for specific variables that are controlled. Explicitly, he suggests the framework of optimal feedback control.

In light of this discussion, the current research documents that the CNS ensures that the task goal of rhythmically bouncing a ball is attained, partly by active error-correcting and partly by more anticipatory processes leading to dynamical stability. It remains to be a fascinating question which processes of the CNS achieve this dynamical stability, a question that should be pursued in more detail in future work.

In sum, the theoretical and empirical results revealed that human actors use a subtle web of active and passive strategies in the performance of a rhythmic task. Humans are sensitive to the dynamical stability of the task but they also apply active modulations when necessary. Given that experienced actors exploit the dynamic stability more but also apply more of these active modulations it can be inferred that this web of open-loop and closed-loop strategies is the result of adaptation and learning processes. As both theoretical and empirical results showed, it is far from evident that variability is the inverse of stability, as often assumed in the absence of more detailed predictions. Similarly, far from evident is that compensatory relations between perceptual and action variables always speak to active error corrections. Dynamic stability within a task can give rise to such compensatory behavior without reflecting active closed-loop strategy. Only the dialogue between model and data allowed more differentiated explanations of such observations.

## APPENDIX A

The model is based on the same deterministic ball-bouncing map previously derived by different studies that revealed its period-doubling route to chaos (Bapat et al. 1986; Guckenheimer and Holmes 1983; Tufillaro et al. 1992). However, all studies on the model's application to the ball-bouncing task have examined only the simplest behavior of this map—fixed-point stability implying periodic bouncing. We therefore review details pertaining only to period-1 behavior.

The map is based on three assumptions: *1*) between bounces the ball follows ballistic flight under the influence of gravity, *2*) the impact is instantaneous with the coefficient of restitution capturing the energy loss at impact (see footnote 1), *3*) the racket movement is assumed to be a pure sine wave with a fixed amplitude and frequency. Figure 1*A* shows the mechanical model with the planar surface and the ball and Fig. 1*B* shows exemplary time series of racket and ball trajectories (with added noise). Based on these assumptions, the deterministic ball-bouncing map can be derived as (A1) (A2)

This is an implicit map with state variables *v*_{k}, the ball velocity just after the *k*th impact, and *t*_{k}, denoting the time of the *k*th impact. Parameters are *g*, the acceleration of gravity; α, the coefficient of restitution; and ω_{r} and *a*_{r}, the angular frequency and amplitude of racket movement, respectively.

The map has a period-1 attractor that is stable when the impact acceleration *AC* fulfills the following constraints (A3) For the selected α values used in the experiment, the stable range for *AC* is between −12.65 m/s^{2} and 0 (α = 0.3) and −9.83 m/s^{2} and 0 (α = 0.9).

The assumptions of ballistic flight and instantaneous impact are obeyed by design in the virtual reality setup, although the assumption of a sinusoidal racket trajectory is not strictly obeyed in experimental data. Human actors move the racket in approximately sinusoidal fashion but with a steeper slope preceding the impact position that has no simple mathematical expression (Sternad et al. 2001). However, because only position and velocity of the racket at the ball contact determine the ball trajectory, an equivalent sinusoid is created that is identical to the actual waveform around the time of the impact. Its equivalent frequency is calculated from the period between bounces; its equivalent amplitude is calculated from the stationary phase of impact θ as (A4)

Because we apply no external perturbations we assume that the stationary phase can be estimated from the mean racket impact phase within a trial. This equivalent amplitude is larger than the observed amplitude but its slope in its upward phase around the impact closely matches the observed racket waveform.

As shown in Dijkstra et al. (2004) the deterministic model predicts different degrees of stability with increasing values of α: higher α values are associated with lower stability. However, this prediction is difficult to test in practice. One would need to administer small perturbations that do not lead to active corrections. However, such perturbations are drowned in performance noise (see the data for the perturbations P − 1 and P + 1 in Figs. 9 and 10 of Wei et al. 2007). It is therefore better to look at autocorrelations that summarize behavior over many small perturbations. This, however, requires a stochastic model to make predictions. Therefore we proceeded to extend the deterministic model with a stochastic element.

The stochastic ball-bouncing map is a straightforward extension with a small caveat for the implicit time map (*Eq*. A*2*): we introduce the intermediate state variable τ_{k+1} representing the noiseless time of the next impact. Using this variable, the stochastic ball-bouncing map reads (A5) (A6) (A7) where ξ_{k} denotes an internal Gaussian white noise process with zero mean and unit variance. The stochastic model has two additional parameters: the noise strengths of the velocity and time dynamics denoted by *q*_{v} and *q*_{t}, respectively. Gaussian white noise was chosen based on the idea that the noise originates from many independent contributions from within the CNS. For example, these contributions could originate at the planning stage, from the transmissions of motor commands, or in the end effectors themselves.

The assumption of additivity of the noise is admittedly ad hoc and mainly chosen for mathematical simplicity. As an analysis of the observed noise strength *q*_{v} shows (Fig. 8), the results seem more consistent with a noise strength that increases with increasing impact velocity, which would suggest a multiplicative noise term. However, such a noise model is left for future analysis.

The model has six parameters for which, in principle, sensitivity analyses may be warranted. However, our model does not have a clear separation of states and parameters. Two of the parameters, *g* and α, are fixed by design in our computer-generated setup. For the two noise strength parameters a partial sensitivity analysis was performed in the theoretical study (Dijkstra et al. 2004) where we showed that the observed noise strengths are such that the linear approximation of *Eq*. B*1* is accurate. The parameters amplitude and frequency of the equivalent sinusoid describe the racket motion of the subject. These parameters vary across subjects and may even vary within a performance and across trials. Therefore a traditional sensitivity analysis is not feasible. Yet, one can view the main thrust of the study, looking for signatures of active control, as a sensitivity analysis of the assumption of a fixed sine wave.

## APPENDIX B

To derive predictions about how noise percolates through the system the correlation structure of the state variables was calculated. To this end the stochastic map was linearized around the period-1 attractor (denoted by ṽ and t̃) and the state of the system was denoted by *y*_{k} and the noise by ε_{k} with the linearized stochastic dynamics (B1) where *J* denotes the Jacobian of the ball-bouncing map (see *Eq. 39* in Bapat et al. 1986 and *Eq. 9* in Dijkstra et al. 2004). The Jacobian quantifies the stability properties of the ball-bouncing map. In particular, the degree of stability can be quantified by the inverse of the eigenvalues (Scheinerman 1995). The derivation of the vector autocovariance function *R*_{y} is standard procedure for this vector autoregressive system of order 1 (Kendall and Ord 1989; Shumway and Stoffer 2001). The autocovariance function at lag-0 of the two state variables, *R*_{y}(0), could be calculated from the Lyapunov equation (B2) where *Q* denotes the autocovariance matrix of the noise ε_{k}: *Q* = diag (*q*_{v}^{2}, *q*_{t}^{2}). For nonzero lag-1 the autocovariance function follows from (B3)

The Jacobian depends on all four parameters of the ball-bouncing map. However, Dijkstra et al. (2004) showed that in many cases the absolute values of the eigenvalues of the Jacobian depend only on α. Exceptions occur close to the boundaries of stability. Specifically, the absolute values of the eigenvalues of the Jacobian equal α. Thus the autocovariance function in *Eq*. B*3* can be approximated as (B4)

There is one more caveat in using the preceding equations for comparing theory and experiment: because state variable *t*_{k}, the time of impact, does not have a constant value at stationary state and continuously increases: *t*_{k+1} = *t*_{k} + *T*, where *T* = 2π/ω_{r}. Although the period *T* is known in the model, it fluctuates in the data due to noise. This created a problem when calculating the predicted autocovariance functions from data. To avoid this problem, we used the time between impacts δ_{k} = *t*_{k} − *t*_{k−1} as state variable. We denoted the new vector state variable by *z*. The autocovariance function based on this new state variable *R*_{z}(*l*) is related to the old one *R*_{y}(*l*) as follows (B5) where *R*^{i,j} denotes matrix element *i*, *j*. Since the (1, 1) matrix element of the autocovariance in unchanged by this change of state variable we used the lag-1 autocovariance of release velocity as quantification of stability. A specific prediction from this relation is that stability decreases as the lag-1 autocovariance of release velocity increases. Importantly, because the lag-1 autocovariance of release velocity equals α, the model predicts a decreasing stability with increasing α. This prediction is consistent with and supports the prediction derived from the deterministic model.

## GRANTS

This research was supported by National Science Foundation Grant BCS-0450218, National Institute of Child Health and Human Development Grant R01-HD-045639, and Office of Naval Research Grant N00014-05-1-0844.

## Footnotes

↵1 The coefficient of restitution can be illustrated with the following example: assuming the ball prior to contact to have a velocity of −2 m/s (relative to the racket) and the coefficient of restitution to be 0.5, the ball would bounce off the racket with a positive velocity of +1 m/s.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society