1 Introduction

When observing people carrying out actions, one may notice that they move at different pace, with some individuals being systematically slow or systematically fast. Recent research has confirmed this observation [1], identifying vigor, an idiosyncratic trait across actions that characterizes the tradeoff between the energy expenditure, accuracy, and time to perform a given goal-oriented movement [2–4], reflecting the value individuals associate with it [5]. At the neural level, vigor may be influenced by reward-driven dopamine secretion in the basal ganglia [6–10], where individuals more sensitive to reward tend to move faster. Furthermore, the relative vigor of an individual within a population is remarkably stable across repeated measurements compared to its variability between people [11]. This stability of vigor allows to predict self-paced human movements across tasks, with various environmental dynamics, as a tradeoff between minimizing a cost of time and the effort required to complete an action with the desired accuracy [2, 3, 12, 13]. Notably, this subjective cost of time can be identified using inverse optimal control techniques from the standard relationship observed between motion amplitude and duration [14].

When people with different levels of vigor perform a task together, how do they combine their motor plan and collaborate efficiently? A pioneering study [15] found that dyads of connected partners performing point-to-point arm movements tend to move faster than each individual alone. Nevertheless, this result may have been influenced by the addition of an explicit reward associated with decreasing movement time given to dyads. Current knowledge suggests that individuals can identify others’ vigor through visual observation [1], but it is unclear whether this is possible when mechanically connected with a partner through an object, such as when moving a table together. In fact, it has been suggested that during point-to-point movements with a prescribed duration, mechanically coupled individuals may move independently as they do not have sufficiently rich information to identify the partner’s motion plan [16]. Importantly, even if they could coordinate their motor plans when removing the time prescription, the remaining question is whether they would do this in a systematic way, with some dyads moving consistently faster than others across conditions. This paper investigates how dyads of connected partners perform point-to-point movements without time constraints, how partners coordinate their motor plans, and how their individual vigor shapes the joint movement.

Several mechanisms may drive the movement coordination, or absence thereof, between mechanically coupled individuals. First, each individual may execute their motor plan independently, causing large interaction efforts and quite arbitrary movement strategies for the dyad, a scenario we refer to as the co-activity (CA) hypothesis [16]. However, one could also expect that, when time is not prescribed, a leader naturally emerges and imposes their preferred motor plan to the dyad. Specifically, a commonly suggested interaction strategy [17] is that the slower (respectively, faster) individual in a dyad could infer and align with the faster (slower) partner’s movement strategy [15, 18, 19], which we refer to as leader-follower hypothesis (LF). More generally, recent studies have shown that connected individuals can use interaction forces to exchange their motion plans and both improve performance in a tracking task [20, 21]. These results suggest that connected individuals may co-adapt to establish a common motor plan using haptic communication. A negotiated motor plan could then emerge, where both partners contribute to setting the dyad’s behavior, based on a weighted average of their nominal vigor in the task. This would generalize the LF hypothesis in that it predicts that both partners similarly contribute to the dyad’s strategy by blending their original intents, a hypothesis we refer to as weighted adaptation (WA). Alternatively, other rules may underlie the dyad’s behavior. Conceivably, uncertainty about a partner’s motor plan could change an individual’s own motor plan and vigor, leading to an interactive adaptation (IA) hypothesis. According to this view, the dyad’s motor plan would not originate from a weighted averaging of the partners’ initial strategies, but from a dynamic change of their strategies in response to the current interaction dynamics, which can be represented in cost functions. According to this view, lasting after effects could be expected. These four hypotheses, schematized in Fig. 1A, yield distinct behavioral predictions that will be tested in the present study.

Figure 1:

2 Results

To understand how dyads’ partners coordinated their motor plans and test the above hypotheses, we investigated point-to-point movements without timing constraint or explicit reward. Relatively large targets were used in order to relax accuracy constraints. The partners carried out these movements while being connected through a virtual elastic band implemented using individual wrist exoskeletons, as illustrated in Fig. 1B. This setup allowed us to systematically investigate the impact of the interaction on movement planning, by varying movement amplitudes and connections’ stiffness, as well as the effort associated with movements using a viscous load. 20 participants (10 dyads) were requested to perform self-paced flexion and extension movements of the right wrist, with amplitudes between 18^° and 90^°, by controlling a cursor on their individual monitors from A₀ towards one of the targets {A₁, A₂, A₃, A₄, A₅} and back to A₀. The individual cursor position was an affine function of the corresponding wrist angle. The two partners of a dyad first performed movements in a solo session so that we could identify their individual motor strategies, in particular their vigor, prior to being connected in a dyadic session with an elastic band of stiffness κ = 0.5 Nm/rad (KL) or κ = 1.6 Nm/rad (KH). To enhance readability, the figures show results obtained with the low-stiffness connection KL (except when the two stiffness and viscous loads led to different behaviors), while complementary figures with large connection stiffness KH and other parameters referred to in the text are provided in the Supplementary materials.

2.1 Solo session

We first analyzed participants’ behavior in the solo session to verify that vigor, as defined in previous studies [2, 11, 22], is an idiosyncratic trait within our sample, and to identify their motor plan as a baseline for comparison and modeling in the dyadic session. As illustrated in Fig. 2A, the average trajectories in the passive condition P1 (with motor off) exhibit a bell-shaped velocity profile, where movement duration increases with amplitude. Submovements near the target appear for the farther targets {A₄, A₅}, as was also observed in self-paced movements of large amplitude carried out with another exoskeleton [3].

Figure 2:

Since vigor is expected to be a stable individual trait across effort levels [11, 22], we first verified that VL and VH imposed sufficient resistance to influence movement duration (ANOVA: F_2,38 = 23.2, p < 10⁻³, η² = 0.19, Supplementary Fig. S.1A). Participants moved significantly slower in VL and VH compared to P1, and slower in VH than in VL (p < 1.3·10⁻³; Cohen’s D > 0.56 for all of VL-P1, VH-P1, VL-VH). These findings confirm that the selected viscous loads required sufficient effort to impact movement duration and analyze the robustness of inter-individual differences.

