The advertisement of the study participants was posted on a web-board of Seoul National University. The inclusion criteria comprised no medical history of neurological disorders, injuries of upper extremities, and vision. In particular, the experimental tasks required the proper acquisition of visual information on the computer screen, all volunteers who passed the Freiburg Visual Acuity and Contrast Test (FrACT; Bach, 1996) were selected for participants (i.e., above 1.0 decimal scale as the normal condition of vision by 20/20 visual acuity scale). A power analysis using G*Power (Faul et al., 2007) was performed to estimate a prior sample size, which suggested recruiting at least nine participants in order to reach an effect size (d) greater than 0.7 with at least 80% power and α = 0.05 as type-I error rate to detect significant differences between the conditions. Twelve volunteers were recruited, and three of them did not satisfy the inclusion criteria. Thus, nine right-hand dominant young males (age, 30.3 ± 2.7 years; height, 167.49 ± 6.53 cm; weight, 69.39 ± 15.73 kg) participated in the experiment. The Seoul National University Institutional Review Board (IRB) approved the use of customized experimental protocols and compatible devices related to the current behavioral tasks (e.g., force transducers and fNIRS devices). After providing information about the study, we requested all the participants to sign a consent form approved by the IRB at Seoul National University (IRB No. 2007/002-028). The original signed consent form was retained in the experimental records, and a copy of the signed consent form was provided to the participants.

Four force transducers (Nano-17, ATI Industrial Automation, Garner, NC) were attached to a customized experimental frame (140 × 90 × 5 mm³) to independently measure the pressing forces (i.e., z-axis forces) produced by each of four fingers (Figure 1B). The transducers’ surfaces were covered with sandpaper to provide sufficient friction to the fingertips. There were four slots in an anterior-posterior direction on the panel to adjust the transducer position to the hand and finger size of individual participants. The mediolateral distance between slots was fixed at 3.0 cm (Figure 1B). The 3.0 cm value was a typical distance reported in the previous hand/finger pressing experiments, which possibly ensured a comfortable hand posture during pressing with fingers (Park and Xu, 2017; Kim et al., 2018; Kong et al., 2019). The experimental frame attached to the transducers was mechanically fixed on an immovable table. Four analog signals from the transducers were digitized with 16-bit analog-digital converter (USB-6225, National Instrument, Austin, TX, USA) via a customized LabVIEW program (LabVIEW 8.0, National Instrument, Austin, TX, USA). The sampling rate of force signals was set to 200 Hz.

Depth-dependent hemodynamic responses in the prefrontal cortex were recorded using a wearable fNIRS device (NIRSIT, OBELAB, Seoul, Korea) at a sampling rate of 8.138 Hz. The fNIRS probes were in contact with the forehead (Figure 1C). The effect of ambient light was minimized, as confirmed by monitoring the channel quality (refer to the data processing section for more details) throughout the experiment. The length of the source-detector (SD) was 30 mm as in previous reports (Scholkmann et al., 2014), and a total of 48-channel signals covering the majority of the prefrontal cortex were measured.

All experimental procedures were performed in accordance with the relevant guidelines and regulations of the IRB. Participants sat in a height-adjustable chair with a computer screen placed at eye level. The right upper arm was positioned in the wrist-forearm brace and strapped with Velcro to prevent excessive forearm and wrist movement during the production of pressing finger forces. A wooden piece was placed underneath the participant’s palm to ensure a fixed finger configuration during the experiment (Figure 1A). Before each trial, the subjects were asked to put their fingertips on the center of the corresponding transducers. Prior to each trial, all transducer signals were set to zero at a relaxed moment; thus, the gravitational effect was excluded, resulting in only active pressing finger forces being measured during data acquisition.

The experiment consisted of auxiliary finger force production tasks as well as the main task as a cyclic–force-production task in an isometric condition. The auxiliary tasks included maximal voluntary contraction (MVC) tasks, which involved single-finger pressing tasks and a multi-finger pressing task, and single-finger ramp-force-production tasks for individual fingers. The MVC pressing forces of the individual fingers or all four fingers (MVC_IMRL) were measured for each participant, and these values were further used to normalize the individualized target force values for the main task. The ramp-force-production tasks were designed to configure an interdependency matrix (i.e., enslaving matrix, E) by measuring unintentional force production by non-task fingers, according to the nearly linear relationship between changes in individual finger force and the total force (Latash et al., 2002a). In the ramp-force-production task, participants were asked to produce a force ramp pattern via one finger force from 5% to 25% of the finger’s MVC for 6 s after maintaining 5% of that MVC for 2 s. Participants were instructed to focus on the force produced by the task finger and on the template displayed on the computer screen and to ignore the force produced by non-task fingers while keeping all non-task fingers on the corresponding sensors.

