Pauline Maurice
email: pauline.maurice@polytechnique.org

Neville Hogan

Dagmar Sternad

Predictability, Force and (Anti-)Resonance in Complex Object Control

Keywords: motor skill, rhythmic movements, object manipulation, prediction, interaction force, impedance

Manipulation of complex objects as in tool use is ubiquitous and has given humans an evolutionary advantage. This study examined the strategies humans choose when manipulating an object with underactuated internal dynamics, such as a cup of coffee. The object's dynamics renders the temporal evolution complex, possibly even chaotic, and difficult to predict. A cart-and-pendulum model, loosely mimicking coffee sloshing in a cup, was implemented in a virtual environment with a haptic interface.

Participants rhythmically manipulated the virtual cup containing a rolling ball; they could choose the oscillation frequency, while the amplitude was prescribed. Three hypotheses were tested: 1) humans decrease interaction forces between hand and object; 2) humans increase the predictability of the object dynamics; 3) humans exploit the resonances of the coupled object-hand system. Analysis revealed that humans chose either a high-frequency strategy with anti-phase cup-and-ball movements or a lowfrequency strategy with in-phase cup-and-ball movements. Counter Hypothesis 1, they did not decrease interaction force; instead, they increased the predictability of the interaction dynamics, quantified by mutual information, supporting Hypothesis 2. To address Hypothesis 3, frequency analysis of the coupled hand-object system revealed two resonance frequencies separated by an anti-resonance frequency. The low-frequency strategy exploited one resonance, while the high-frequency strategy afforded more choice, consistent with the frequency response of the coupled system; both strategies avoided the anti-resonance.

Hence, humans did not prioritize small interaction force, but rather strategies that rendered interactions predictable. These findings highlight that physical interactions with complex objects pose control challenges not present in unconstrained movements.

Introduction

Using tools has been essential in human evolution, and a large variety of tools now enhance and augment our daily actions. Tool-supported actions range from the simple swinging of a hammer and cutting meat with a knife to more complex or exotic actions, such as eating spaghetti and cracking a whip. The latter tasks are challenging and require practice because the objects themselves, spaghetti and whip, are flexible hence underactuated, i.e. have internal degrees of freedom that are not directly controlled by the user.

Another seemingly mundane example is carrying a cup of coffee: the human manipulates the cup that, in turn, exerts a force on the coffee that exerts forces back on the cup and the hand. Complex interaction forces arise between the hand, the cup and the coffee. Despite this complexity, humans are extremely skilled at interacting with such underactuated objects. Our understanding of how humans achieve such dexterity is still limited and becomes an ever-growing barrier to current developments in prosthesis control, brain-machine interfaces and robotic rehabilitation.

Despite the abundant literature on the control of goal-directed upper-limb movements, most studies have focused on free movements without physical interaction, such as reaching and pointing [START_REF] Flash | The coordination of arm movements: an experimentally confirmed mathematical model[END_REF]; Bhushan and Shadmehr 1999; [START_REF] Krakauer | Independent learning of internal models for kinematic and dynamic control of reaching[END_REF][START_REF] Sabes | The planning and control of reaching movements[END_REF], or interactions with rigid objects, such as grasping with isometric grip forces [START_REF] Flanagan | The role of internal models in motion planning and control: evidence from grip force adjustments during movements of hand-held loads[END_REF][START_REF] Fu | Coordination between digit forces and positions: interactions between anticipatory and feedback control[END_REF]. The control of "complex objects", which we define as objects with underactuated internal dynamics, i.e. non-rigid objects, has been largely ignored. The few studies that examined the control of complex objects have focused on the two classic control models of balancing a pole and manipulating a linear mass-spring system. For balancing a pole one needs to stabilize an inherently unstable inverted pendulum. Based on kinematic measurements and mathematical modeling, different mechanisms have been suggested, such as intermittent, continuous or predictive control, with forward or inverse models [START_REF] Mehta | Forward models in visuomotor control[END_REF]; [START_REF] Gawthrop | Human stick balancing: an intermittent control explanation[END_REF][START_REF] Insperger | Acceleration feedback improves balancing against reflex delay[END_REF]. Another set of studies on the inverted pendulum system focused on noise and delays to distinguish between the continuous vs. intermittent nature of control [START_REF] Cluff | Dynamical structure of hand trajectories during pole balancing[END_REF][START_REF] Milton | The delayed and noisy nervous system: implications for neural control[END_REF][START_REF] Milton | Delayed pursuit-escape as a model for virtual stick balancing[END_REF]. A linear mass-spring system has served as a model to examine optimization criteria in human control, such as generalized kinematic smoothness (Dingwell et al. 2014), effort and accuracy [START_REF] Nagengast | Optimal control predicts human performance on objects with internal degrees of freedom[END_REF], or minimum acceleration with constraints on the center of mass [START_REF] Leib | Minimum acceleration with constraints of center of mass: a unified model for arm movements and object manipulation[END_REF]. Two studies compared the contributions of visual and haptic feedback and their results highlighted the essential role of haptic feedback over visual feedback in controlling the object [START_REF] Huang | Visual and haptic feedback contribute to tuning and online control during object manipulation[END_REF]; [START_REF] Danion | The role of haptic feedback when manipulating nonrigid objects[END_REF]). Lastly, another set of studies looked at the compression of a buckling spring, modeling the buckling behavior with a subcritical pitchfork bifurcation of the nonlinear dynamic system, including integration of multi-sensory information with different time delays [START_REF] Venkadesan | Manipulating the edge of instability[END_REF]; Mosier et al. 2001).

All these studies examined point-to-point movements, or short sequences of discrete movements, in which the full complexity of the system's dynamics may not yet be fully manifest. A more extended continuous interaction may reveal more of the challenges arising from complex underactuated dynamics.

For instance, when a system is near an anti-resonance frequency, its evolution is very sensitive to small changes in the input, rendering the system's behavior chaotic, and essentially unpredictable in the longer term. Such small perturbations readily arise from the fact that human movements are intrinsically variable. This presents a problem for the widely-held assumption that humans rely on internal models of the manipulated object to select and execute a movement policy [START_REF] Flanagan | Control strategies in object manipulation tasks[END_REF]Dingwell et al. 2012[START_REF] Danion | The role of haptic feedback when manipulating nonrigid objects[END_REF]. How can humans learn an internal model of a complex underactuated object that has a potentially unpredictable temporal evolution? How can humans control the behavior of such objects? Relying on feedback control is largely insufficient for the manipulation of objects with complex dynamics due to neural transmission delay. Despite these challenges, humans skillfully manipulate complex objects of all degrees of complexity. How humans achieve this is an open question.

Extending previous work by Sternad and colleagues (Hasson et al. 2012a;[START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF]; [START_REF] Sternad | Control of intermittent and continuous objects[END_REF]Hasson 2016, Bazzi et al. 2018), this paper investigates continuous manipulation of an underactuated object with nonlinear internal dynamics. The task of moving a bowl-shaped cup with a ball inside was implemented in a virtual environment, using a cart-and-pendulum model to mimic the ball rolling in the moving cup. Notably, one of our previous studies demonstrated that the continuous evolution of this system shows features of deterministic chaos [START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF]). Using mathematical modeling and simulation of the task dynamics, this previous study examined the strategy that humans adopt when manipulating this complex object in continuous rhythmic fashion. Moving at an imposed frequency, participants chose movement amplitudes that made the interaction easier to predict.

Counter to expectation, interaction force and smoothness were not minimized.

The present study examined the same task, but extended the question in two ways. First, rather than imposing a frequency for the oscillatory movement, the present study prescribed the movement amplitude, leaving frequency free to choose. The task of choosing a frequency gave rise to new behaviors and new questions, because the resonance structure of the system may now play a significant role in the choice of strategy. Second, we extended the modeling of human control by including the mechanical impedance of the hand. The previous study on the same system only considered the dynamics of the cartand-pendulum system [START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF]). However, the object is in continuous interaction with the human, whose neuromechanical properties are likely to influence the cart-and-pendulum dynamics.