Next, as shown in Fig. 2B, movement time increased linearly with amplitude for each participant (r² = 0.92 ± 0.09) and for the average participant in our sample (r² = 0.99). We leveraged the extracted affine amplitude-duration relationships to compute a vigor score v_i for each participant relative to the population [3, 22, 23]:

where is the average movement duration of the population for the k^th movement amplitude in the relevant passive condition (P1 for the solo session, KL or KH for the dyadic session), and T_i(A_k) is the average movement duration of participant i with amplitude A_k in the considered condition. The vigor score is positive, and a participant moving at the same speed as the population average has a unit vigor.

A decomposition of variance revealed that inter-individual variability of the vigor scores through the P1, VL and VH conditions accounted for 81.2% of the total variance, despite the possible effects of the viscous load on the total variance. This result is consistent with prior experiments testing the timing stability over multiple sessions [11], rather than by varying the effort in the present study. The vigor scores for P1, VL, VH exhibit distributions similar to previous studies [3, 22] (Supplementary Fig. S.1B) and correlate across conditions (Fig. 2C; p < 10⁻³; r > 0.71, Pearson tests for P1-VL, P1-VH, VL-VH), confirming the robustness of the individual vigor score across effort levels.

In summary, the solo session data — characterized by affine amplitude-duration relationships, greater inter- than intra-individual variability, and stable vigor scores across conditions with varied effort levels — confirm that vigor is an idiosyncratic trait in our sample of participants, consistent with prior findings [11, 22].

2.2 Dyadic session

Coordination of motor plans within dyads

To analyze coordination, we identified the slow and fast partner in each dyad according to their vigor in P1, thereby forming the slow partners group and fast partners group. First, we analyzed whether reaction times in these two groups were different in P1, as it could naturally lead to the emergence of a leader- follower behavior, the leader likely being the one starting earlier. However, we found no significant difference (p > 0.45, reaction time about 225 ms on average), thus the reaction time did not depend on the group in our task. Furthermore, no clear trend of adaptation was observed, suggesting that participants were able to rapidly set a compromise whilst moving towards the initial target (see Supplementary Fig. S.3).

To test the CA hypothesis, we modeled the fast and slow partners as planning trajectories independently with their respective preferred movement duration. Connecting these planned trajectories via an elastic band of stiffness KL or KH, we computed the CA trajectories and interaction forces by assuming participants optimally tracked their desired trajectory with minimum effort and error (see Methods). The averaged partners’ trajectory (solid black line in Fig. 3A) exhibits a non-smooth slope change, due to the submovement and stabilization of the fast partner’s within the target (solid green curve in Fig. 3A), and produced a relatively large interaction torque that increased with vigor difference between the partners.

Figure 3:

The data of connected participants in our experiment did not correspond to these predictions. First, averaged participants’ trajectories in the connected session exhibited bell-shaped velocity profiles without any visible non-smoothness (Fig. 3B), which was further confirmed as participants’ trajectories were even smoother when connected than in solo performances (analyzed with spectral smoothness metrics [24], see Supplementary Fig. S.6). Second, the absolute interaction torque was low across all conditions, irrespective of target, viscosity, or connection stiffness (Fig. 3C, Supplementary Figs. S.5C & S.7–S.16). It was increasing slightly with movement amplitude (main effects of the condition and target (F_4,36 > 18.4, p < 1.3·10⁻³, η² > 0.24 for KL and KH, confirmed by pairwise comparisons, p < 0.03, Cohen’s D > 0.83 for all conditions and targets) but was always less than 0.1 Nm (p < 10⁻³, Cohen’s D > 3.73 in the six connected conditions), and at maximum close to 0.05 Nm for target A₅ with KL. Interestingly, interaction torques were higher for KL than KH, likely due to deteriorated haptic communication with a more compliant connection [25]. Furthermore, interaction torque showed no significant correlation with differences in individual vigor within dyads (Fig. 3D, Supplementary Fig. S.5D). In contrast to these results, CA predicted average interaction torques for A₅ of 0.1 Nm with the KL connection (i.e. around twice higher than observed) and 0.23 Nm with KH (i.e. around ten times higher than observed). These findings rule out the CA hypothesis. They rather suggest that dyadic partners closely coordinated their motor commands, and that at least one of them adapted their individual motor plan.

Characteristics of individual movements within connected dyads

As shown in Fig. 3E (and Supplementary Fig. S.5B), the movement duration of connected participants increased linearly with amplitude, both with low and high connection stiffness, for each participant (KL: r² = 0.85 0.21; KH: r² = 0.92 ± 0.09) as well as for the average participant (KL: r² = 0.98; KH: r² = 0.99). Furthermore, the individual movement duration of connected participants was affected by VL and VH loads (F_2,38 > 85.7, p < 10⁻³, η² > 0.52 for both KL and KH conditions, as illustrated by Fig. 3F and Supplementary Fig. S.5E). As in solo trials, dyads moved slower when ν > 0 (p < 10⁻³, Cohen’s D > 1.72 in all cases) and with VH compared to VL (p < 10⁻³, Cohen’s D > 1.03 for both KL and KH). Despite potential difficulties in coordinating with the partner, movement duration was not affected by connection stiffness, independently of the required level of effort. Moreover, movement duration was not significantly different between dyads’ partners (Fig. 3G and Supplementary Fig. S.4, no main effect, p = 0.86). Altogether, these results indicate that the viscous loads were sufficient to impact the dyads’ behavior, and dyads were able to plan and coordinate their movements efficiently in all conditions.

A striking observation is that, when using KL and KH without viscous load, connected participants from both the slow and fast groups moved faster than in the solo session (F_2,36 = 42.2, p < 10⁻³, η² = 0.24), with a greater change in movement duration for the slow partners (F_1,19 = 11.3, p < 10⁻³, η² = 0.04, see Fig. 3H, Supplementary Fig. S.5F). Post-hoc comparisons confirmed this effect (p < 10⁻³, Cohen’s D > 0.4 for both groups). With the resistive load, the slow partners still saved time compared to solo trials (Fig. 3H, Supplementary Fig. S.5G–I, F_2,18 = 6.72, p = 6.6·10⁻³, η² = 0.11), whereas the fast partners maintained a similar duration (Supplementary Fig. S.5H,I). Post-hoc comparisons confirmed that the slow partners saved time for all the connected conditions (p < 0.04, Cohen’s D > 0.57 in all cases). Importantly, regardless of whether dyads saved time or not, it did not affect their movement smoothness or accuracy (see Supplementary Fig. S.6 and associated analyses).

In summary, dyads moved faster than lone individuals when there was no viscous load. This increased movement speed was a strong marker of unloaded dyadic movements, with, for instance, the slow partners moving 30% faster for movements to target A₅ and the fast partners 20% faster (medians in Fig. 3H). With viscous load, the slow group saved time while the fast group kept their movement duration constant (medians in Fig. 3H). The combination of all these results on interaction efforts and time supports the existence of a shared motor plan in dyad’s members, which plan is sometimes clearly differing from the lone motor plan of both participants.

Dyadic vigor

As connected partners coordinated their movements and exhibited affine amplitude-duration relationships, we computed vigor scores based on the average movement duration of the two partners of each dyad, using Eq. 1 and the same approach as for individuals.

As for the solo session, we analyzed the proportion of variance in vigor explained by inter-dyad differences throughout the connected session. When grouping all the connected conditions, we found that intra-dyad (that is variability of vigor of a dyad across blocks) variability explained 63% of the total variance. However, when analyzing the blocks with and without viscous load separately, inter-dyad variability explained 99.4% (KL and KH) and 72.3% (in the four loaded conditions) of the variance. This means that the viscous loads changed how vigor scores are spread, thereby increasing intra-individual variability. However, when comparing conditions with similar dynamics, the variability in dyadic vigor scores was primarily accounted for by differences between dyads, as for lone movements. Therefore, each dyad seemed to behave as a single unit with its own preferred vigor.

To confirm this analysis, we then examined whether the relative dyadic vigor scores across the dyads changed with load and connection stiffness using Pearson correlation coefficients between pairs of conditions. All 15 Pearson correlations were significant (p < 2.3·10⁻³, Pearson-r > 0.64 in all cases), showing that the order of dyads according to their vigor scores was robust to evolving dynamics and effort levels (see Fig. 3I and Supplementary Fig. S.5J for examples of correlations). This confirms that vigor scores mainly vary across dyads rather than between conditions.

Therefore, we conclude that goal-oriented movements of mechanically connected individuals are characterized by a dyadic vigor, supported by: (i) coordinated motor plans between partners, evidenced by minimal interaction effort and movement duration differences; (ii) kinematics similar to individual movements, with overall higher performance (time, smoothness and accuracy); (iii) affine amplitude-duration relationships; (iv) larger inter-dyad than intra-dyad variability; and (v) vigor scores robustly correlated throughout all conditions.

Emergence of dyadic vigor from the partners’ individual vigor

Connected movements duration were either below (unloaded conditions) or close (loaded conditions) to the movements durations of the fast partners (Fig. 3H and Supplementary Fig. S.5G–I), suggesting that the faster partner primarily drives the dyad and selects dyadic vigor, thereby acting as a leader. This rules out a LF hypothesis where the leader would be the slower partner. On the converse, the LF hypothesis with the fast partner being the leader, where the slow individual acts as a follower simply aligning their motor plan to the leader based on interaction forces, remains possible although it would difficultly explain the observed increase of movement velocity in unloaded conditions. Alternatively, the fast partner may be the leader, but could still marginally account for the vigor of the slow partner to set the dyadic vigor, for instance based on a convex sum with a high weight on their own vigor. According to this WA hypothesis, both partners’ individual vigor should be related to the final dyadic vigor as follows,

where v^d is the dyadic vigor, v^f and v^s are the fast and slow partner’s vigor, and α allows to adjust the importance given to each partner. Specifically, (i) α = 1 and α = 0 respectively fall back to the LF hypothesis with a fast and a slow leader, (ii) α ∈ (0.5, 1) corresponds to a fast leader accounting for the slow partner’s vigor, or the converse with α ∈ (0, 0.5), and (iii) α = 0.5 corresponds to an equal contribution of both partners. Given the results obtained on movement durations, one can assume α should be close to 1 in our dataset.

To test the validity of this hypothesis, we used a linear mixed model (LMM) that examined how the dyadic vigor was predicted by the vigor of the fast and/or slow partner throughout conditions (see Eq. 4 in Methods). Surprisingly, this analysis revealed that the dyad’s vigor was related to the fast partners’ vigor only in KLVL (r = 0.14, p = 0.03) and not in any other condition (KL & KH: r < 0.14, p > 0.22; KHVL: r = 0.13, p = 0.15; KLVH & KHVH: r < 0.1, p > 0.24). On the contrary, the slow partners’ vigor was a significant predictor of the dyadic vigor in five out of six conditions (KL & KH: r > 0.624, p < 0.001; KLVL & KHVL: r > 0.24, p < 0.009; KLVH: r = 0.2, p = 0.028).

This surprising finding invalidates the LF hypothesis, which implies that the vigor of a dyad should be directly related to the individual vigor of the fast partner and not the vigor of the slow partner. Furthermore, the LMM analysis also shows that the dyadic vigor does not depend on mixed parameters such as the difference in vigor between partners or a convex sum of both partners’ vigor as in Eq. 2, as they imply that the vigor of the fast partner should be a significant predictor of dyadic vigor. This invalidates the WA hypothesis, which assumes an essentially similar contribution of the two partners to the dyadic vigor.

We are left with the IA hypothesis, according to which both participants contribute to dyadic vigor by adapting their own strategy to the new dynamical context. While dyad’s individual may not be able to perfectly infer the partner’s motor plan, they may identify the distribution of the partner’s movement duration and use it to plan their motion under this timing uncertainty. In particular, both participants may partially identify each other’s motion plan and timing and optimize their behavior accordingly so as to minimize interaction efforts variability, as they otherwise generate instability. This uncertainty induced by the partner’s timing, the effects of which accumulate through time, could cause an individual to move faster than initially planned. This is because increasing speed have been shown to reduce the detrimental effects of uncertainty on performance [26]. This dynamic change of vigor could lead to a new vigor of the dyad, which could in principle capture the observation that dyads move even faster than the faster partner in unloaded conditions.

Computational modeling of coordination mechanism

To test this premise, we designed a computational model representing as accurately as possible this dynamic adaptation related to the partner’s uncertainty. Importantly, this model had to reproduce all the quantitative and qualitative observations reported above. In particular, it was based on the above-reported observations to (i) minimize variations of interaction efforts, as dyads were shown to coordinate their movements to allow the emergence of dyadic vigor, (ii) predict that dyads move faster than either partners in solo when there is no load, (iii) with the dyads’ movement times close to those of the fast partner in solo in presence of load, (iii) predict a dyadic vigor that depends on the vigor of the slow partner and is not significantly influenced by the fast partner, and (iv) result in overall small dyadic movement duration prediction errors across the 30 tested conditions (5 targets × 2 stiffness levels × 3 viscous loads).

We simulated a dyad that would adapt their behavior relative to an inferred distribution of the slow partner’s movements (see Methods for details), while integrating costs of the fast partner in the planning. Therefore, this model assumes an asymmetric role for the fast and slow partners, where the fast partner costs are used by the dyad to plan the movements and the uncertainty arising from the interaction with the slow partner acts at a higher level through the interaction torque. First, we extract an optimal control model of effort and time characterizing the participants’ average behavior from data of the solo session (as explained in the Methods) [14]. Using the identified cost function, we can predict the movement duration and mechanical work of the average lone participant in the VL and VH conditions based on their P1 behavior (see Supplementary Fig. S.2). This identified individual motor plan then serves as a basis for simulating the interaction within dyads. Specifically, the effort cost identified for an average subject is used to build the model associated with IA, while the time cost is identified separately for the slow and fast group (see Methods for values of fitted parameters). To model the uncertain component of the interaction, the slow partner is represented by a distribution of minimum jerk trajectories of stochastic duration with average µ_k and deviation σ_k for target k, based on their average vigor. The dyads’ common plan can then be computed through the minimization of an expected cost function, including the common effort cost and the individual cost of time of the fast partner, as well as the variations of interaction torque :

where 𝔼 is the expectation with respect to the slow partner’s distribution, the cost function including the common cost of effort and the cost of time of the fast participant (see Eq. 9 in Methods), Q_τ a weight and T_k the free final time minimizing a positive function of the interaction torque variation , associated with the distribution of the slow partner’s movements duration. The computation and optimization of this cost function are described in the Methods.

The model’s predictions are summarized in Fig. 4. They correspond well to the experimental data in terms of the shape of the velocity profile, which is bell-shaped and smooth, and in terms of the predicted movement duration. The median prediction errors across the 30 simulated conditions were small (KL: 28.7 ms; KH: 50.8 ms) with most errors below 80 ms except for A₁ with viscosity, A₂ with KL and viscosity, and A₃–A₅ for KH without viscosity. Importantly, all these predictions were obtained with the same weighting Q_τ of the cost function in Eq. 3 (or Eq. 16 in the Methods for details).

Figure 4:

How does the IA model work? Co-adaptation of the two partners yields an average dyadic movement duration close to the average preferred duration of the fast group, but the slow partner ultimately determines which dyad is faster or slower relative to this average behavior. Importantly, the uncertainty associated with the distribution of the slower partner’s motor plan may explain the faster movement in KL and KH. Specifically, the partners would try to minimize the accumulation of errors in expected variations of interaction efforts to stabilize the interaction, resulting in faster movements as this cost naturally increases with time (see uncertainty term in Eq. 16 in Methods). In this context, the slower partner would only be willing to incorporate this acceleration up to a certain point depending on their individual vigor, thereby determining the final dyadic vigor. Interestingly, this would also explain why this global acceleration is not observed in presence of a viscous load. Although such efforts may increase the deviation of the inferred partner’s distribution by masking interaction efforts, they also naturally increase stability by dissipating energy. Consequently, they automatically reduce the variations in interaction efforts, and their associated detrimental effects on dyadic dynamics, thereby reducing the cost of increased movement duration. In our model, this important property naturally emerges from the compromise between time and effort, as identified for individual behavior, and uncertainty, due to the interaction with the partner, when controlling the dyad’s dynamics.

Supplementary Figs. S.17,S.18 present the results of a sensitivity analysis varying the model’s three parameters: weight of interaction torque variation minimization and the slow partner’s movement time distribution average and deviation . This analysis shows that (i) the model remains stable beyond the durations and parameters examined in the present study, and (ii) dyadic movement duration decreases as µ_k decreases for large ranges of parameters, consistent with the observed correlation between dyadic vigor and the slow participant’s vigor. We also verified the effect of changing the cost of time of the fast partner, i.e. their vigor, on predicted movement duration (see Supplementary Fig. S.19). This analysis shows that modifying the fast partner’s vigor only marginally impacts the optimal predicted time, with changes of 10 ms and 35 ms for a 30% increase in their cost of time for the A₁ and A₅ targets respectively. These minor impacts are consistent with the fast partner’s vigor not being a predictor of dyadic vigor. IA is therefore the only hypothesis that captures all the qualitative features of the experimental data, and the associated model provides accurate movement duration predictions, thereby complying with all the requested conditions for computational modeling summarized above. Furthermore, as the impact of the fast partner’s vigor on the predicted time is small, the predicted motion plan can be considered as representing the dyad’s emerging common plan, rather than a plan that would result from a fast leader that would just account for the slow partner’s behavior.

2.3 After effects of dyadic interaction

Did dyadic practice influence individual movement timing? To address this question we compared the movements of slow and faster partners before (P1) and after (P2) dyadic practice. Fig. 5A shows that averaged trajectories across the population were notably different in P2 compared to P1. In particular, the average movement durations decreased, while the velocity profiles remained similar to those with both KL and KH, without terminal submovement. Furthermore, the average amplitude-duration relationship in P2 changed for the slow partners but not for the fast (Fig. 5B).

Figure 5:

A one-way ANOVA confirmed behavioral differences between groups and conditions (SS> 0.3, p < 10⁻³, η² > 0.08 in both cases). In P1, the slow partners moved significantly slower than the fast partners (p < 10⁻³, Cohen’s D = 0.72) but not in P2 (p > 0.72). Furthermore, the group of slow partners changed their movement durations between P1 and P2 (p < 10³, Cohen’s D = 0.8) but not the group of fast partners (p = 0.41). Movement amplitude had a main effect on time saving (F_4,76 = 115.8, p < 10⁻³, η² = 0.36), which was confirmed across the population for targets {A₃, A₄, A₅} (p < 0.027, Cohen’s D > 0.51 in all cases, Fig. 5C). The vigor analysis (Fig. 5D, Supplementary Fig. S.20A) confirmed this adaptation: while the slow partners’ vigor increased between P1 and P2 (p < 0.001, Cohen’s D = 2.73), the fast partners’ vigor did not change (p = 0.59).

Interestingly, despite the change of motor plan of the slow partners, the P1 and P2 vigor scores over the whole population remained correlated (p = 3.2·10⁻³, r = 0.62, see Supplementary Fig. S.20B). Furthermore, the vigors of the group of slow partners in P2 were correlated to their behavior in P1 (p = 0.001, r = 0.87). Surprisingly, the vigors of the group of fast partners in P2 was correlated with the vigor of their slower partner in P1 (p = 0.017, r = 0.73), but not with their own vigors in P1 (p = 0.14, r = 0.5). This highlights a structural change in the repartition of vigor scores in the fast partners, reflecting their adaptation to the slow partners in each dyad.

Therefore, dyadic practice led to an adaptation of the average vigor in the slow partners, but their relative repartition was not modified. That is, fast participants within the group of slow partners remained fast after dyadic practice, and slow participants remained slow. In contrast, the average vigor scores of the fast partners remained unchanged in P2 relative to P1. However, there was a structural change in their repartition according to vigor, which matched the repartition of their slow partner after dyadic practice (Fig. 5E). Importantly, at the dyad level the adaptation is sometimes larger for the slow participant and sometimes larger for the fast participant. For instance, the dyad composed of subjects S5 and S6 shows a large decrease of vigor for the fast (S6) and an increase of a lower magnitude for the slow (S5) participant after connected practice (see Supplementary Fig. S.20).

3 Discussion

The movements of individuals are characterized by their vigor, revealing how prone individuals are to perform actions fast or slowly in various conditions. Vigor is a fundamental idiosyncratic trait of volitional movements that relates to energy, time/reward, and performance. How do individuals combine their vigor to coordinate movements when physically collaborating on joint actions? We addressed this question by investigating how partners connected by a virtual elastic band performed goal-directed hand movements without time constraints. Unlike previous human interaction studies, where timing was externally prescribed, either by a tracking task [18, 19, 27], explicit instructions [16], or the introduction of an explicit reward when connected [15], our study allowed for natural and unbiased timing emergence.

We first observed that connected partners used minimal interaction force and had negligible time differences in movement execution. This suggests that partners did not move independently, as proposed by [16], but instead adopted a common motor plan. Importantly, this implies that sensory exchanges occur between the partners even in simple point-to-point movements not only for complex scenarios as was shown previously [27–30]. Second, in conditions without viscous load, dyadic movements were strikingly faster and smoother than the movements of either of their individuals, without compromising accuracy. With a load the group of slower partners moved as fast as the faster partners, while the group of faster partners essentially maintained their pace. These findings confirm and extend the benefits of dyadic control observed in collaborative actions [15, 27]. Furthermore, dyads exhibited consistently slower or faster movement duration compared to the average over dyads, which was shown for different levels of effort and connection stiffness. As dyads’ partners coordinated their actions with these systematic features, this allowed us to define dyadic vigor similarly to individual vigor [11, 22].

We found a striking result regarding this dyadic vigor, which we expected to primarily correlate with the fast partners’ vigor in a leader-follower scheme [17]. While the connected participants had a movement duration close to that of the group of fast partners, dyadic vigor was instead predicted by the slower partner’s vigor. This unexpected finding was consistently observed across varying stiffness and resistance conditions, where the slower partners’ vigor was a significant predictor of dyadic vigor. This reorganization of the fast partners’ vigor scores highlights the slow partners’ pivotal role in determining the dyad’s movement pace. While the slower partners move faster connected than alone, the faster partners adapt their effort accordingly to match the slow partners’ influence. Note that the magnitude of the adaptation can be large for both of them, sometimes larger for the faster.

This finding challenges the leader-follower framework often used to describe human-human and human-robot interactions, where the leader would typically be the individual exerting greater motor force or being ahead in the movement [15, 18, 19]. Our results show that this interpretation is misleading. It also challenges the weighted averaging hypothesis suggesting that the dyad vigor should be systematically in-between the slow and fast partners’ initial vigor, with both contributing to its setting. Rather, our results demonstrate that a dyad motor plan emerges from a co-adaptation where the slower partner’s vigor is the final determinant. Instead of a fixed leader-follower relationship or weighted adaptation, human-human interactions appear to depend on the context rather than on predefined roles or a simple mix of the initial individual strategies [31]. This should be considered when designing human-robot interaction controllers [32], contrary to the numerous control methods relying on a leader-follower strategy to provide assistance [17].

Following connected training, strong after effects were observed in both fast and slow partners. For slow partners, a large increase in the individual vigor was observed. Interestingly, the vigor scores of the fast partners were reorganized depending on their partner, while the order of the slow partners’ scores was unchanged by dyadic practice, similarly to dyadic vigor during the connected session. This highlights how both partners adapted. The slow partners increased vigor, converging to the average vigor of the fast partners, while the fast partners reorganized their vigor score corresponding to the repartition of the slow partners. These effects persisted across 100 simple reaching movements, whereas after-effects of force fields or visuomotor rotations are typically re-optimized within a few trials [33, 34], and effects of previous attempts at modifying individual vigor tended to vanish through unconstrained practice [35]. This persistence indicates a potential change in the time perception for both groups after practicing with a partner.

These findings suggest joint training as a potential therapeutic approach to restore declining movement vigor, such as in individuals with Parkinson’s disease [36, 37]. In this context, future works should examine whether the vigor changes observed after dyadic practice are retained over longer periods or with intensive movement sessions, or whether they end up vanishing as was observed in studies where participants were instructed to move faster using visual cues [35]. Fundamental changes in vigor may be assessed by analyzing the dopamine levels in the basal ganglia [38, 39] (e.g. using positron emission tomography [40]), as dopamine plays a major role in vigor regulation [6, 7, 9, 10, 41, 42].

We developed a computational model based on stochastic nonlinear optimal control with free final time, thereby combining several computational approaches [2, 3, 43], in which dyads plan movements based on the slower partner’s average vigor under uncertainty, as was observed during connected movements. This model successfully predicts the key qualitative features of our data, in contrast to hypotheses that assume i) independent motor plans for slow and fast partners (CA), ii) a vigor distribution decided by a leader and accepted by a follower (LF) or iii) a negotiation based on a weighting of their initial strategies in the task (WA). These results further show the role of the slow partner and the critical impact of uncertainty with respect to the partner’s plan in the emergence of dyadic vigor. These findings argue in favor of a dynamic adaptation of an individual’s motor strategy caused by the presence of the other partner (and its inherent induced uncertainty).