The main task of the current study was cyclic-force-production using all four fingers. Participants were asked to produce a smooth sine-wave total finger force (F_TOT) while wearing the fNIRS device (Figure 1A). Note that the inherent time delay of the hemodynamics response was about 10 s (Fazli et al., 2012), the data acquisition of the fNIRS measure for each trial lasted until 10 s after the force measurement ended. The participants were asked to follow the real-time F_TOT on the template displayed on the computer screen while minimizing any unnecessary behaviors such as head motion and speaking that possibly affected the fNIRS measure. The fraction of the time window display on the screen was 10 s, from −5 s to + 5 s with respect to 0 s, that represented the real-time moment of actual force production. The force templates consisted of a 10 s steady force line with a 15% MVC_IMRL value, followed by an 80 s sinusoidal force template with three frequency conditions, including 0.1 Hz, 1 Hz, and 2 Hz and the same amplitude (i.e., from 10% to 20% of the MVC_IMRL value). The 10 s steady-state phase was designed to elicit a stable baseline hemodynamic state before cyclic force production (Bajaj et al., 2014). Each participant was allowed to practice for about 30 min, which possibly minimized the learning effect during the experiment. After the orientation session and sufficient rest, each participant performed three trials; they were allowed to rest for 5 min after the completion of a single trial. The three frequency conditions were block-randomized across all the participants.

A standard description of parametric statistics was used, and the data presented as means and standard errors. A parametric repeated-measures analysis of variance (ANOVA) with between-subject factors of frequency (three levels: 0.1, 1, and 2 Hz) and phase (two levels: F_DW and F_UP) was performed on the RP, (CV_peak and CV_trough), averaged ΔHbO at each of 48 channels, and FC_N (Hypotheses 1 and 3). All factors were selected for the particular statistical tests. Mauchly’s sphericity test was used to confirm or reject the assumptions of sphericity. Greenhouse-Geisser corrections were used when the sphericity assumption was rejected. The statistical power for all comparisons from the pool of nine participants was computed, and in all cases, the power was >0.7. For the pair-wise comparison, the changes in the magnitude of the estimated variables were presented using post calculation of the effect size (partial eta-squared, ηp2) and the lower and upper 95% confidence intervals (CIs). In addition, the normality of the measured data was assessed by the coefficients of curvature (skewness) and elongation (kurtosis) upon the statistical significance. The level of significance for all statistical tests was set at p < 0.05.

To test the influence of frequency on the time-series variables from the UCM computation, including V_UCM(t), V_ORT(t), and ΔV(t), one-dimensional statistical parametric maps (SPM) with repeated-measured ANOVAs (SPM{F}, Pataky et al., 2016) were performed for the F_UP and F_DW phases separately (hypothesis 2). SPM analyses were conducted using open-source software¹ (Pataky, 2012). A critical threshold for the SPM analyses was computed based on the random field theory (Worsley et al., 2004), which was set at α = 0.05. When the SPM curve crossed the critical threshold, the difference was deemed to be a significant effect of frequency at particular time points of the curves. In case of significant effects of frequency, follow-up SPM{t} paired t-tests were performed. In particular, linear regression models were applied to each sample of the normalized data set (i.e., 100 samples) to generate the SPM{F} curve across three frequency conditions and the SPM{t} curve for the pairs of two frequency conditions.

The group means and standard deviations of the RP determined using the pooled data over F_UP and F_DW phases were 2.68° ± 0.10°, 2.91° ± 0.16°, and 3.16° ± 0.20° under the conditions of 0.1 Hz, 1 Hz, and 2 Hz, respectively. These results imply that the two sets of cyclic force data were close to in-phase synchronization for all three conditions. Two-way repeated-measured ANOVAs performed with the factors of frequency (three levels: 0.1, 1, and 2 Hz) and phase (two levels: F_DW and F_UP) showed no significant main effects and factor interaction effects on the RP. The average CV values across subjects were negative for all three frequency conditions, representing force undershoots in most cases. In particular, the amount of force undershoots increased with frequency, especially for the CV_peak. These observations were supported by the Bayesian one-way repeated-measures ANOVAs, separately on CV_peak and CV_trough, with frequency as a factor, which confirmed a significant effect of this factor on CV_peak (F_(1,8) = 27.19, p < 0.01, ηp2 = 0.90), but not on CV_trough. Post hoc pairwise comparison further confirmed that the average CV_peak was −0.17 at 0.1 Hz [CI: (−0.23 −0.10), p < 0.05], > −0.42 at 1 Hz [CI: (−0.56 −0.29), p < 0.05], and −0.58 at 2 Hz [CI: (−0.73 −0.43), p < 0.05]. Lastly, the coefficients of skewness and kurtosis were −0.79 and 1.14, respectively.