Therefore, this study introduced a simplified model of hand mechanical impedance interacting with the cart-and-pendulum system.

Several studies on unconstrained movements have demonstrated that humans tend to move in a way that minimizes physical effort (e.g. Alexander 2000; [START_REF] Prilutsky | Optimization-based models of muscle coordination[END_REF]. Extending these findings to the manipulation of complex underactuated objects, our first hypothesis is that humans seek to minimize the effort, or specifically the interaction force (Hypothesis 1). We assessed this hypothesis by quantifying the root-mean-squared value of the interaction force between the object and the hand. However, while demonstrated for free movements, this principle may become less prominent when the manipulated object presents additional challenges, specifically when it develops increasingly erratic behavior that becomes hard or impossible to predict. Therefore, we also tested the hypothesis that humans adopt strategies that make the hand-object interaction more predictable (Hypothesis 2). When interactions are predictable it is easier for humans to anticipate the object motion and hence the force arising from the object's internal dynamics. Anticipating this "perturbing" force, subjects can directly generate the appropriate interaction force to achieve the desired movement. Conversely, unpredictable object behavior requires continuous correction and adaptation of the hand movement, which may be tiring, both physiologically and cognitively. Predictability of the object dynamics may therefore obviate computational effort and afford simpler internal models to guide feedforward control. We assessed predictability by quantifying mutual information between the hand-cup interaction force and the object kinematics.

Addressing Hypotheses 1 and 2 rendered insight into human movement strategies (what do humans optimize), but they did not inform how humans achieved these strategies. Such explanation required closer analysis of the object dynamics. Numerous studies on rhythmic movements have provided evidence that resonance properties of the limbs or the object influence behavior. For example, in walking, the preferred stepping frequency maps onto the resonance frequency of the leg modeled as a simple pendulum [START_REF] Holt | The force-driven harmonic oscillator as a model for human locomotion[END_REF]). A study of infants in a "jolly jumper" showed that infants tune into the resonance frequency of the jolly jumper [START_REF] Goldfield | Infant bouncing: The assembly and tuning of action systems[END_REF]). Rhythmically swinging hand-held pendulums of different mass and length has demonstrated that humans have a tendency to oscillate at the natural frequency of the hand-pendulum system [START_REF] Yu | Task-effector asymmetries in a rhythmic continuation task[END_REF]. One main advantage of moving at the resonance frequency is its energetic efficiency: in oscillatory systems at resonance, the ratio between the amplitude of the movement output and the force input is maximal. Another feature of oscillating at resonance has been shown by [START_REF] Goodman | Advantages of rhythmic movements at resonance: minimal active degrees of freedom, minimal noise, and maximal predictability[END_REF] in a study on rhythmic limb movements. Time series analysis using phase space embedding revealed that the trajectories became more predictable when oscillating at resonance. However, that study focused on pendular limb movements, and the applicability of its findings to the manipulation of underactuated objects is unclear. We therefore tested an additional hypothesis that in complex underactuated object control, humans exploit the resonance structure of the manipulated object (Hypothesis 3). As the analyses showed, the manipulated object together with the hand not only had one, but two resonance frequencies separated by an anti-resonance frequency, a structure that will aid in interpreting the results.

In the experiment, participants manipulated a virtual cart-and-pendulum system at their preferred frequency with the movement amplitude prescribed. To evaluate the strategies that humans adopted we mathematically examined the cart-and-pendulum system coupled to a simple model of hand impedance. This model-based analysis allowed us to assess alternative execution strategies, i.e. different values of frequency and hand impedance that could be used to perform the task. Interaction forces and the degree of predictability were calculated both experimentally and in simulation. Comparison of human behavior with the mathematically derived results showed that participants did not minimize interaction force, but favored strategies with high predictability. In addition, frequency analysis of the coupled object-hand system showed that the degree of predictability was closely related to the resonance and anti-resonance frequencies of the system.

Behavioral Experiment

Participants

Ten young adults with no self-reported neuromuscular pathology volunteered for the experiment (mean age = 24.3±1.8 yrs). All participants performed the task with their dominant hand. They were naive to the purpose of the study and gave written informed consent before the experiment. All procedures were approved by the Northeastern University Institutional Review Board.

The Virtual Task

To test the three hypotheses, a virtual task mimicking the manipulation of a bowl-shaped cup with a ball inside was developed. Importantly, this system is underactuated, since moving the cup causes movements of the ball, which simultaneously exerts forces on the cup: the person moving the cup has to take into account these indirectly-controlled forces to obtain the desired movement of the cup. A simplified model of a cup-and-ball was simulated in a virtual environment with visual and haptic feedback via a robotic manipulandum. Participants were asked to move this virtual cup rhythmically between two specified targets, but were allowed to choose their preferred frequency.

The Mechanical Model