Our model achieved minimal movement duration prediction error across all 30 tested conditions, spanning five movement amplitudes and six load conditions. Furthermore, we conservatively assumed that participants used a single cost function across all conditions, even though humans can adjust strategies based on task demands [44]. A sensitivity analysis confirmed that increasing the slow partner’s average movement duration leads to an increase in the optimal dyadic movement duration, whereas changing the fast partner’s cost of time has only marginal effects on the predicted dyadic vigor, aligning with our findings on the correlation between slow partners vigor and dyadic vigor.

In that process, we have shown that the faster individual in a dyad relies on knowledge of the partner’s movements distribution, thereby providing a plausible mechanism to explain the faster movements without viscous load, the preserved or decreased average movement duration of the fast partners in presence of viscous load, and the surprising dependence of dyadic vigor on the slower partner. This model may be used to develop robotic systems interacting with humans in a manner similar to human partners, with applications including collaborative robots for manufacturing [45], rehabilitation robotics [46], exoskeletons [47] and supernumerary robotic limbs [48].

4 Methods

4.1 Participants and experimental setup

The experimental protocol was approved by the Imperial College’s Science, Engineering and Technology Research Ethics Committee (SETREC number 7112072). N = 20 right-handed young adults (4 females) with no known motor or cognitive impairment participated in the experiment in 10 dyads of two partners. Their biographical data were: age 27.5 ± 5.5 years old, weight 73.7 ± 14.5 kg, height 175 ± 8.5 cm and right hand length 19.1 ± 1.2 cm (from the wrist flexion/extension rotation axis to the farthest fingertip). Each participant was informed about the experiment, and signed an informed consent form before participating in it.