The average time profiles of ΔV with standard deviations across subjects for the F_UP and F_DW phases are presented in Figure 3. For both F_UP and F_DW phases, a U-shape profile of ΔVs over time was observed under the 1 Hz and 2 Hz conditions, while ΔV was consistently higher under the 0.1 Hz condition than under the other two conditions (Figure 3A). SPM analysis with repeated-measures ANOVA performed separately on the F_UP and F_DW phases showed significant effects of frequency on ΔV over the entire time period (Figure 3B). In particular, the significant effect of frequency was prominent in the middle of the time period for both the F_UP and F_DW phases. Post hoc pairwise SPM t-tests confirmed ΔV at 0.1 Hz was greater than that at 1 Hz and 2 Hz for both F_UP and F_DW conditions, especially when SPM{F}crossed the critical value (Figure 3C).

Normalized time functions (i.e., 100-time points) of the two variance components, V_UCM(t) and V_ORT(t), were computed. The time profiles between the two variances were different, such that V_UCM(t) and V_ORT(t) were related to the magnitude of the total force (|F_TOT|) and the absolute value of the force derivative (|dF_TOT/dt|, i.e., inverted U-shape), respectively. To facilitate comparisons between the three frequency conditions, the condition with a larger frequency showed higher V_ORT values for both F_UP and F_DW phases (Figure 5A). However, little difference in V_UCM magnitude was observed between frequency conditions (Figure 4A).

The SPM analysis with repeated-measures ANOVA showed significant effects of frequency on both V_UCM (Figures 4B,C) and V_ORT (Figures 5B,C), with a stronger effect on the latter than on the former. In particular, a strong effect of frequency on V_ORT was observed in the middle portion of the period for both F_UP and F_DW phases, where the peak values of the force derivative (dF/dt) were observed (Figure 5B). For the SPM{F} for V_UCM, a relatively strong frequency effect was observed when the force magnitude was small, such as in the latter portion of the F_DW and the early portion of the F_UP phases. Post hoc pairwise SPM t-tests confirmed that V_ORT trend for both F_UP and F_DW conditions was as follows: 0.1 Hz > 1 Hz 2 Hz (Figure 5C).

There were no significant effects on frequency and phase on ΔHbO for all 48 channels; on average, ΔHbO was around 0.3 μM (Figure 6A). Higher ΔHbO values were observed in the left hemisphere than in the right hemisphere, and the middle portion of the left hemisphere showed stronger ΔHbO for all three frequency conditions (Figure 6B). In addition, the FC_N under each frequency condition was determined by one-sample t-tests with the criteria of p ≤ 0.05 and r ≥ 0.7. Note that a prior analysis showed no significant difference between FC_N for F_UP and F_DW phases; thus, the data presented in Figure 7 and Table 1 were quantified using pooled data across F_UP and F_DW phases for each frequency condition. The number of significant pairs (i.e., FC_N) decreased with an increase in frequency (128, 50, and 61 under the 0.1 Hz, 1 Hz, and 2 Hz conditions, respectively; Table 1). Similarly, both intra- and inter-hemispheric FC_N decreased with increasing frequency, where the intra-hemisphere FC_N was larger than the inter-hemisphere FC_N (Table 1).

The main premise behind the concept of synergy is a crucial feature of all intentional movements performed by redundant systems that possibly represent the system stability (Latash, 2008a). Specifically, stable performance is coined to address the flexible patterns while compensating for errors among a redundant or, to be specific, abundant set of elements considering the mechanics of a particular task, whereby the solution family is organized in such a way that the actions of elements are varied to maintain the performance variables in an unchanged (or stabilized) state (Reschechtko et al., 2014; Ambike et al., 2016). One of the straightforward motivations of the current study is the curiosity as to whether the strategy and consequences of the organization of the solution family (i.e., combinations of four finger forces) change with different frequencies during cyclic force production. The current results of the inconsistent peak or trough force production (i.e., higher CV) were accompanied by a small synergy index at relatively high-frequency conditions, which is in line with the previous findings of the frequency dependence of various motor behaviors from the perspective of variability in muscle activation (Lewis et al., 2001; Huang et al., 2021), sensory information processing (Almeida et al., 2002), etc.