Similar to (Hasson et al. 2012a(Hasson et al. , 2012b;;[START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF][START_REF] Sternad | Predictability and robustness in the manipulation of dynamically complex objects[END_REF], the cupand-ball system was modeled as a ball sliding in a semi-circular cup (Fig 1A). The cup motion was limited to one direction in the horizontal plane, without any friction. Under the assumption that the ball does not roll, but only slides without friction between the cup and ball, the cup-and-ball system was mathematically equivalent to an undamped pendulum attached to a moving cart (Fig 1B). The ball corresponded to the pendulum bob, the cup's horizontal position corresponded to the cart position, and the arc of the cup corresponded to the pendulum's semi-circular path. With this simple model, the full dynamics of the task could be computed more easily, without sacrificing the essential elements of the dynamics: underactuated and nonlinear. Hence, the equations of the cart-and-pendulum motion are

(𝑚 𝑐 + 𝑚 𝑝 ) 𝑋 ̈= 𝑚 𝑝 𝑑 [ 𝜃 2 ̇sin 𝜃 -𝜃 ̈cos 𝜃] + 𝐹 𝑖𝑛𝑡𝑒𝑟 = 𝐹 𝑏𝑎𝑙𝑙 + 𝐹 𝑖𝑛𝑡𝑒𝑟 𝜃 ̈= - 𝑋 d cos 𝜃 - 𝑔 𝑑 sin 𝜃 ( 1 )
where X is the cart position, θ is the pendulum angle, Finter is the force applied by the human on the cart, and Fball is the force applied by the pendulum (the ball in the conceptual model) on the cart. Parameters of the system are the mass of the cart mc, mass of the pendulum mp, the pendulum length d, and the gravitational acceleration g. The following values were used: mc = 2.40 kg, mp = 0.60 kg, d = 0.45 m. These values were chosen because they rendered resonance and anti-resonance frequencies of the system that were well within human motor capacities and within reach of participants. The cart and pendulum masses were chosen to make the object light enough to avoid fatigue. The ratio of cart and pendulum masses was set to make the underactuated internal dynamics a prominent feature, i.e. participants clearly felt the forces generated by the ball. For lighter ball masses, the cart-and-ball system approximated a rigid object.

Apparatus and Data Acquisition

The dynamics of the cup-and-ball system were simulated in a virtual environment (Fig 2). Participants were seated on an adjustable chair in front of a screen and interacted with the virtual environment via a 3-degree-of-freedom robotic manipulandum (HapticMaster®, Motekforce, Amsterdam, Netherlands) (Van der Linde and Lammertse 2003). The force applied by the participants on the handle of the robotic arm (Finter in Eq 1) controlled the position of the virtual cup (X in Eq 1). The movements of the robotic arm were restricted to horizontal translations parallel to the participant's frontal plane to ensure a onedimensional motion of the cup as in the model. Participants felt the interaction force (system inertia and ball force Fball in Eq 1) via the force feedback provided by the robotic manipulandum. A custom-written C++ program based on the HapticAPI (Moog FCS Control Systems) computed the ball kinematics and controlled the virtual display as well as the force feedback. of the task in the virtual environment: the robotic manipulandum provided haptic feedback of the mechanical interaction with the object, while the behavior of the system was displayed online on the back-projection screen.

The physical model used the distances shown on the figure, while the distances displayed on the screen were multiplied by a factor of 4 for visibility. The cup displayed was 7.5 times smaller than the physical arc determined by the length d of the pendulum. B: A participant using the HapticMaster to interact with the simulated cup-andball system. The position of the cup was controlled by the position of the end-effector of the robot.

The cup and ball movements were displayed on a 2.40 m × 2.40 m back-projection screen located 2.15 m in front of the participants. The display consisted of two green rectangular targets on a horizontal line delimiting the displacement of the cup; a yellow semi-circle represented the cup and a small white circle represented the ball (Fig 2). Although the cup was only displayed as a semi-circle, there was no restriction on the ball angle and the pendular rotations could exceed 90° without the ball escaping the cup. The visual translation of the cup was 4.0 times the physical displacement of the manipulandum. The cup displayed on the screen was 7.5 times smaller than the physical dimension of the cup (set by the pendulum length d), in order to have plausible dimensions and fit the display. The force applied by the participants on the robotic arm (Finter), the cup kinematics (position X, velocity 𝑋 ̇, and acceleration 𝑋 ̈) and the computed ball kinematics (angular position θ, angular velocity 𝜃 ̇, and angular acceleration 𝜃 ̈) were recorded at 120 Hz.

Experimental Task and Instructions

Participants were asked to move the cup rhythmically between two targets located at a horizontal distance of 16.5 cm from one another (physical distance between the center of each target, Fig 2A). Participants were instructed to place the cup within the target rectangle at each excursion, so movement amplitude was prescribed. However, the scaled cup was 3 cm wide, while each target was 4.5 cm wide; the peakto-peak excursion of the physical cup oscillation could therefore range from 15 to 18 cm and still satisfy the task. This tolerance gave participants some leeway to develop their preferred motion. Further, participants were told that they could freely choose their frequency of oscillation and that they could change it throughout the experiment to arrive at their most preferred frequency. Even though participants did not receive explicit restrictions on the movement frequency, a demonstration of the task by the experimenter and the emphasis to "move rhythmically" discouraged them from extremely slow movements. Note that people do not necessarily prefer to move as slowly as possible, even though this may save effort (Van der Wel et al. 2010, [START_REF] Park | Moving slowly is hard for humans: Limitations of dynamic primitives[END_REF]. No instruction was given regarding the position of the ball within the cup, but participants were informed that the ball could not escape the cup (i.e. the behavior was that of a pendulumattached with a stringrather than that of a loose ball).

However, due to the haptic feedback provided by the manipulandum, participants could not ignore the movement of the ball: the ball movement affected the cart movement, as in a real system, and participants felt and saw it. Note that this experimental design intentionally refrained from specifying a single optimal task performance, but rather aimed to give insight into what participants preferred to do, especially after some exploration and practice.

The experiment consisted of 5 blocks of 10 trials each. Each trial lasted 45 s. The trials within a block were separated by a 15 s pause, and the blocks were separated by a break of several minutes. At the beginning of each trial, the cup was positioned at the center of the left target, and the ball rested at the bottom of the cup.

Data Analysis

As the task could be achieved by multiple solutions, i.e. it had redundancy, we distinguished between execution and the outcome or result of the movement. Performance was quantified by variables that fully described the kinematics of the system, i.e. amplitude and frequency of cart and pendulum, while the outcome was quantified by the task or result variables interaction force, predictability and resonance.

Result variables are metrics that explicitly tested the hypotheses.

Task Performance and Kinematic Variables:

The task instructions elicited trajectories close to a sinusoid, therefore the movements of the cart (cup) were characterized by the amplitude Ak and the frequency fk of each cycle k (i.e. each back-and-forth movement). The cart amplitude Ak was defined as the half-distance between the minimum and the maximum of the cart position during cycle k. The cart period Tk was defined as the time between two successive maxima of the cart position; the oscillation frequency was fk = 1/Tk. In addition, we quantified the relative phase between the cart and pendulum movements by computing the time lag that maximized the cross-correlation between the time-series of the cart position and pendulum angle. The resulting time lag was then converted into relative phase.

In order to detect the extrema in the cart position, the difference between successive data points, i.e. velocity, was computed. Extrema were detected as those values where the sign changed. In order to ensure robust detection of the cart extrema, the cart position data were smoothed with a zero-phase-lag, fourth-order, low-pass Butterworth filter with a 3 Hz cut-off frequency. Note that this smoothing was used only for detecting the extrema.

Result Variables: Hypothesis 1 -Minimize Interaction Force:

The net force required to perform the task was estimated by the root mean square of the continuous interaction force RMSF

𝑅𝑀𝑆𝐹(𝐹 𝑖𝑛𝑡𝑒𝑟 ) = 1 𝑇 ∫ 𝐹 𝑖𝑛𝑡𝑒𝑟 2 (𝑡)𝑑𝑡 𝑇 0 (2)
where T is the duration of the trial. Note that this hypothesis is about the hand-cart interaction force and not the overall force exerted by the participants. In particular, muscular effort was not evaluated.

Hypothesis 2 -Maximize Predictability: Predictability is a mathematical concept that can be operationalized in several ways. We opted to characterize the degree of predictability of the object dynamics by the mutual information between the input and the output of the system, i.e. the cart trajectory and the interaction force Finter. Mutual information is an information-theoretic metric that quantifies the statistical dependency between two variables, and thereby quantifies how much knowing one of the variables reduces the uncertainty about the other. High mutual information indicates a small degree of uncertainty (Cover and Thomas 2012). In the present context, mutual information quantifies the degree to which the long-term evolution of the interaction force can be expected, i.e. predicted, if the cart trajectory is known. Unlike cross-correlation, which is limited to linear relations between variables, mutual information assesses both linear and nonlinear dependency. It is therefore more suitable for this nonlinear system. In particular, mutual information has been commonly used to quantify predictability of weather and climate, which are modeled by chaotic dynamical systems (DelSole 2004; Kleeman

2011).

The cart trajectory, which was close to sinusoidal, was represented by its phase in state space φ(t) = arctan(𝑋 ´/ (2 f π X)). The interaction force Finter(t) was used as defined above. The predictability measure MI was therefore

𝑀𝐼(𝜑, 𝐹 𝑖𝑛𝑡𝑒𝑟 ) = ∬ 𝑝(𝜑, 𝐹 𝑖𝑛𝑡𝑒𝑟 )𝑙𝑛 [ 𝑝(𝜑,𝐹 𝑖𝑛𝑡𝑒𝑟 ) 𝑝(𝜑)𝑝(𝐹 𝑖𝑛𝑡𝑒𝑟 ) ] 𝑑𝜑𝑑𝐹 𝑖𝑛𝑡𝑒𝑟 ( 3 )
where p denotes the probability density functions for φ(t) and Finter(t). Mutual information is a dimensionless quantity, and its unit depends on the base of the logarithm that is used. Here, the natural logarithm was used, and the unit of mutual information is the nat.

Hypothesis 3 -Exploit Resonance: Determining the resonance structure of the system requires analytical or numerical analysis of the system dynamics and cannot be inferred from the behavioral data alone.

Therefore, Hypothesis 3 will be addressed later in the modeling and simulation section.

Data Processing: For all kinematic and result variables, only the data between t = 20 s and t = 40 s of each trial were analyzed to eliminate transients at the beginning and end of the trial. As the experimental data were compared with model simulations described below, trials that significantly deviated from periodicity needed to be excluded as the model assumed periodicity. Hence, when the standard deviation of the oscillation frequency exceeded 10% of its mean, the trial was excluded as this indicated significant deviation from the instructed periodic movements. Similarly, a trial was excluded if the mean cart excursion was smaller than 12 cm or larger than 21 cm, as it did not satisfy the instructed excursion (15 to 18 cm), even allowing an additional 3 cm of tolerance. These relatively stringent inclusion criteria were adopted in post-processing only to enable meaningful comparison with the simulation study reported below (the simulation assumed constant movement frequency within a given amplitude range).

They were not success/failure criteria for the participants. One participant's majority of trials did not satisfy these criteria and his entire data were eliminated from subsequent analysis. From the remaining 450 trials of 9 participants, only 17 trials did not meet these criteria. These 17 trials were not at the beginning of the experiment, but distributed across early and late trials. This indicated that the task did not require practice, and performing with periodicity was not a challenge per se.

The data processing and analyses were performed with MATLAB® (The Mathworks Inc., Natick, MA) and Gnumeric. The numerical values of the interaction force and predictability estimates for each experimental trial were computed with Matlab from the experimental trajectories. Mutual information was calculated with the Matlab MIToolbox-2.1.2. Statistical comparisons were performed using t-tests since the measures were normally distributed (confirmed by Kolmogorov-Smirnov tests).

Results

Task Performance and Kinematic Variables: As a first overview of participants' performance, Fig 3 shows the frequencies fk adopted by participants plotted as a histogram. To obtain a sufficiently large number of data, each cycle, i.e. one back-and-forth movement, was a data point. Two distinct strategies were observed: frequencies were concentrated either between 0.4 and 0.7 Hz (low-frequency strategy) or between 0.9 and 1.8 Hz (high-frequency strategy). The low frequencies were densely concentrated with a sharp peak at around 0.65 Hz, while the higher frequencies were distributed more broadly. These two strategies were separated by a gap between 0.7 and 0.9 Hz: only very few oscillations had a frequency within this range. Four participants adopted the low-frequency strategy, and four participants chose the high-frequency strategy. One participant used low frequencies for the first 35 trials, and then switched to high frequencies; his first 35 trials were therefore put in the low-frequency strategy, and the subsequent trials in the high-frequency strategy. All others were consistent in their choice throughout their 50 trials, excluding the very first trials that were exploration. Note that the x-axis is in log scale. Hypothesis 1 -Interaction Force: The root mean square interaction force RMSF increased from 2.57 ± 0.56 N to 5.49 ± 0.10 N in the low-frequency group, and from 5.48 ± 1.59 N to 9.09 ± 0.38 N in the highfrequency group between early and late trials. The increase was significant in both groups (p < 0.001).

This evolution suggests that participants did not minimize interaction force, counter to Hypothesis 1.

Instead, with practice they increased the exerted interaction force. Further, 5 out of the 9 participants chose the high-frequency strategy which was associated with significantly higher RMSF values. If minimization of interaction forces had been the criterion, all participants should have converged to the low-frequency strategy.

Hypothesis 2 -Predictability: Mutual information MI between the interaction force and the cart kinematics of the low-frequency group increased from 1.25 ± 0.05 nat in the first 5 trials to 1.44 ± 0.06 nat in the last 5 trials. In the high-frequency group, mutual information increased from 1.36 ± 0.08 nat to 1.53 ± 0.03 nat between early and late trials. The increase was significant in both groups (p < 0.003)

supporting Hypothesis 2 that participants sought to increase predictability of the system they interacted with. Note that though the increase in MI seemed modest, the maximum achievable value of MI was around 1.8 nat (for achievable oscillation frequencies). Therefore, the observed relative increases were important.

Simulations and Analysis of the Result Space

The results of the behavioral experiment provided support for Hypothesis 2 that humans strive to increase the predictability of the interaction when manipulating an inherently erratic or unpredictable system.

Conversely, the interaction force was not minimized in this interactive task (counter to Hypothesis 1). To further evaluate these findings and to test Hypothesis 3, we compared the strategies adopted by participants with possible alternative executions to shed light on priorities in human control. To this end, model simulations were performed to compute the result variables for alternative executions that could have achieved the task.

A Coupled Model

In a previous study, the task dynamics was analyzed by considering the behavior of the cart-andpendulum system alone without including the controlling hand [START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF]). However, this uncoupled model only partly replicated our experimental data (see Appendix A). We therefore extended the model to include the continuous coupling between the cart and the hand.

Mechanical Model and Forward Dynamics:

To capture the dynamics of the task more accurately, the cart-and-pendulum system was coupled to the hand dynamics (Fig 7). The hand dynamics was represented by an ideal force generator (force Finput) in parallel with a spring (stiffness K) and a damper (damping coefficient B). 𝐹 𝑖𝑛𝑝𝑢𝑡 (𝑡) was the force required to follow a desired trajectory (𝑋 𝑑𝑒𝑠 (𝑡), 𝑋 ´𝑑𝑒𝑠 (𝑡)). If the full dynamics of the taskincluding the pendulum forcewere perfectly anticipated, participants would be able to generate an input force Finput allowing the cart to exactly follow the desired trajectory Xdes(t). In reality, however, it was unlikely that participants learnt the perfect model due to the pendulum force acting as a perturbation. Therefore the motion due to the generated input force Finput(t) did not exactly track the desired cart trajectory, so that the actual cart trajectory X differed from Xdes. The spring and damperwhich were a simplified model of hand impedancethen served to resist this perturbation.

Note that this model represented the impedance at the level of the limb: the stiffness K and damping B corresponded to limb features and not to properties of the involved muscles. The equations of motion of the coupled model are

(𝑚 𝑐 + 𝑚 𝑝 ) 𝑋 ̈= 𝑚 𝑝 𝑑 [ 𝜃 2 ̇sin 𝜃 -𝜃 ̈cos 𝜃] + 𝐹 𝑖𝑛𝑡𝑒𝑟 = 𝐹 𝑏𝑎𝑙𝑙 + 𝐹 𝑖𝑛𝑡𝑒𝑟 𝜃 ̈= - 𝑋 d cos 𝜃 - 𝑔 𝑑 sin 𝜃 (4)
𝐹 𝑖𝑛𝑡𝑒𝑟 = 𝐹 𝑖𝑛𝑝𝑢𝑡 -𝐾(𝑋 -𝑋 𝑑𝑒𝑠 ) -𝐵(𝑋 ̇-𝑋 ̇𝑑𝑒𝑠 )

Given the task instructions, the desired trajectory was a sinusoid Xdes(t) = A sin(2 π f t + π/2). The coupled model was simulated with forward dynamics, i.e. computing the system state variables 𝑋(𝑡), 𝑋 ´(𝑡), 𝜃(𝑡), 𝜃 ´(𝑡) and interaction force Finter(t) from a known Finput(t). Since Finput(t) could not be measured experimentally, it was chosen to match the force required to manipulate a rigid object of similar mass, i.e. 𝐹 𝑖𝑛𝑝𝑢𝑡 (𝑡) = (𝑚 𝑐 + 𝑚 𝑝 ) 𝑋 ̈𝑑𝑒𝑠 (𝑡). Humans can manipulate rigid objects very accurately, suggesting that they have a good model of the task dynamics. The hand impedance parameters K and B were considered constant during a trial.

Execution Variables: To evaluate the three hypotheses, one must first define a "strategy": a strategy was defined by the set of execution variables that participants directly controlled and that fully determined the task outcome (and hence referred to as result variables). While the cart oscillation amplitude A was prescribed in the experiment, participants could freely choose three variables of the coupled model: the movement frequency f, the hand stiffness K and the damping B, referred to as execution variables.

Unlike the movement frequency f, the experimental hand stiffness and damping could not be measured directly, but had to be estimated to afford forward simulations. To this end, an optimization was conducted which aimed to estimate the values of K and B for which the simulated cart and pendulum trajectories best resembled the experimental trajectories. The optimization process and the cost criterion

C are detailed in Appendix B.

Simulation of Result Variables and Hypothesis Testing:

As for the behavioral experiment, the simulation tested the hypotheses by evaluating the result variables root mean squared interaction force RMSF (Eq. 2) and mutual information MI between the cart kinematics and the interaction force (Eq. 3).

To obtain the space of all executions spanned by execution variables f, K and B forward dynamics simulation of the coupled model were run to generate the profiles of the cup kinematics φ(t) and the interaction force Finter(t). Using Matlab-Simulink, the simulation time was 45 s, but only data from 20 ≤ t ≤ 40 s were analyzed to eliminate transients. The two result variables MI and RMSF were then calculated with Matlab as for the experimental data. These results then served to test Hypotheses 1 and 2.

To evaluate Hypothesis 3 (exploit resonance), a frequency response analysis of the coupled model was conducted in Matlab. Due to the nonlinearity of the coupled cart-and-pendulum plus human hand system, classic frequency response tools could not be used. However, the system could be linearized assuming small pendulum angles. Although this approximation was not valid for all frequencies, the linear analysis allowed further insight into the behavior of the system. In the frequency response analysis, only one of the execution variables, the movement frequency f, was varied, while the hand stiffness K and damping B were fixed to typical values: one corresponding to the mean values of K and B adopted by participants in the low-frequency group, and the other to the mean values in the high-frequency group (see Appendix B for the identification procedure of experimental values of K and B). ; therefore, the two hypotheses of interaction force minimization and predictability maximization are mutually exclusive. Conversely, for frequencies around 0.64 Hz and higher than 1.20 Hz, predictability was high, but interaction force was high as well. This fast change in MI was likely due to the resonance structure of the system detailed below. The more intricate variation of MI at higher frequencies might be due to chaotic behavior. The data suggest that participants adopted strategies with relatively high MI or high predictability.

Simulation Results of the Coupled Model

Additionally, MI was computed from the time series of the simulated data and is presented in Fig 6D . MI increased from 1.11 ± 0.05 nat in the early (first 5) trials to 1.30 ± 0.03 nat in the late (last 5) trials in the low-frequency group (p = 0.003). In the high-frequency group, the simulated MI increased from 1.21 ± 0.07 nat to 1.29 ± 0.02 nat (p = 0.02). Again, note that the maximum value of MI was about 1.8 nat. Surprisingly at first sight, the second peaks at the higher frequencies are hardly noticeable. This arises from the fact that the simulation assumed that subjects generated a sinusoidal predictive force Finput(t) intended to produce the desired cart motion Xdes(t). This predictive force was based on an incomplete model of the object dynamics which considered only its lowest-frequency mode of behavior, i.e. as though the pendulum and the cart moved as one body 𝐹 𝑖𝑛𝑝𝑢𝑡 = (𝑚 𝑐 + 𝑚 𝑝 )𝑋 ´𝑑𝑒𝑠 . This imperfect predictive force only partially compensated for object dynamics, which was nevertheless sufficient to counteract the object's resonances, especially at the higher frequencies. Mathematically, the predictive force introduced complex-valued zeros near the complex-valued poles that describe the high-frequency resonance. These zeros tended to cancel or 'mask' the effect of the adjacent poles, converting a sharp resonant peak into a broad region of nearly-constant magnitude (see Footnote 1).

Comparing this progression

Importantly, the response of cup displacement for both systems shows a sharp valley, indicating the antiresonance at 0.74 Hz between the two resonances. Note that the anti-resonance frequency is identical in system A and B, i.e. independent of the values of K and B. The phase plots in Fig 10 display the relative phase between the input force and the cart movement (red line), and the relative phase between the input force and the pendulum movement (blue line). Comparison between these two curves highlights that for low frequencies the cart and pendulum are in-phase, while for frequencies higher than the anti-resonance frequency, cart and pendulum motions are anti-phase. In addition, the relative phase between the input force and the cart movement (red line) reveals that for frequencies outside the two resonance frequencies, the cart movement is anti-phase with the input force. Conversely, over a small interval between the two resonance frequencies, the relative phase between the input force and the cart movement is changing.

For comparison of the model's resonant peaks with the experimental data, the distributions of the observed frequencies in participants are shown in grey ( Fig 10). For the low-frequency group (System A) the peak in the distribution is very close to the system's resonance peak. For the high-frequency group, participants show a very broad distribution that matches with the smeared-out resonance peak of System B. Comparison between Fig 9 and 10 reveals that the two resonance frequencies of the system coincided with areas of high MI. This suggests that the behavior of the system is easily predictable when oscillating at a resonance frequency. Conversely, the anti-resonance frequency coincides with a region of low MI, therefore the behavior of the system is hard to predict when oscillating at or around the anti-resonance frequency. These results are consistent with Hypothesis 3.

Discussion

This study examined strategies that humans adopt when manipulating objects with underactuated internal dynamics. To date, the majority of research in motor neuroscience has examined unconstrained movements in highly controlled experimental tasks to render interpretable data; only relatively few studies have examined control of complex objects. However, everyday behavior is full of complex manipulations that set humans apart from primates and other animals. The present study focused on continuous physical interaction with a cart-and-pendulum system, representing the simplified dynamics of a moving a cup of coffee. Participants had to move with a prescribed amplitude, but could choose their preferred frequency. Importantly, in continuous interaction with the complex object, the dynamics of this system is underactuated and can exhibit erratic and unpredictable behavior. Such unpredictable dynamics poses significant challenge to any internal model guiding the goal-directed manipulation.

Using both behavioral data and numerical analysis of the cart-and-pendulum system coupled to a model of hand impedance, we tested three hypotheses: humans minimize the interaction force required to move the system (Hypothesis 1); alternatively, they maximize predictability of the system behavior (Hypothesis 1); and/or they exploit the resonance structure of the system (Hypothesis 3). Interaction force between hand and cart was quantified by its root mean squared value. Predictability was operationalized by the mutual information between the kinematics of the cart and the interaction force. Exploiting resonance was tested by comparing the chosen frequencies with the resonance structure of the system. Results of the experiment showed that participants increased, not decreased, the interaction force (counter to Hypothesis 1), while they also increased predictability of the system with practice (consistent with Hypothesis 2). Half the participants chose a strategy that had significantly higher interaction forces, while affording similarly high degree of predictability.

The results of the simulations gave further support that, among alternative strategies (defined by values of movement frequency and hand impedance that humans could adopt), participants chose strategies with high predictability, but not with low interaction force. These results corroborate and generalize those obtained by [START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF] in a similar experiment that prescribed movement frequency, but left amplitude free to choose. In addition, frequency response analysis of the linearized coupled system showed that participants chose movement frequencies close to the resonance frequencies of the system, while avoiding the anti-resonance frequency (consistent with Hypothesis 3). These findings demonstrate that predictability is a control priority in complex underactuated object manipulation, which takes precedence over principles such as interaction force minimization. The fact that results support both Hypothesis 2 and Hypothesis 3 suggests that predictability may be explained by the resonance structure of the system. Therefore, manipulation of underactuated objects cannot be understood simply by extending principles of free movements or rigid object manipulation; underactuated object manipulation constitutes a different class of tasks with different control challenges.

Assumptions of the Coupled Model

To provide an entry to a quantitative understanding of this complex task, an essential element in our approach was simulation of the task dynamics with only minimal assumptions about the controller. We therefore coupled a simplified model of hand impedance to the cart-and-pendulum system. This coupled model approximated the experimental data more accurately than a previous model with the cart-andpendulum alone (Appendix A). However, as this model went beyond the physics of the task alone and included the human controller, certain assumptions had to be made.

Invariance of Input Force:

One first assumption was that the input force (Eq 4) was equal to the force required to move a rigid object of the same mass as the cart-and-pendulum system; further, the amplitude, frequency, and phase of this input force was the same sinusoidal signal during and across trials. While this is a reasonable initial assumption, it is likely that humans learned to adapt their input force, based on the perceived interaction force and/or the cart displacement. As the simulation kept the input force invariant, the desired cart trajectory was not always accurately tracked, especially when the hand impedance was low. A plausible next modeling step would be to modulate the amplitude of the sinusoidal input force based on the difference between the actual and desired cart amplitude. Even though it is relatively straightforward to include such an adaptation of the input force, this would evidently make the model more complex and not necessarily help to understand the data.

Invariance of Hand Impedance:

A second simplifying assumption was that the hand impedance was constant throughout one trial. Given the task instruction and the virtual display, the amplitude of the cart movement was the main concern for participants, while the actual trajectory between the two targets was secondary. Therefore, it could be speculated that participants may increase their arm impedance close to the targets to ensure accuracy in the amplitude, but decrease impedance during translation between targets. A sinusoidally changing impedance might therefore better match experimental data. However, as with the modulation of input force, the potential gain in realism would be at the cost of more parameters to identify. Therefore, constant impedance and constant input force is a reasonable compromise between accurate replication of experimental data and transparency of the model.

Predictability, Muscular Effort and Antagonist Co-Contraction

The simulations reveal that high predictability and low interaction force are non-overlapping strategies and the data provide evidence that it is predictability that determines the choice of control strategy. The finding that humans do not try to minimize interaction force may seem to run counter to many studies on unconstrained movements that have shown that humans favor energy-or effort-efficient strategies [START_REF] Nelson | Physical principles for economies of skilled movements[END_REF][START_REF] Alexander | Energy-minimizing choices of muscles and patterns of movement[END_REF][START_REF] Prilutsky | Optimization-based models of muscle coordination[END_REF]. It should be pointed out that our force criterion only quantified the net external force, i.e. interaction force. While this external force increased, it might be that higher predictability had a secondary effect on decreasing internal muscular effort: when the system dynamics is erratic, it is difficult to anticipate and preempt the perturbing force of the pendulum by feedforward control. The user may then rely on his/her hand impedance to reject these perturbations and maintain the desired cart trajectory. This requires increasing the impedance through co-activation of antagonist muscles, which results in higher muscular effort without any consequences on the net external force. Conversely, predictable object dynamics may enable participants to anticipate the perturbing interaction force, and thereby reduce effort due to co-contraction. Predictability can therefore afford a way to minimize the overall muscular effort.

The strongest evidence that force minimization was not an objective was that half of the participants chose the high-frequency strategy associated with higher forces than the low-frequency strategy (Fig 6).

If effort were the main concern, all participants should have chosen the lower frequency and lower impedance (Appendix B). As mutual information was similar in both frequency groups, the lowfrequency solution would have decreased the overall effort and reconciled the predictability and interaction force objectives. However, one point to note is that the task required only relatively low forces, which may be one reason why optimizing effort was not a priority. Testing the same experiment with different masses for the cart-and-pendulum system is a direction for future work.

Predictability, Error Correction and Computational Cost

Another factor that may have influenced participants' choices was that the low-frequency strategy was close to the boundary of the low predictability zone (starting around 0.7 Hz in Fig 9 ), compared to the high-frequency solution that was more robust or tolerant to variation in frequency. With the lowfrequency strategy, small variations could easily lead to erratic behavior and perturbations that require correction. If such error corrections were executed by the CNS, then the computational cost would increase. Computational effort has been recognized and included as a cost in several optimization studies (Todorov and Predictability of the interactive dynamics of complex object manipulation may again be a manifestation of human controllers seeking to simplify the control task.

Resonance/Anti-Resonance Structure, Effort and Predictability

Did participants choose to move at resonance peaks to reduce effort? As Fig 10A showed, participants who moved the cart and pendulum in phase could take advantage of the low-frequency resonance to reduce effort, but had to exert precise control of frequency to avoid the nearby anti-resonance frequency.

Participants who chose the anti-phase strategy expended more muscular effort due to the higher frequency of anti-phase motion and to the elevated stiffness and damping they exhibited. However, the anti-phase motion was available over a much broader range of frequencies (Fig 10B ) and therefore required much less precise control of frequency. Further, they were far away from the anti-resonance frequency or dynamic zero at 0.74 Hz.

Did participants prefer certain cup frequencies because they were associated with specific relative phases between the cart and the pendulum movements or between the input force and the cart movement?

Several studies on rhythmic bimanual coordination have shown that humans prefer in-phase and antiphase relations between two limbs over other phase relations [START_REF] Kelso | Phase transitions and critical behavior in human bimanual coordination[END_REF][START_REF] Schöner | A synthetic theory of environmentally-specified and learned patterns of movement coordination[END_REF][START_REF] Sternad | Average phase difference theory and 1:1 phase entrainment in interlimb coordination[END_REF][START_REF] Sternad | synaptic and synergetic coupling: An evaluation through in-phase and antiphase rhythmic movements[END_REF]. In the present experiment, participants also oscillated the cart either in-phase and movement is what one gets in the very common situation of manipulating a rigid object. However, if relative phase was the only concern, participants' data points would be spread over all the frequencies with anti-phase coupling, and not grouped over a narrow frequency range. These observations support that potential phase preferences alone do not account for our observations.

Why did participants avoid the anti-resonance frequency? At anti-resonance, the force generated by the pendulum movement (Fball in Eq 4) exactly opposes the interaction force exerted by the human (Finter in Eq 4), resulting in zero displacement of the cart. In addition, near the anti-resonance frequency the relation between cart motion and input force undergoes a large and rapid, almost discontinuous, phase shift, whereas the relation between pendulum motion and input force does not (phase plot in Fig 10).

Around the anti-resonance frequency, the oscillations of the cart and pendulum desynchronize very quickly and small variations result in large changes in the direction of the perturbing force due to pendulum motion. This makes the compensatory input force that should be applied to obtain the desired cart movement hard or impossible to predict. The results clearly showed that subjects consistently avoided the anti-resonance frequency and, implicitly, favored predictability.

A Task-Dynamic Approach, Internal Models and Predictability

Most computational studies on movement control start with a hypothesis about the human controller. For example, several studies of the pole-balancing task investigated specific hypotheses about the neural control system, ranging from different control models to the role of noise or sensory feedback (Mehta In contrast, our task-dynamic approach shifted the emphasis to first understand the task and its affordance, while minimizing assumptions about human neuromotor control [START_REF] Sternad | Control of intermittent and continuous objects[END_REF]. Starting with a mathematical model of the task and analysis of its dynamics, the solution space can be derived and human solutions can be evaluated. To make this mathematical approach transparent a simplified model is advantageous. Here, we reduced the fluid dynamics of the coffee to a single degree of freedom. As with any virtual implementation, this may raise the question whether the problem has become too simple and results will generalize to the real cup of coffee. Recently, two theoretical studies have indeed analyzed the cup of coffee system in its full physical complexity (Mayer & Kretchetnikov 2012[START_REF] Han | A Study on the Coffee Spilling Phenomena in the Low Impulse Regime[END_REF]). Comparison of these and our studies may reveal the advantages and disadvantages of the realistic versus computationally simplified approach.

Our task-based approach does not contradict, but complement controller-based approaches. When for example [START_REF] Nagengast | Optimal control predicts human performance on objects with internal degrees of freedom[END_REF] studied optimal control for the manipulation of a virtual mass-springdamper system, they assumed that participants had complete knowledge of the system dynamics.

Similarly, [START_REF] Dingwell | Manipulating objects with internal degrees of freedom: evidence for model-based control[END_REF][START_REF] Dingwell | Experimentally confirmed mathematical model for human control of a non-rigid object[END_REF] showed that participants manipulating a linear mass-spring system displayed behavior compatible with learning an internal model of the object dynamics. However, underactuated objects, such as our cup-and-pendulum system pose a significant challenge due to their possibly unpredictable dynamics leading to an apparent absence of correlation between the human action and the resulting behavior of the system. Increasing the predictability of object dynamics might therefore be a way to increase the chance of acquiring an internal model.

Footnotes

Appendix A: Limitations of a Model without Hand Impedance

In a previous study, the dynamics of the cup-and-ball task was analyzed by looking at the behavior of the cart-and-pendulum system alone without the controlling hand [START_REF] Nasseroleslami | Rhythmic manipulation of objects with complex dynamics: predictability over chaos[END_REF]). This uncoupled model is depicted in Fig A1 and the motion of the system is described solely by Eq 1. It is straightforward to simulate this uncoupled model using inverse dynamics calculations: if the cart trajectory X(t) and initial conditions of the cart and pendulum (𝑋 0 , 𝑋 ´0, 𝜃 0 , 𝜃 ´0) are given, the pendulum trajectory θ(t) and the interaction forces Finter(t) can be computed using Eq 1 and a numerical integration scheme for θ. This uncoupled model has the advantage that it does not require any assumptions about control by the human (contrary to the coupled model). The only assumption is about the movement of the cart, which could reasonably be modeled by a sinusoid X(t) = A sin (2 π f t + π/2) given the task instructions. To quantify the divergence, the root-mean-square errors (RMS) between the experimental and simulated trajectories were computed. Table A1 summarizes RMS error for each quantity X, Ẋ, θ, θ̇ and Finter, expressed as percent of its respective maximum value in the corresponding experimental trial. In the high-frequency group, the RMS error was small and fairly consistent across variables (median RMS error around 10% of the variable maximum experimental value), indicating a reasonably good match between the experimental and simulated profiles. This uncoupled model was therefore a competent representation of the cup-and-ball task for the high-frequency strategy. With the low-frequency strategy, however, the RMS error varied greatly and reached up to 30 % of the maximum value for the experimental pendulum angle and angular velocity (and interaction force to a lesser extent). These discrepancies between experimental and simulated data demonstrate that the uncoupled model did not represent the execution strategies adopted by the low-frequency group sufficiently accurately. A likely reason for the divergence between experimental and simulated data is the assumption of a perfectly sinusoidal cart trajectory in the simulations, whereas experimental trajectories exhibited small deviations from this ideal shape. Given the sensitivity of the cart-and-pendulum dynamics to initial conditions, small changes in the participant's movement could lead to significant changes in the system evolution. These deviations of the experimental cart trajectories from a perfect sinusoid could have two main causes: the intrinsic variability of human movements, and the perturbations caused by the internal dynamics of the object. The first cause results from the ever-present human variability: even if the object was rigid, or if there were no object at all, humans are unable to repeat the same exact movements. While present in both frequency strategies, this variability could have different consequences, since the sensitivity of the system to initial conditions is not constant. trials were evidently not captured in the simulation as constant desired cart amplitude/frequency was assumed. The simulations used the average experimental values of A and f across all cycles of the trial (20 ≤ t ≤ 40 s) to define the desired trajectory Xdes(t) = A sin(2 π f t + π/2) and the input force Finput(t) = (mc + mp) Ẍdes(t). However, the average amplitude and frequency were only representative of the experimental trial if they did not vary significantly throughout the trial. This motivated the stringent inclusion criteria in the analysis of the behavioral data.

where the superscripts s and e stand for simulation and experimental, respectively. Only the data within 20 ≤ t ≤ 40 s were included to avoid confounding by transients (both for experimental and simulated trials).

While the movement frequency f was fixed in the simulations, experimental frequencies were not exactly constant within trials. Such variations of the experimental frequency created a temporal offset between the experimental and simulated trajectories, which could lead to high RMS errors even when the two profiles were similar. To limit this artifact, C was computed cycle by cycle, i.e. the RMS errors were computed for each cycle k by time-aligning the experimental and simulated trajectories of cycle k.

Subsequently, they were averaged over all cycles.

Across all trials, the median cost C measured for the best impedance fit of each trial was 0.104 with an interquartile range of 0.051. Table B1 gives the ratio between the RMS error between experimental and simulated time-series and the maximum experimental value of the corresponding trial for the state variables (𝑋, 𝑋 ´, 𝜃, 𝜃 ´) as well as for the interaction force Finter. The median value of the RMS error was between 9 and 13% of the maximum value, depending on the variable. Importantly, the error was consistently low in both groups, unlike for the uncoupled model above (see Table A1 in Appendix A). The values of hand impedance were different between groups. The comparison of stiffness and damping values between the two frequency groups was performed with a Wilcoxon signed rank test because the data were not normally distributed. Both the stiffness K and damping B were significantly lower in the low-frequency group, with p = 10e-10 and p = 10e-13 respectively. This is consistent with the known fact that, for a similar task accuracy, limb stiffness usually increases with movement speed.

These results are the basis for characterizing experimental trials with hand impedance. The coupled model with optimized K and B reproduced experimental trajectory and force time-series much more accurately than the uncoupled model (especially for the low-frequency group), thus confirming its better competence to analyze the experimental task.

Fig 1 .

 1 Fig 1. Model of the task. A: Conceptual model of the cup-and-ball system. B: Mechanical model of cup-and-ball

Fig 2 .

 2 Fig 2. Experimental set-up of the ball-and-cup task using virtual reality and force feedback. A: Rendering

Fig 3 .

 3 Fig 3. Distribution of frequencies adopted by all participants when manipulating the virtual cup-and-ball

Fig 4

 4 Fig 4 depicts a low-and a high-frequency strategy with exemplary time series of the cart and pendulum

Fig 4 .

 4 Fig 4. Experimental cart and pendulum trajectories. Representative trajectories of the cart (top panel) and

Fig 5

 5 Fig 5 shows how the kinematic variables A, f and the relative phase between the cart and pendulum

Fig 5 .

 5 Fig 5. Evolution across trials of the experimental kinematic variables. A: Amplitude A of the cart oscillations.

Fig 6 .

 6 Fig 6. Evolution across trials of the result variables. Evolution of the experimental (A, C) and simulated (B, D)

Fig 7 .

 7 Fig 7. Model used to analyze the dynamics of the task in simulation. Forward dynamics of the cart-and-

Figs

  Figs 8A and 9A display the 3D execution space spanned by frequency f, stiffness K and damping B. For

Fig 8 .

 8 Fig 8. 3D plot and 2D contour map of RMSF in the space of the execution variables. A: 3D plot of the root

Fig 9 .Hypothesis 1 -Hypothesis 2 -

 912 Fig 9. 3D plots and 2D contour map of MI in the space of the execution variables. A: 3D plot of the mutual

  with the experimental values (Fig 6C) shows that both the time course and the magnitudes of the MI simulated values were close to the experimental values, supporting the adequacy of the coupled model and the estimated values of stiffness and damping. This simulation result strengthens the experimental results that predictability was increased with practice. Hypothesis 3 -Resonance: One essential feature of the task dynamics is its resonance structure: the coupled system has two resonance peaks and one anti-resonance frequency or dynamic zero between the two resonance frequencies. Fig 10 displays Bode magnitude and phase plots of the linearized coupled model for two representative values of hand impedance. System A was simulated with K = 100 N/m and B = 10 N.s/m, values that were typical for the low-frequency group. System B with K = 200 N/m and B = 15 N.s/m was typical for the high-frequency group. As the responses of the two systems reveal, the resonance peaks depend on the values of K and B. The panels for pendulum angle show one clear resonant peak at 0.68 Hz for system A and at 0.71 Hz for system B.

Fig 10 :

 10 Fig 10: Bode amplitude and phase plots of the linearized coupled model, for different values of hand

(

  at low frequencies) or anti-phase (at high frequencies) with the ball movements and avoided intermediate relative phases at the anti-resonance frequency. However, this observation does not imply that participants chose strategies for their relative phase values. Except at anti-resonance, the task dynamics did not allow other relative phases as the frequency response plots show (Fig 10). The entire frequency range below 0.65 Hz corresponds to in-phase coupling, but participants of the low-frequency group nevertheless all converged to a narrow area of high predictability (Fig9). Similarly, the highfrequency group favored those subsets of the frequency range with high predictability. In addition, a large set of frequencies outside of the two resonance frequencies correspond to anti-phase coupling between the input force and the cart movement (red line inFig 10). It is reasonable to think that participants may prefer this anti-phase coupling between what they predict (input force) and what they actually obtain (cart movement) over any other relative phase. Indeed, anti-phase coupling between force

and

  Schaal 2002; Venkadesan et al. 2007; Milton 2011; Milton et al. 2013; Gawthrop et al. 2013; Insperger et al. 2013).

Footnote 1 :

 1 With K = 100 N/m and B = 10 N.s/m, the high-frequency poles are -1.87 +/-6.72i and the zeros are -1.67 +/-5.53i (in rad/s). With K = 200 N/m and B = 15 N.s/m, the high-frequency poles are -3.06 +/-8.95i and the zeros are -2.50 +/-7.77i (in Hz).

Fig A1 .

 A1 Fig A1. Model of the dynamics of the task. Inverse dynamics model of the cart-and-pendulum system alone.

Fig

  Fig A2 displays one representative example of cart and pendulum trajectories X(t) and θ(t) and the

  Jordan 2002;[START_REF] Ronsse | Optimal control of a hybrid rhythmic-discrete task: the bouncing ball revisited[END_REF]). Yet in these modeling approaches, computational cost terms have remained unspecified placeholders for unaccounted factors contributing to human control choices. A series of studies by Sternad and colleagues have argued that the human controller may exploit the stability properties of a task to avoid computationally expensive corrections[START_REF] Sternad | Control of intermittent and continuous objects[END_REF]). Using the task of rhythmically bouncing a ball with a paddle, several experiments provided robust evidence that human subjects learned to attain dynamic stability, such that small errors passively decayed, obviating the need for explicit corrections[START_REF] Schaal | One-handed juggling: A dynamical approach to a rhythmic task[END_REF][START_REF] Sternad | Dynamics of a bouncing ball in human performance[END_REF][START_REF] De Rugy | Actively tracking 'passive' stability in a ball bouncing task[END_REF]. When applying larger perturbations, additional corrections were evidenced, although the signature of dynamic stability was still visible[START_REF] Siegler | Passive vs. active control of rhythmic ball bouncing: the role of visual information[END_REF][START_REF] Wei | Passive stability and active control in a rhythmic task[END_REF][START_REF] Wei | Stability and variability: Indicators for passive stability and active control in a rhythmic task[END_REF]. In a similar spirit, mathematical and empirical studies of a throwing task showed that humans seek solutions that are tolerant to error and noise, therefore requiring fewer corrections[START_REF] Sternad | Neuromotor noise, sensitivity to error and signal;dependent noise in trial-to-trial learning[END_REF][START_REF] Sternad | Acquisition of novel and complex motor skills: stable solutions where intrinsic noise matters less[END_REF][START_REF] Cohen | Variability in motor learning: relocating, channeling and reducing noise[END_REF]).