Each of the dyad’s partner interacted with an HRX-1 robotic interface (HumanRobotiX, London, UK) using their right wrists, while seated in front of their respective monitor (Fig. 1B). Participants had to control a 0.6 cm diameter cursor towards a 2.7 × 2.7 cm² target via wrist flexion/extension. These dimensions ensured consistent movement amplitudes across participants with minimal accuracy demands. The two HRX-1 were controlled and synchronized at 1 kHz using a TI PiccoloF28069M LaunchPad. Wrist flexion/extension angles were recorded at 300 Hz with an optical encoder (MILE 512-6400, 6400 counts per turn). Angles (18^°, 36^°, 54^°, 72^°, 90^°) were linearly mapped to targets (A₁, A₂, A₃, A₄, A₅) on the monitor, with 9.3 cm between the centers of consecutive targets (Fig. 1B).

4.2 Experimental protocols

Participants completed two experimental sessions (Fig. 1C):

• In a first solo session, the two participants performed goal-directed wrist flexion/extension movements independently. These movements were used to: (i) verify that our population sample followed a vigor law, (ii) estimate their individual vigor and cost of time for computational modeling, and (iii) establish a baseline to analyze after-effects of the human-human interaction.

• In the subsequent dyadic session, they were coupled via a virtual elastic band implemented through the robotic interfaces. The data enabled us to investigate movement coordination between the partners and its dependence on individual vigor.