From a computational point of view, the synergy index refers to the relative magnitudes of V_UCM and V_ORT with respect to the total variance (V_TOT). In other words, less flexible combinations of elements (V_UCM) could be associated with a larger synergy index if the error variance (V_ORT) is significantly small. The current results showed that a strong statistical effect of frequency conditions was mainly observed in V_ORT, while the statistical effect on V_UCM was relatively weak, resulting in small synergy indices under high-frequency conditions. These results reflect the fact that the elements (i.e., force modes, m) had a dominant positive covariance at higher frequencies. The dominant positive co-variation between elements possibly causes inconsistent peak and trough force values, whereby performance errors increase accordingly. On the contrary, a strong negative co-variation between elements would serve as an error compensation strategy as to the total force stabilization, which was dominantly shown under the small frequency condition (0.1 Hz; Tseng et al., 2006).

However, we hesitate to strongly claim that the low values of synergy indices under high-frequency conditions represent weakened force stabilization. In Figure 5A, it was clearly observed that the time patterns of the orthogonal complement of variance (V_ORT) fit well with the rate of force changes (dF/dt, not shown in Figure), and this phenomenon is well described by the mathematical model, that is, the linear model proposed by Gutman (Gutman and Hagander, 1985; Latash et al., 2002b) and other experimental outcomes reported previously (Friedman et al., 2009). Because the direction of force changes results only from the combinations of elements within the subspace of the orthogonal complement to the null space (i.e., UCM space), the variance within the orthogonal subspace strongly depends on the scaling of the rate of force change. Thus, it is plausible that the increment of frequency naturally requires changes in the V_ORT, which results in less accurate (or more variable) force values at the peak and trough of the force profile. Contrary to our expectations, the results of the RP were close to the in-phase synchronization under all frequency conditions, which suggests that the phase difference between the produced and constrained force values was not the primary source of error variance. Of course, a less performance accuracy in the high-frequency condition is also caused by the inertial or biomechanical factors such as the bandwidth of muscle activation (Van Boxtel et al., 1998), electromechanical delay (Cavanagh and Komi, 1979), etc. In particular, the consideration of the inertial properties of body segments would be critical to avoid spurious interpretation of behavioral changes with frequency in dynamics. However, the relative contributions of the central and peripheral components to the performance accuracy remain unknown, even with the current results; thus, future research will have to ascertain the relationship between the inertial/biomechanical properties of the human body and neural response with various movement frequencies. Nevertheless, we have to find an alternative approach to answer the following question: What are the neural activities that cause changes in synergy indices for multi-finger force production with different frequencies of cyclic force production?

In the current study, we measured the hemodynamic responses in the prefrontal cortex while acquiring individual finger forces during cyclic force production. We quantified two indices of prefrontal hemodynamics, including ΔHbO and functional connectivity. Note that ΔHbO is indicative of the amount of localized oxygen consumption, and the functional connectivity represents significant connections (i.e., concurrent changes in oxygen consumption) between anatomically separated brain regions within the prefrontal structure. The current results showed non-parallel changes in the indices with the frequencies of cyclic force production, that is, O₂Hb concentration remained unchanged, whereas the functional connectivity (i.e., the number of significant connections) decreased with frequency. In particular, it might seem counterintuitive that the ΔHbO was not frequency-dependent, and these results are in stark contrast to previous findings regarding the frequency-dependency of cortical activations (i.e., the higher the frequency, the larger the cortical activation; Jenkins et al., 1997; Turner et al., 1998; Lutz et al., 2004), which presumably increased the oxygen concentration (Scholkmann et al., 2014). Oxygen is definitely a source of energy or cost for brain activities (McKenna et al., 2012). This difference between the results of the current and previous studies may have been caused by dissimilar experimental tasks and the focused region within the brain. Nevertheless, the unchanged index of oxygen consumption, particularly in the prefrontal cortex, is unlikely the “cost” for the current experimental tasks. On the contrary, concurrent prefrontal activation seems to be the frequency-dependent index for the current motor tasks, consistent with previous findings (Cordes et al., 2001; Salvador et al., 2008).