Table A1 . RMS error between experimental and simulated trajectories and force time-series for the uncoupled model.

 A1 Ratio of RMS error between experimental and simulated data normalized by the maximum value for the cart and pendulum trajectories and interaction force in both subject groups. The simulated data were obtained from inverse dynamics simulation of the uncoupled model. The median and interquartile range were computed over all 433 valid trials.

			Low-frequency group High-frequency group
			Median	IQR	Median	IQR
	𝑟𝑚𝑠(𝑋 𝑒 -𝑋 𝑠 ) ‖𝑋 𝑒 ‖ ∞		0.10	0.04	0.08	0.02
	𝑟𝑚𝑠(𝑋 ´𝑒 -𝑋 ´𝑠) ‖𝑋 ´𝑒‖ ∞		0.13	0.06	0.08	0.03
	𝑟𝑚𝑠(𝜃 𝑒 -𝜃 𝑠 ) ‖𝜃 𝑒 ‖ ∞		0.29	0.52	0.13	0.09
	𝑟𝑚𝑠(𝜃 ´𝑒 -𝜃 ´𝑠) ‖𝜃 ´𝑒‖ ∞		0.31	0.39	0.11	0.07
	𝑟𝑚𝑠(𝐹 𝑖𝑛𝑡𝑒𝑟 𝑒 ‖𝐹 𝑖𝑛𝑡𝑒𝑟 -𝐹 𝑖𝑛𝑡𝑒𝑟 𝑠 𝑒 ‖ ∞	)	0.22	0.29	0.12	0.04

  All combinations of 10 ≤ K ≤ 350 N/m (step size 2 N/m) and 3 ≤ B ≤ 50 N.s/m (step size 1 N.s/m) were tested to find the best fit. The difference between the experimental and simulated trajectories was quantified by the cost C of the normalized root mean square errors of the four quantities X(t), Ẋ(t), θ(t),

	θ(t)
	𝐶 =	1 4	[ 𝑟𝑚𝑠(𝑋 𝑒 -𝑋 𝑠 ) ‖𝑋 𝑒 ‖ ∞	+	𝑟𝑚𝑠(𝑋 ´𝑒 -𝑋 ´𝑠) ‖𝑋 ´𝑒‖ ∞	+	𝑟𝑚𝑠(𝜃 𝑒 -𝜃 𝑠 ) ‖𝜃 𝑒 ‖ ∞	+	𝑟𝑚𝑠(𝜃 ´𝑒 -𝜃 ´𝑠) ‖𝜃 ´𝑒‖ ∞ ]