The participants performed back and forth reaching at preferred speed, starting with a wrist extension movement from A₀ to one of the targets {A₁, A₂, A₃, A₄, A₅}, followed by a flexion movement, back to A₀. After each movement, the timing of the next movement start was drawn from a uniform distribution to prevent anticipation and rhythmic behavior. After holding at A₀ for 2 s ± 0.15 s, A₀ was switched off and a target was switched on for 2 ± 0.15 s, then A₀ was switched on again. This cycle was then repeated.

Solo session

Participants completed 3 blocks of 100 trials (10 extension and 10 flexion movements per amplitude: {18^°, 36^°, 54^°, 72^°, 90^°}) in a fixed pseudo-random sequence. A two-minute break was provided between blocks to prevent fatigue. The first block was performed in “transparent” mode (robot motor off). The next two blocks introduced resistive viscous loads (ν ∈ {0.075, 0.15}Nm s/rad in random order) to examine vigor robustness under varying effort costs. A final washout block was introduced to avoid carrying after-effects of the two visco-resistive blocks into the connected session. This block was shorter than the other blocks, consisting of 40 reaching movements with 4 flexions and 4 extensions per amplitude. The same amplitudes and targets were used as for the other blocks but their order was not randomized, i.e. the blocks consisted of 4 series of movements towards A₁, then A₂, …, then A₅.