The current motor tasks required the subjects to use visual feedback. In other words, the values of cyclic force trajectories were prescribed through the template on the screen, and the subjects were asked to produce a total finger force such that the produced force with a cursor on the screen chased a prescribed target trajectory. For the low-frequency condition (e.g., 0.1 Hz condition), the subjects had enough time to utilize the feedback of visual or proprioceptive information, whereby it is possible that feedback-based force production is a critical factor in achieving a better performance accuracy under low-frequency than under high-frequency condition. The feedback-based error compensation strategy between finger forces may be less (or not) used during high-frequency conditions because of the limited time to use and process sensory feedback signals, which may result in a relatively large deviation and variance (i.e., error variance, V_ORT). Further, this speculation is supported by the dominant positive co-variation of finger forces (Goodman et al., 2005) under high-frequency conditions. By combining the results of hemodynamic responses and patterns of finger force co-variation, the error variance in a particular task could be linked to the functional connectivity of the prefrontal cortex, and not the magnitude of oxygen saturation. Indeed, it has been reported that the metabolism of brain activation with other brain imaging techniques (e.g., positron emission tomography, magnetic resonance imaging, etc.) showed a non-linear relationship with the degree of connectivity (Neubauer and Fink, 2009; Tomasi et al., 2014; Tomasi and Volkow, 2019) and that concurrent brain activation is caused by a combined effect of anatomical connectivity and synaptic efficiency (Buckner, 2010). Furthermore, unchanged oxygen saturation in the prefrontal cortex was observed during continuous movements with different properties (Jenkins et al., 1997; Kim et al., 2005), which is in line with the current results. However, the methodologies for the acquisition of brain activities in the current and previous studies are different, and functional connectivity is assumed to be an indirect quantification of brain activation. A possible interpretation, however, of the increase in functional connectivity, combined with unchanged oxygen consumption and small error variance for the stabilization of total finger forces in the low-frequency condition, is that the behavioral stability could be associated with the efficient organization of prefrontal activation.

This interpretation is in line with theories and computational models used for controlling a multi-element system. The computational model proposed by Todorov and Jordan (2002) (i.e., stochastic optimal control model) claims that the controller’s effort is not dedicated to the deviation of elements that do not interfere with critical mechanics but is more concerned with error correction than with the required mechanics. In other words, the energy/cost function may be associated with minimizing performance error rather than organizing variable combinations of solutions. This idea is conceptually similar to the UCM hypothesis based on the principle of motor abundance, which describes a crucial feature of the co-variation patterns between elements (Todorov and Jordan, 2003; Martin et al., 2009). In turn, the experimental outcomes based on the UCM hypothesis and optimal control model provide evidence supporting the current claim, both of which consistently suggest the minimization of a fraction of error variance as the controller’s strategies to govern the multi-element system. In this regard, the current findings show that the changes in the concurrent prefrontal cortex activation (i.e., functional connectivity) are affected by the frequency of cyclic actions, which is associated with the modulation of error variance (V_ORT), but not with the variance corresponding to the required mechanics (V_UCM), that is, the total finger forces.

To the best of our knowledge, the current study is the first work to uncover the effect of frequency on the multi-finger synergies in relation to the hemodynamic response in the prefrontal cortex, which is interconnected with diverse brain areas, including the cortical, subcortical, and brainstem regions (Miller and Cummings, 2017). In particular, the experimental findings in this study would contribute to the current knowledge on the neural origin of finger-force combinations for the stabilization of salient performance variables. However, this study has several limitations. First, the measurement of hemodynamics using fNIRS yields is more beneficial in terms of portability, cost-efficiency, and absence of hazardous effects of radiation. On the other hand, quantifying hemodynamic indices through fNIRS provides an indirect measurement of neuronal activity. Second, fNIRS may not be suitable for the direct monitoring of subcortical activity, even though the motor cortico-basal ganglia circuit originates from the prefrontal cortex (Fuster, 2015). As the formation of synergy for stable performance is assumed to be a function of subcortical structures (i.e., trans-thalamic loop; Rispal-Padel et al., 1981; Park et al., 2012b; Latash and Huang, 2015), future studies will have to ascertain the validity of fNIRS measures by comparing them with those of other brain imaging data during behavioral tasks. Also, the apparent drawback of the participant recruitment was that the gender of all participants was male with a small sample size. Although the results of previous studies claimed that the stability indices are not a gender-dependent (i.e., strength-dependent) quantity if the participants have no medical history of the peripheral and neurological disorder (Zhang et al., 2007; Friedman et al., 2009), the extremely swayed recruitment as to the gender may cause uncertainty and hamper the interpretation and generalization of the current messages. Lastly, given the diverse factors related to the coordinative action in humans, future studies will have to investigate whether the current message would be valid with age-related changes, gender, and other biological indices.