Table B1 : RMS error between experimental and simulated trajectories and force time-series for the coupled model.

 B1 Ratio between root mean square error RMS between experimental and simulated data and the maximum value for the cart and pendulum trajectories and interaction force. The results are separated for the two frequency groups. The simulated data were obtained with forward simulation of the coupled model, using the optimized values of K and B for each trial (i.e. the values for which the cost C was minimum).

			Low-frequency group	High-frequency group
			Median	IQR	Median	IQR
	𝑟𝑚𝑠(𝑋 𝑒 -𝑋 𝑠 ) ‖𝑋 𝑒 ‖ ∞		0.09	0.04	0.09	0.02
	𝑟𝑚𝑠(𝑋 ´𝑒 -𝑋 ´𝑠) ‖𝑋 ´𝑒‖ ∞		0.11	0.05	0.08	0.03
	𝑟𝑚𝑠(𝜃 𝑒 -𝜃 𝑠 ) ‖𝜃 𝑒 ‖ ∞		0.11	0.07	0.12	0.07
	𝑟𝑚𝑠(𝜃 ´𝑒 -𝜃 ´𝑠) ‖𝜃 ´𝑒‖ ∞		0.12	0.07	0.10	0.05
	𝑟𝑚𝑠(𝐹 𝑖𝑛𝑡𝑒𝑟 𝑒 ‖𝐹 𝑖𝑛𝑡𝑒𝑟 -𝐹 𝑖𝑛𝑡𝑒𝑟 𝑠 𝑒 ‖ ∞	)	0.13	0.05	0.13	0.04