Dyadic session

Dyads of participants completed 6 blocks of 100 trials while coupled via a virtual visco-elastic band. The first two blocks used an elastic band of stiffness κ∈ {0.5, 1.6} Nm/rad in randomized order. The next four blocks tested all κ and ν combinations in random order to explore time-effort tradeoffs. A seventh block without coupling (i.e. κ = ν ≜ 0) was used to assess carry over effects of the interaction with the partner on individual vigor.

4.3 Data processing

Wrist kinematics

The wrist velocity and acceleration were computed by numerical differentiation, after low pass filtering of position data using a forward-backward fifth order Butterworth filter with 5 Hz cut-off frequency. Wrist flexion and extension movements were first grossly segmented based on the timings of the targets, then a threshold of 5% of peak velocity was used to compute the start and end of each movement. Leftwards and rightwards movements were pooled together as their vigor was similar, consistently with expectations [3, 22].

The amplitude-duration relationship of each participant was obtained by i) averaging movement durations and amplitudes per target, ii) performing an affine least-square regression (using the LinearRegression function of scikitlearn.linear_model). The amplitude-duration relationship of a group of participants was obtained by: i) averaging across the group the average movements durations and amplitudes of each participant, for each target separately, then ii) performing an affine regression.

All the movement durations were obtained from affine fittings of the amplitude-duration relationships. The vigor scores in the VL and VH conditions were computed using the average movement duration of the population for each target from P1, which allowed to illustrate the general decrease in vigor with a resistive torque. The same approach was adopted for the dyadic session, using the population averaged movement duration in KL and KH to compute the vigor scores in {KLVL, KLVH} and {KHVL, KHVH}, respectively. As a consequence, the vigor scores were equally distributed around 1 only in the P1, KL, and KH conditions. To analyze the variance, the vigor scores were computed based on the average movement duration of the population in each condition.

The reaction time of participants was computed as the time difference between the apparition of a target and the first instant where the wrist angular velocity exceeded 0.01 rad/s.

4.4 Statistical analysis

The statistical analyses were performed using custom Python 3.8 scripts, the statsmodel package [49], and the Pingouin package [50]. For all metrics, the results were first averaged in each participant, condition, and movement amplitude before performing statistical analyses. Normality (Shapiro-Wilk test [51]) and sphericity (Mauchly’s test [52]) of the distribution of the residuals were verified. Main effects of the condition and amplitude were then assessed using repeated measurements ANOVA, and one-way ANOVA for comparisons between the fast and slow partners. Greenhouse-Geisser corrections were applied if sphericity was violated (ϵ < 0.75). The significance level of the ANOVA was set at 5%.

In case of a significant main effect, post-hoc pairwise comparisons were performed using paired t-tests within subjects, and unpaired t-tests between the fast and slow partners. The significance level of post-hoc comparisons was set at 5%. All the significant comparisons are reported with the Cohen’s D value to illustrate the effect size, where D < 0.4 would be considered as small and reported as a weak result.

The prediction of the dyads’ vigor in the connected conditions was performed using a linear mixed model attempting to correlate the vigor of the dyad to the vigor of its fastest and slowest participants as

where v^d is the vigor of dyads in one of the connected conditions {KL, KLV L, KLV H, KH, KHV L, KHV H}, v^f is the vigor of the dyad’s faster partner in the P1 condition and v^s of the slower partner.

4.5 Computational modeling and simulations

4.5.1 Solo session

Data collected in solo setting were used to model the dyad behavior and predict the dyadic behavior from individual movement characteristics. The solo session included three levels of viscous resistance ν ∈ {0, 0.075, 0.15} Nm s/rad (referred to as P1, VL and VH, respectively), resulting in the task dynamics

where is the wrist angular acceleration, the angular velocity, τ [Nm] the human wrist torque, I = 9.2·10⁻³ Nm s² the population average hand+robot inertia, and D = 0.03 Nm s is the wrist damping selected in the range of values identified during passive movements in [53, 54].

Following the inverse optimal control procedure of [3, 14], we first estimated the cost of time of the average participant based on the population average amplitude-duration relationship (see Fig. 2D) in P1 without viscous load. The estimated amplitude-duration relationship for the population was

Moving the hand from an initial angle q₀ to q_e was formulated as an optimal control problem, using the minimum torque change framework [55]. Utilizing the state and the motor command (as in [3]), the cost of effort to minimize for amplitude k was

where is the final time of the average participant for amplitude k computed from Eq.(6). Q_k = diag(15.5/A_k, 0, 0) and β = 0.95, determined by trial and error, allow to adjust the compromise between minimizing the torque change and the distance to the target. As detailed in [14], the human cost of time G(t) can be identified by fitting and integrating a sigmoid curve on the Hamiltonians of the fixed final time problems from Eq. 7, yielding

where p = (p₁, p₂, p₃, p₄) = (1.5713, 6.9723, 0.4497, 1) was identified for the average participant in our sample. Using the same effort formulation as in Eq. 7, it is then possible to predict the movement duration and cost of effort for ν ∈ {0.075, 0.15} Nm s/rad by minimizing

in free final time, i.e. when the movement duration is optimized simultaneously to the effort cost. Here V_k is the cost to minimize for amplitude A_k and T_k ∈ [0.01, 10] s the optimal final time for amplitude k. Note that minimizing V_k with ν ≜ 0 leads to , and to the same predicted motor behaviors as the fixed final time optimal control problems minimizing Eq. 7. To compute the effort cost associated with a constant time strategy when subjected to the resistance of the exoskeleton, we also simulated the case with ν ∈ {0.075, 0.15} Nm s/rad.

The simulations of the solo session were performed using the Matlab (R2021a, MathWorks) version of GPOPS-II [56–58], which is based on an orthogonal collocation method relying on SNOPT to solve the nonlinear programming problem [59].

4.5.2 Dyadic session

The task of the dyadic session included two levels of connection stiffness κ ∈ {0.5, 1.6} Nm/rad, and three levels of viscous resistance as for the solo session, yielding the coupled dynamics of the fast and slow partners

where the superscripts f and s indicate the fast and slow participants, respectively, τ ^f and τ ^s are their wrist torques, I is the total inertia of the wrist+exoskeleton, D the total damping as in the solo session, and are the wrist angular acceleration, speed and position, respectively.

To simulate movements according to the CA hypothesis, we split the participants in the groups of the faster and slower partners in each dyad according to their vigor in P1. Furthermore, to simulate the IA hypothesis, assuming a coordination of the partners’ motor plans, we had to identify the cost of time of the fast and slow partners, which was done based on their respective averaged amplitude-duration relationships

Importantly, to simulate IA we identified the cost of time of the fast (p = (3.9615, 5.8958, 0.4459, 1)) and the slow partners (p = (1.0184·10⁷, 11.2958, 1.9980, 1.4447·10⁻⁵)) separately, which allowed us to simulate the effects of the average partners’ individual vigor on the dyad’s motor plan (see Eq. 8 for the definition). Then, the simulated hypotheses induced different state variables, conditions and numerical methods to predict movements under the coupled dynamics as described below.

Co-activity (CA)

assumes that the participants do not coordinate their movements, with each partner trying to perform their predefined desired trajectory (as suggested in [16]). For the sake of simplicity, we set these predefined plans using simple minimum jerk plans with durations identified in the lone session for the fast and slow participant separately. Optimal control commands and to track the obtained desired trajectories for the fast and slow partners for target k, respectively, are computed using a linear quadratic regulator. We use the state and the cost matrices Q = diag(100, 0.1, 100, 0.1) and R = 𝕀₂. The fast partners’ trajectory is extended with null values to match the duration of the slow partners’ plan. The resulting coupled movements are Euler integrated using these control inputs in Eq. 10.

Interactive adaptation (IA)

Here we assume that dyads’ partners collaboratively optimize their movements based on the solo behavior of the slow partner, who “drives” the co-adaptation. Our model describes the dyad’s motion plan using a single cost function. While one would intuitively seek for multiple agents modeling, to our knowledge there is no tractable solution to simulate multiple agents under uncertainty, in free final time, and not based on usual quadratic cost functions. Both partners contribute to selecting the final movement control by optimizing the common motor plan based on a prediction of the slow partner’s motion, which is uncertain for the fast partner. For the sake of simplicity, we represent the slow partner’s trajectory as a distribution of minimal jerk trajectories [60]

where the movement duration is a random variable with mean µ_k and standard deviation σ_k. Note that q₀ and the {A_k} are considered common and known to both partners, and no assumption is made on ξ_k’s distribution beyond its deviation. The task dynamics inferred by the fast partner are

where is the interaction torque for target k. Since ξ_k is a random variable, the uncertain interaction torque change can be approximated based on a linearization of the velocity of the slow partner around µ_k:

yielding

The expected cost function of the dyad can the be expressed as a tradeoff between the fast partner’s cost function when moving alone (including their cost of time) and the uncertain interaction torque change with respect to the behavior of the slow participant, reflecting a strategy of stabilization of the interaction efforts. We consider that participants would minimize a cost related to the expected interaction torque and its uncertainty, while neglecting the cross-covariance terms, expressed as in the main text:

This leads to the expected cost

where is the dyadic motor command, is the variance of ξ_k and the velocity component of the state covariance matrix. This stochastic optimization problem can be approximated under a deterministic form and solved within the stochastic optimal open-loop control framework of [43] (see in particular Eq. 3), which provided codes on which we based our simulations, using Julia [61], JuMP [62] and Ipopt [63], where the problem was transcribed into an NLP problem using a trapezoidal scheme [64].

In the simulation results of Fig. 4B corresponding to the IA hypothesis, the dyads were considered to use the same cost function for all amplitudes and conditions with Q_τ = 3. The uncertainty with respect to the slow partner’s behavior was assumed to increase with movement amplitude and resistive viscosity. Specifically, we used values uniformly spread within σ_k ∈ [0.15, 0.25] s for KL and KH, within σ_k ∈ [0.2, 0.3] s for KLVL and KHVL, and within σ_k ∈ [0.32, 0.4] s for KLVH and KHVH.

Data availability

Data are publicly available on Zenodo.

Acknowledgements

We thank Jingwen Zhao for her help during the development of the experimental setup. We thank Claudia Clopath for her comments on prior versions of the manuscript. We thank Dario Farina for his help in funding the present study. This study was funded by the European Research Council (ERC Synergy, “Natural BionicS”, grant number: 810346), and by the European Commission (FETOPEN H2020, “NIMA”, grant number: 899626 & Leadership in enabling and industrial technologies - Information and Communication Technologies (ICT), “ReHyb”, grant number: 871767).

Additional information

Contribution statement

• Conceptualization – all authors

• Methodology – all authors

• Software – DV, BB

• Validation – DV

• Formal analysis – DV, BB

• Investigation – DV

• Resources – EB

• Data Curation – DV

• Visualization – DV, EB

• Supervision – EB

• Project administration – EB

• Funding acquisition – EB

• Writing - Original Draft – DV

• Writing - Review & Editing – all authors

Funding

EC | European Research Council (ERC) (810346)

• Etienne Burdet

EC | H2020 | PRIORITY 'Excellent science' | H2020 Future and Emerging Technologies (FET) (899626)

• Etienne Burdet

EC | H2020 | IL | LEIT | H2020 LEIT Information and Communication Technologies (ICT) (871767)

• Etienne Burdet

Additional files

SupplementaryFile