Acknowledgments

The authors wish to thank Fei Ye for her help with the data collection and first analysis.

Grants

Dagmar Sternad was supported by the National Institutes of Health R01-HD087089, R01-HD081346, and R21-DC013095 and the National Science Foundation NSF-NRI 1637854 and NSF-EAGER-1548514. Neville Hogan was supported by NIH R01-HD087089, NSF-NRI 1637814, NSF-EAGER-1548501 and by the Eric P. and Evelyn E. Newman fund. Pauline Maurice was supported in part by the European Union's Horizon 2020 Research and Innovation Program under Grant Agreement No. 731540.

The second causethe perturbation forces created by the pendulum movementsaffected the cart trajectory because the human hand is not an ideal position generator. Unexpected pendulum forces disrupted hand and hence cart movement. Though this is again true for both frequency strategies, the cart trajectory was likely less perturbed in the high-frequency strategy, because hand movements were faster, which is often associated with a higher hand impedance; higher impedance would result in better resistance to external perturbations and lower RMS error (Table A1).

Furthermore, the interaction force Finter results from two different forces (Eq 1): one is the cart-andpendulum inertial force Finertia = (mc + mp) Ẍ, and the other is the pendulum force Fball. The average ratio between the RMS pendulum force and the RMS inertial force (computed for 20 ≤ t ≤ 40 s) was 0.70 ± 0.16 in the low-frequency group and 0.32 ± 0.05 in the high-frequency group (averaged across all trials of all participants in each of the two groups). Relative to the expected force (i.e. required to accelerate the total system inertia, similar to the manipulation of a rigid object), the magnitude of the unexpected perturbation (the pendulum force) was thus much higher in the low-frequency group and was therefore less likely to be resisted. Hence, the current study included the effect of hand impedance on the dynamics of the cart-and-pendulum system.

Appendix B: Estimation of Hand Impedance in the Coupled Model

Unlike the movement frequency f, the experimental hand stiffness K and damping B could not be measured directly, but had to be estimated from the human data. To this end, an optimization was conducted which aimed at finding the values of K and B for which the simulated cart and pendulum trajectories most resembled the experimental trajectories. For each combination of K and B a 45 s forward dynamics simulation of the coupled model was performed and compared with the corresponding experimental trial. The continuous variations in the cart amplitude and/or frequency in the experimental