## Abstract

We can easily learn and perform a variety of movements that fundamentally require complex neuromuscular control. Many empirical findings have demonstrated that a wide range of complex muscle activation patterns could be well captured by the combination of a few functional modules, the so-called muscle synergies. Modularity represented by muscle synergies would simplify the control of a redundant neuromuscular system. However, how the reduction of neuromuscular redundancy through a modular controller contributes to sensorimotor learning remains unclear. To clarify such roles, we constructed a simple neural network model of the motor control system that included three intermediate layers representing neurons in the primary motor cortex, spinal interneurons organized into modules and motoneurons controlling upper-arm muscles. After a model learning period to generate the desired shoulder and/or elbow joint torques, we compared the adaptation to a novel rotational perturbation between modular and non-modular models. A series of simulations demonstrated that the modules reduced the effect of the bias in the distribution of muscle pulling directions, as well as in the distribution of torques associated with individual cortical neurons, which led to a more rapid adaptation to multi-directional force generation. These results suggest that modularity is crucial not only for reducing musculoskeletal redundancy but also for overcoming mechanical bias due to the musculoskeletal geometry allowing for faster adaptation to certain external environments.

## 1. Introduction

We can easily learn and perform desired movements, such as reaching for and grasping a glass or throwing a ball at a target. The central nervous system (CNS) uses an extraordinary number of variables to determine the proper solution for generating appropriate muscle activation patterns and for controlling movements [1]. Previous empirical findings have demonstrated that complex activation of a large number of muscles could be adequately captured by the combination of a few modularly coordinated muscle activation patterns, the so-called muscle synergies [2–6]. Muscle synergies have been estimated based on the linear decomposition of muscle activation patterns in various tasks [6–10]. Moreover, several findings have shown that similar muscle synergies were extracted across behaviours with different biomechanical demands, indicating that muscle synergies underlying a complex motor repertoire were shared across different motor tasks [11–14]. These studies demonstrated the robustness of modular architectures in muscle synergies for the generation of movements with reduced control parameters. Accordingly, modularity would reduce the number of parameters that the neuromuscular system must control [1]. When we learn to move in a novel environment, a modular controller will reduce the computational effort by reducing the dimensionality required to explore an appropriate motor solution. How, then, does the reduction of neuronal redundancy through a modular controller contribute to motor learning?

A previous study demonstrated that modular architectures in muscle synergies were reorganized only when the required muscle activation patterns were not represented by the combination of the existing muscle synergies during adaptation to a complicated environment [15]. This study provided important evidence that extracted muscle synergies might reflect the constraints of neuromuscular control. In addition, a previous simulation study showed that a modular controller could generalize to different but related tasks faster than a non-modular controller [16]. However, how a hard-wired modular controller contributes to adaptation to a novel environment remains unclear. To examine this question, it is necessary to compare the motor learning performed in the presence and absence of the modules. Accordingly, we performed a simulation using a neural network model with and without modules. Neural network models have been widely used to clarify neural properties [17–19]. In the present study, we constructed a descending neural network model with three intermediate layers representing neurons in the primary motor cortex (M1), spinal modules and motoneurons controlling individual muscles. Recent studies have provided evidence of modular controllers in the neuronal circuitry [15,20–23]. Based on these studies, we regarded the modules as interneuronal circuitry in the spinal cord [2,22,23]. The neurons of each layer have descending connections based on the empirical findings of the relationship between neurons in M1 and muscle synergies [24,25] and of the correlation between the activations of interneurons (regarded as the modules) and certain muscles [22,23,26]. It should be noted that the feed-forward neural network model in the present study is a large simplification of the real complexity of the motor control system, but evidence [17–19] suggests that such simplifications can provide meaningful insight, with predictive power, into the working of the real system.

Therefore, the purpose of the present study was to clarify the functional roles of modularity during motor learning based on a simulation of redundant motor control. To this end, a descending neural network model with modules was first constructed. Then, the speed when adapting to a visuomotor transformation [27,28] was compared in the presence and absence of the modules included in this model. This study provided the novel finding that modular control is important not only for reducing controlling variables but also for coordinating neuronal circuits. Furthermore, our findings are important in further testing the muscle synergy hypothesis, which is still controversial in the field of neuroscience [29].

## 2. Material and methods

### 2.1. Construction of the neural network model including modules

The descending neural network model including three intermediate layers (modular model; figure 1*a*) was constructed. The neural network performed the isometric torque production task with a two-joint system (shoulder and elbow) described in experimental [30] and simulation studies [18]. The elbow and shoulder joints were flexed to 90° from full extension and 30° from the frontal plane, respectively. The input layer was the desired shoulder flexion (or extension) and/or elbow flexion (or extension) torques, **τ***=* (**τ**_{1}, *τ*_{2})* ^{T}*. The 12 desired torques were uniformly distributed in 30° increments to cover the entire shoulder and elbow joint torque plane with the same intensity that was set to 1 Nm (figure 1

*b*). The first intermediate layer contained 1000 M1 neurons (

*N*

_{neu}= 1000). Each neuron (for example, the

**th neuron) received the desired shoulder flexion and/or elbow flexion torque vectors from the input layer with a synaptic weight (). Hence, the activation of each M1 neuron () obeys cosine tuning [17,18,31]. These M1 neurons were connected to an intermediate layer that indicated the modules and was regarded as a spinal interneuronal system [2,20,22,23]. Thus, these connections were considered the corticospinal pathway [32]. The number of modules () was determined according to procedures described below. The innervation of modules () by M1 neurons was uniformly distributed on the surface of a hypersphere in the dimensional space, the radius of which was 2/**

*i**N*

_{neu}. Each module controlled 18 muscles (

*N*

_{mus}= 18; deltoid anterior [DeltA], pectoralis major [PectMaj], coracobrachialis, subscapularis, supraspinatus, biceps brachii short head [BicShort], biceps brachii long head [BicLong], brachialis [Brac], brachioradialis [BracRad], deltoid posterior [DeltP], deltoid middle [DeltM], latissimus dorsi [LatDorsi], infraspinatus [InfraSp], teres major [TerMaj], teres minor [TerMin], triceps brachii long head [TriLong], triceps brachii lateral head [TriLat] and triceps brachii medial head [TriMed]) spanning the elbow and/or shoulder joints with a synaptic weight () in the third intermediate layer, i.e. the muscle layer. Accordingly, the activations of modules and muscles were determined by the function of the torque directions, indicating that they approximately obeyed cosine tuning [8,31,33]. The

*W*_{mod}was assumed to be positive. The positive constraint was applied in previous studies in which muscle synergies were estimated using a non-negative matrix factorization algorithm because muscles are pull-only actuators [2,6,34]. The activations of each module () and muscle () were determined by the following equations: 2.1and 2.2where the operator represents that for and for . The output layer was the resultant shoulder and/or elbow joint torques,

**= (**

*T*

*T*_{1},

*T*_{2})

*, generated by the combination of the mechanical directions (MDs) of muscles (see below).*

^{T}

*W*_{inp}and

*W*_{mod}were modified using an error back-propagation algorithm [35] to minimize the sum of the squared error () between output torque (

**) and the randomly selected 12 desired input torques (**

*T***): 2.3where**

*τ**t*represents the trial (

*t*= 1, 2, … ,

*N*

_{trial};

*N*

_{trial}: maximum number of trials).

*W*_{inp}and

*W*_{mod}were updated according to the following equations: 2.4and 2.5where

*α*

_{tor}and

*α*

_{mod}represent the learning rate of neurons and modules, respectively (

*α*

_{inp}=

*N*

_{neu}/50,

*α*

_{mod}=

*N*

_{mod}/8), and

*E*means the sum of squared error, . The third term is a decay term that prevents (or ) from drifting without bounds [18]. and represent the decay rate (). This decay term also minimized the sum of the squared activation of M1 neurons [18] as well as modules. As noted above, the

**was assumed to be positive. The operator was the same as that described above.**

*W*_{mod}The network was trained to produce the desired joint torques over 30 000 trials across the desired joint torques (*N*_{trial} = 30 000). In each trial, one of 12 equally distributed targets was randomly presented. This continuous learning was then repeated 30 times to change the initial value. Learning speed, *b*, was calculated by fitting the exponential function to [19]. The values for learning performance, such as learning speed and final error after learning, were averaged in the 30 sets of simulation.

### 2.2. Mechanical directions

We calculated the muscle MD (; figure 1*c*) based on the physiological parameters of each muscle used in a previous study [18], such as the cross-sectional area, *S _{m}*; pennation angle,

*α*[36]; muscle tension,

_{m}*γ*[37] (62 N cm

^{−2}); moment arm, [18]; and the following equation: 2.6

Each was then normalized with the average norm of the for all muscles. We also defined MDs of a neuron and a module as and , respectively.

### 2.3. Initial condition of the neural network model

The initial synaptic weights of *W*_{inp} and *W*_{mod} were set to random values as follows: and , where and represent white noise with a zero mean and indicates their standard deviation. The standard deviations and were set to 0.5 except for one verification, in which five different initial weighting matrices of M1 neurons were generated using a different standard deviation, i.e. 0.5, 1.0, 1.5, 2.0 and 2.5, to evaluate the state of the neural network across the different initial synaptic weights.

### 2.4. Determination of the number of modules

The best number of modules () was determined according to Akaike's information criterion (AIC) [38] and the Bayesian information criterion (BIC) [39]. To this end, the number of modules (*N*_{mod}) was varied from 1 to 30 across continuous learning. The sum of the squared error between the input and output torques was calculated for all desired torques across a different number of modules (*N*_{mod} = 1, 2, … , *N*_{mus}). Based on the errors, the values of the AIC and the BIC were calculated to measure the fitness between the model and the desired torques. Then, we selected the number of modules at which the values of the AIC and the BIC were lowest. The model based on this number of modules adequately fitted the desired torques and minimized redundancy in the modular controller.

### 2.5. Comparison of the learning performance between models with and without modules

Learning performance in the presence and absence of the modularity was evaluated by comparing the model with modules to the previously constructed model, which did not include the layer of modules (figure 1*d*; non-modular model [18]), when the two sets of the neural network model adapted to the change in an external environment. To this end, the tasks were simulated by assuming adaptation to a rotational perturbation, i.e. a force field or a visuomotor rotation. When the state of the neural network was approximately constant after learning (*t* = 30 000), the output torque was exposed to a rotational perturbation , where *ϕ* represents the rotational angle (*ϕ* = 30°). Hence, the resultant output torque, , was determined using the following equation: . During the adaptation to the rotational perturbation, the synaptic weight, *W*_{inp}, was modified in both neural network models [40] to minimize the error between the resultant output torque (*T*_{final}) and 12 randomly selected desired input torques (** τ**). In a particular simulation, however, we tested the effect of modifying both

*W*_{inp}and

*W*_{mod}. The adaptation trial of the rotational perturbation was conducted for 5000 trials. To equally modify the mechanical contribution of M1 neurons across all trials in both models, the innervation weight in the non-modular model, , was compensated for. At first, we calculated the sums of the mechanical contribution of neurons (neuron MDs), i.e. , after initial training in the modular model. Then, we divided by . Depending on the solution, was modified so that the sums of the magnitude of neuron MDs were equal to each model, i.e. .

### 2.6. Comparison of the learning performance between models with uniformly distributed and biased muscle mechanical directions

The learning performance of the rotational perturbation was compared between the modular and non-modular models with equal numbers of dimensions controlling the muscles. To this end, we constructed the models with two types of muscle MDs, ** M_{uniform4}** and

**, which were uniformly distributed and biased to the first and third quadrants on the torque plane, respectively.**

*M*_{bias4}**was composed of the following four-unit vectors: 2.7**

*M*_{uniform4}** M_{bias4}** was calculated as

**, where ().**

*R*_{bias}*M*_{uniform4}### 2.7. Non-modular model with lower bias of distribution in 18-muscle mechanical directions than original muscle mechanical directions

To test whether the bias of distribution in muscle MDs affects the learning performance against the rotational perturbation, we performed the simulation using the model with 18-muscle MDs the bias of the distribution of which was lower than the original, i.e. as shown in figure 1*c*. Various patterns of the low-biased muscle MDs were produced by randomly rotating each muscle MD. Then, we grouped them in two groups; one was the bias of muscle MD distribution which was lower than the original and higher than the module-MD distribution, and the other was the bias of muscle MD distribution which was lower than the module-MD distribution in the intrinsic condition with 18 muscles. In each group, 30 patterns of the low-biased muscle MDs were constructed.

### 2.8. Influence of different numbers of initial trainings on learning performance against rotational perturbation

We verified whether the different number of initial training trials affected the learning performance in the adaptation period. To this end, the neural network was trained with three different numbers of training trials, i.e. *t* = 5000, 10 000 and 20 000, in addition to the original, *t* = 30 000. Then, we compared the learning performance against the rotational perturbation between models with and without modules. The adaptation trial of the rotational perturbation was conducted for 5000 trials.

### 2.9. Calculation of preferred direction

To quantify the directional dependence of neuron activations on the torque plane, multiple linear regression was used to fit the following equation to the activations of each neuron, , after learning:
2.8where *θ* is the direction on the torque plane, *c*(*θ*) is the activation of neurons across each direction and is an offset parameter [8]. The amplitude (*r*) and preferred direction (PD; *θ*^{PD}) of neuron activation were then estimated using the parameters and as follows:
2.9and
2.10

The goodness-of-fit was quantified with an *R*^{2} value, and its significance was tested with an *F* test.

### 2.10. Analysis of preferred direction and mechanical direction distribution

The significance of the bimodal distribution in neuron PDs was tested based on the Rayleigh test, which detects unimodal deviation from uniformity [41]. The PDs () were first multiplied by 2 and then transformed to unit vectors, i.e. (), so that the Rayleigh test was used for bimodal distributions where the two modes were assumed to be oriented at 180° from each other [42]. The direction () and length (*R*) of the resultant vector, which represented the direction and strength of the PD bias, respectively, were calculated as follows:
2.11and
2.12where *N*_{neu} is the number of neurons for 30 sets of simulation. To compare them with experimental data, we randomly extracted the *N*_{neu} raw PDs from the literature [30] (178 neurons for an isometric task from figure 9a; the duplication of neurons was allowed) in each resampling. The 95% confidence interval was then estimated using a bootstrapping procedure with 10 000 times resampling [18].

The distribution of neuron MDs was quantified by the angle of the longitudinal axis (). To estimate the bias in the distribution of neuron, module and muscle MDs, eigen values, *σ*_{1} and *σ*_{2}, were calculated from the variance–covariance matrix of each MD distribution. Then, we defined the bias of the MD distribution as (*σ*_{1} > *σ*_{2}). These parameters were calculated across each of the 30 simulations.

### 2.11. Simulation in the two-dimensional extrinsic space model

We simulated the initial phase of goal-directed reaching movements on a two-dimensional extrinsic task plane. The model produced a linear acceleration of the fingertip on a horizontal plane (*xy* coordinate; electronic supplementary material, figure S5) by replacing the input and output torques with a desired and produced linear acceleration of the fingertip, respectively. A two-dimensional forward-kinematics model for a two-joint system (shoulder and elbow) was given by
2.13where and are the shoulder and elbow joint flexion/extension angles, respectively. In addition, *l*_{1} and *l*_{2} represent the length of the proximal (an upper arm) and distal (a forearm and hand) links, respectively. Based on the forward-kinematics model, the Jacobian matrix () was computed as:
2.14

The muscle MDs in linear acceleration space () were calculated as the following equation:
2.15where is the inertia matrix of the two-joint system. All morphological data were based on the parameters used in the previous study [18]. The model based on four modules was selected, which adequately fitted the desired linear acceleration of the fingertip and minimized redundancy in the modular controller. The learning rate of neurons (*α*_{neu}) and modules (*α*_{mod}) were set to *N*_{neu}/500 and *N*_{mod}/80 for the initial training period. The network learned over 50 000 trials across the desired fingertip accelerations. After the model learned to produce the desired fingertip acceleration, the output was exposed to a rotational perturbation, ** R**. During the adaptation to the rotational perturbation, the synaptic weight,

*W*_{inp}, was modified over 5000 trials. Then, we compared learning performance between the models with and without modules. In the rotational period, the learning rate of neurons (

*α*

_{neu}) was set to

*N*

_{neu}/50.

### 2.12. Simulation in three-dimensional extrinsic space model

We also simulated the initial phase of goal-directed reaching movements to each of 33 targets distributed on the surface of a sphere in a three-dimensional extrinsic space with the 4 degrees-of-freedom joint system (3 degrees of freedom for the shoulder joint and 1 degree-of-freedom for the elbow joint; electronic supplementary material, figure S7A,B). We set the shoulder flexion angle at 30°, the internally rotated shoulder angle at 12°, and the elbow flexion angle at 80° [18]. The model produced a linear acceleration of the fingertip in three-dimensional space (*xyz* coordinate; electronic supplementary material, figure S7A) by replacing the input and output torques with a desired and produced linear acceleration of the fingertip, respectively. The Jacobian matrix () was computed as:
2.16where , , and are the shoulder joint flexion/extension, internal/external rotation, elevation/depression and elbow joint flexion/extension angles, respectively. The muscle MDs in linear acceleration space () were calculated using the following equation:
2.17where is the inertia matrix of the 4 degrees-of-freedom joint system. The inertia matrix, , was computed using previously described equations [43]. All morphological data were calculated based on the musculoskeletal parameters used in the previous study [18]. The model based on six modules was selected, which adequately fitted the desired linear acceleration of the fingertip and minimized redundancy in the modular controller. The learning rate of neurons (*α*_{neu}) and modules (*α*_{mod}) were set to *N*_{neu}/50 and *N*_{mod}/80. The network learned over 50 000 trials across the desired fingertip accelerations. After the model learned to produce the desired fingertip acceleration, the output was exposed to a rotational perturbation as follows:
2.18where *ϕ* represents the rotational angle (*ϕ* = 30°). During the adaptation to the rotational perturbation, the synaptic weight, *W*_{inp}, was modified over 10 000 trials. Then, we compared learning performance between the models with and without modules.

### 2.13. Similarity of muscle-weighting vectors from modules across different models

The similarity of muscle-weighting vectors from the module to muscle layer was calculated across the different models based on the inner product of the vectors [6]. Because each weighting vector was normalized by its norm, the inner product of the two vectors represented the cosine of the angle between the vectors. Thus, the inner product closer to 1 indicated a greater similarity in the directions of the two vectors. We compared the all pairs of vectors between the models. To verify the significance of the similarity, we also calculated the inner product across each pair of muscle-weighting vectors in which the order of muscles was shuffled. The confidence interval was estimated using a bootstrapping procedure with 1000 times resampling across each pair of vectors.

### 2.14. Models with different distributions of neuron mechanical directions

To compare the learning performance across different distributions of the neuron MD, a simple descending neural network model including one intermediate layer as the M1 layer, i.e. a non-modular model without a muscle layer, was considered. In this model, the neuron MDs were biased to the first and third quadrants with different degrees of the rotational bias, i.e. (). was composed of eight-unit vectors uniformly distributed on the torque plane; i.e. 2.19

) was uniformly distributed on the surface of a hypersphere in the eight-dimensional space, the radius of which was 2/*N*_{neu}.

### 2.15. Statistical tests

For the statistical tests that compared the learning performance and the parameters indicating the distribution of M1 neurons between two different types of models, such as modular and non-modular models, we used the Wilcoxon signed-rank test. The effect of the models based on the different number of modules to the learning performance was tested using Friedman's test. Friedman's test was also performed to test the difference of learning performance in modular and non-modular models with the different number of initial training trials.

## 3. Results

### 3.1. Neural network model with modular controllers

A descending neural network model learned to produce the desired joint torques by modifying the synaptic weights from the input to neuron layers and from the module to muscle layers to minimize the sum of squared error between the input and output torques. According to the AIC and the BIC, we selected the model based on four modules () that best fitted the desired torques in all calculated models (see Material and methods and electronic supplementary material, figure S1). During learning of the neural network model, the error gradually decreased (electronic supplementary material, figure S2). After the learning period (*t* = 30 000), the preferred direction (PD) of neurons was distributed along a bimodal axis (; figure 2*a,b*) as observed in an experimental study [30]. Conversely, the MD of M1 neurons denoted the larger contribution of the first and third quadrants (; figure 2*c*). In the muscle layer, a misalignment of the PD and MD of muscles was observed, as shown in a previous experimental study [44]. Although the MD of each muscle was directed to shoulder flexion and/or elbow flexion or to shoulder extension and/or elbow extension (figure 1*c*), the PD of mono-articular muscles was rotated to the shoulder extension and elbow flexion directions or shoulder flexion and elbow extension directions on the torque plane (figure 2*d*).

In contrast, both the PDs and MDs of the modules were distributed to all quadrants on the torque plane (figure 3). The module **W**_{1}, which consisted of mono-articular shoulder flexors, bi-articular flexors and mono-articular elbow flexors, was activated around the direction of the shoulder and elbow flexion (*θ*_{PD} = 76.8°).*θ*_{PD} = The MD of module **W**_{1} was directed to the same quadrant on the torque plane (*θ*_{MD} = 60.9°). The module **W**_{2} constructed by bi-articular flexors, mono-articular elbow flexors, mono-articular shoulder extensors and the bi-articular extensor was activated on the second quadrant on the torque plane (*θ*_{PD} = 150.6°), which produced shoulder extension and elbow flexion torques (*θ*_{MD} = 165.3°). The module **W**_{3} was dominant for bi-articular and mono-articular elbow extensors, the PD and MD of which were in the same quadrant relevant to shoulder and elbow extension on the torque plane (*θ*_{PD} = 254.6°; *θ*_{MD} = 240.0°). The module **W**_{4}, which was mainly composed of mono-articular shoulder flexors and elbow extensors, was activated in the fourth quadrant on the torque plane (*θ*_{PD} = 330.5°) to generate shoulder flexion and elbow extension torques (*θ*_{MD} = 346.9°). These results suggest that the bias of the muscle MD due to the musculoskeletal configuration was compensated by modularity as the neuronal constraint.

### 3.2. Comparison between models with and without a modular controller

The main purpose of this study was to reveal the functional roles of modularity during motor learning. To this end, we first compared the neural state between models with and without a modular controller after an initial training period (*t* = 30 000). There was no difference between the bimodal axes of neuron PD distributions (*p* = 0.2536; 130.2° [s.d. 2.63] and 131.0° [s.d. 7.30] in modular and non-modular models, respectively), whereas the strength of the PD bias (*R*) was greater in the modular model than in the non-modular model (*p* = 6.16 × 10^{−4}; *R* *=* 0.167 [s.d. 0.02] and 0.208 [s.d. 0.05] in modular and non-modular models, respectively). The value of *R* in both models was within the 95% confidence interval (0.149–0.330) based on the experimental data [30]. In the case of the neuron MD distribution, the angles of the bimodal axes were different (*p* = 0.0036; [s.d. 1.98] and 41.8 [s.d. 2.06] in models with and without modules, respectively), although these were in the same quadrant. The bias of the neuron MD distribution was higher in the non-modular model than in the modular model (*p* = 1.73 × 10^{−6}; 1.96 [s.d. 0.11] and 2.42 [s.d. 0.17] in modular and non-modular models, respectively). We then compared the learning performance against the rotational perturbation between the models with and without the modules. In the rotational period, only the synaptic weights of M1 neurons () were updated, i.e. adapting the sensorimotor transformation [27,28], because an empirical study demonstrated that the modular structure represented by muscle synergies did not change during adaptation to a visuomotor rotation [40]. It is noted that the variance of the distribution of the neuron MDs in the non-modular model was comparable to that in the modular model (see ‘*Comparison of the learning performance between models with and without modules*' in Material and methods): the total changes of the output torques were equal for both models when the same amount of was updated. The errors between the desired and produced torques in the models with and without modules are shown in figure 4*a*. Learning speed was calculated by fitting the exponential function to the error curve across each trial. The learning speed averaged across all targets was increased in the modular model compared with the non-modular model (figure 4*a,b* and table 1; *p* = 1.73 × 10^{−6}). In addition, the error after learning (*t* = 5000) was smaller in the modular model than in the non-modular model (figure 4*c*; table 1; *p* = 1.73 × 10^{−6}). When we compared the learning performance across each target, the modular model learned faster for six out of 12 targets (figure 2*d*; *p* = 0.098, 0.206, 0.299, 0.018, 5.75 × 10^{−6}, 0.262, 0.254, 1.80 × 10^{−5}, 1.92 × 10^{−6}, 1.97 × 10^{−5}, 0.262, 0.035 for targets 1−12 in figure 1*b*, respectively) and led to the smaller errors for all targets compared with the non-modular model (figure 2*e*; *p* = 1.73 × 10^{−6} for targets 1−12 except for target 5, and 1.92 × 10^{−6} for target 5 in figure 1*b*). In the beginning of the adaptation period, the total activation of neurons and muscles was increased from the initial activation level before the adaptation in both models (figure 4*d,e*). Subsequently, the activations immediately converged to being constant in the modular model (*t* = 100), whereas they were modified for a long period and became constant in the non-modular model. The total activation of both neurons and muscles was lower in the modular model (figure 4*d,e*). The total activation of modules showed a pattern similar to the activation of neurons and muscles in the modular model (figure 4*f*). During adaptation to the rotational perturbation, the PDs of the modules were rotated in the opposite direction to that of the perturbation due to the reformation of the neuron weights (figure 4*g*).

In order to verify whether the assumption that only was updated during the adaptation to the rotational perturbation affected the results, we compared the two models in which the synaptic weights of modules () were updated or not in addition to (table 1; electronic supplementary material, figure S3A,B). There was a difference in the learning speed and the error after learning (*t* = 5000) between the two models (*p* = 1.73 × 10^{−6} and *p* = 1.73 × 10^{−6}, respectively). Importantly, corresponding to the previous result noted above [40], was similar between the two models even though both the and were updated in one model (electronic supplementary material, figures S3D,E; the values of similarity were 0.97, 0.98, 0.97 and 0.98 in *W*_{1–4}, respectively; *p* < 0.01).

### 3.3. Influence of the number of initial training trials on learning performance against rotational perturbation

To verify whether the duration of the initial training period affected the learning performance in the adaptation period, we trained the neural network with three different numbers of training trials in addition to the original number of training trials, *t* = 30 000. Then, we compared the learning performance against the rotational perturbation between models with and without modules. As observed above, learning speed remained faster in the modular model than in the non-modular model irrespective of the duration of training (*p* = 1.40034 × 10^{−31}; table 1; electronic supplementary material, figure S4). This result indicates that the improvement of learning speed in the modular model was not affected by the network state due to the number of training trials after the errors had converged. Furthermore, errors after learning were smaller in the modular model than in the non-modular model (*p* = 8.90867 × 10^{−38}).

### 3.4. Simulation in two-dimensional extrinsic space model

We simulated the initial phase of goal-directed reaching movements in a two-dimensional extrinsic task plane (electronic supplementary material, figures S5 and S6). The model produced a linear acceleration of the fingertip to each of 12 directions on a horizontal plane. We compared the neural state between models with and without a modular controller after an initial training period (*t* = 50 000). There was no difference between the bimodal axes of neuron PD distributions (*p* = 0.1204; 126.9° [s.d. 1.08] and 126.4° [s.d. 1.25] in modular and non-modular models, respectively), whereas the strength of the PD bias (*R*) was greater in the non-modular model than in the modular model (*p* = 1.73 × 10^{−6}; *R* *=* 0.376 [s.d. 0.013] and 0.577 [s.d. 0.020] in modular and non-modular models, respectively). In the case of the neuron MD distribution, the angles of the bimodal axes were different (*p* = 0.0053; = 37.3 [s.d. 0.50] and 37.8 [s.d. 0.54] in models with and without modules, respectively). The bias of the neuron MD distribution was higher in the non-modular model than in the modular model (*p* = 1.73 × 10^{−6}; 5.19 [s.d. 0.18] and 20.91 [s.d. 1.26] in modular and non-modular models, respectively). All muscle-weighting vectors from the module to muscle layer were similar to the vectors in the intrinsic model (figure 3*a*), in which the representation of the number, i.e. **W**_{1–4}, was corresponding to that in the extrinsic model (electronic supplementary material, figure S6A; *p* < 0.05; the values of similarity were 0.82, 0.99, 0.81 and 0.99 in **W**_{1–4}, respectively). We then compared the learning performance against the rotational perturbation between the models with and without the modules. The learning speed was higher for all targets in the modular model compared with the non-modular model (*p* = 1.73 × 10^{−6}; table 1; electronic supplementary material, figure S6D,E). In addition, the error after learning (*t* = 5000) was smaller in the modular model than in the non-modular model (*p* = 1.73 × 10^{−6}; table 1; electronic supplementary material, figure S6D,F).

### 3.5. Simulation in the three-dimensional extrinsic space model

We also performed the simulation for the initial phase of goal-directed reaching movements to each of 33 targets distributed on the surface of a sphere in the three-dimensional extrinsic space with a 4 degrees-of-freedom joint system (3 degrees of freedom for the shoulder joint and 1 degree-of-freedom for the elbow joint; electronic supplementary material, figures S7 and S8). Muscle-weighting vectors, , and (electronic supplementary material, figure S8A), from the module to muscle layer were similar to the vectors, ** W_{1}**,

**and**

*W*_{2}**, in the intrinsic model, respectively (figure 3**

*W*_{3}*a*;

*p*< 0.05; the values of similarity were 0.79, 0.90 and 0.72 in , respectively). The learning performance during adaptation to the rotational perturbation was compared between the models with and without the modules. The learning speed was increased in the modular model compared with the non-modular model (

*p*= 5.6061 × 10

^{−6}; table 1; electronic supplementary material, figure S8C,D). In addition, the error after learning (

*t*= 10 000) was smaller in the modular model than in the non-modular model (

*p*= 1.3706 × 10

^{−6}; table 1; electronic supplementary material, figure S8C,E).

### 3.6. Functional roles of modularity in motor learning

We next investigated why the learning speed increased and the total task error decreased in the model assuming a modular controller. During the rotational period, because only the synaptic weight, ** W_{inp}**, was updated, the contribution of individual M1 neurons to the output torques would be an important factor that affected the learning performance. Because the dimensions of the mechanical contributions were different between the models with and without modules, i.e. they converged to four and 18 dimensions, respectively, we first tested the effect of the dimensionality reduction due to a modular controller using the models with four muscle MDs,

**(‘uniform’ condition) and**

*M*_{uniform4}**(‘bias’ condition; see Material and methods,**

*M*_{bias4}*Comparison of the learning performance between models with uniformly distributed and biased muscle MDs*). As a result of the simulation to learn to produce the desired torques against the rotational perturbation, the difference in the learning speed and final error between modular and non-modular models was substantially larger for the ‘bias’ condition than for the ‘uniform’ condition although even for the ‘uniform’ condition the reduction of the subtle bias in the initial M1 distribution due to the modules led to the higher learning speed and the smaller final error in the modular model than in the non-modular model (

*p*= 1.73 × 10

^{−6},

*p*= 1.73 × 10

^{−6}, respectively; table 1; figure 5). This result suggested that the bias in the muscle layer largely affected the learning performance. Indeed, the model with 18 low-biased muscle MDs learned faster against the rotational perturbation than the model with the 18 original muscle MDs (

*p*= 0.0041; table 1; see Material and methods,

*Non-modular model with lower bias of distribution in 18-muscle MDs than original muscle MDs*). If the bias of the muscle MD distribution was lower than the bias of the module-MD distribution (figure 3

*c*), i.e. 1.43, there was no difference in the learning speed between the modular model and non-modular model with low-biased muscle MDs (

*p*= 0.0643; table 1). Furthermore, the neuron MDs reflected the bias of muscle MDs: the distribution of the neuron MDs showed a steeper bias in the non-modular model than in the modular model (

*p*= 1.73 × 10

^{−6}; 1.96 [s.d. 0.11] and 2.42 [s.d. 0.17] in modular and non-modular models, respectively) as was the case for the model including 18 muscles. Therefore, the learning performance would relate to the bias of the neuron MD distribution rather than the dimensionality reduction due to the modular controller.

Additionally, we compared the learning speed and final error across the model based on the different numbers of modules during the rotational condition (electronic supplementary material, figure S9). There was no difference in the learning speed, whereas the final error tended to increase with the number of modules. These results indicated that the degree of dimensionality reduction due to the modular controller did not affect the learning speed during adaptation to the rotational perturbation in this study.

### 3.7. Relationship between learning performance and bias in neuron MDs

According to our simulation, the bias in the neuron MDs may be the crucial factor that affected the learning performance, including the learning speed and task error. To verify this prediction, we compared the learning performance of the neural network model among the four different degrees of bias in neuron MDs (; figure 6*a*). We observed that the learning speed decreased as the bias and the task error increased (figure 6*b,d,e*; table 1). In addition, both before and after adaptation to the rotational perturbation (*t* = 5000), the desired joint torques were produced with a low total neuron activation when the distribution of neuron MDs approached being uniform (figure 6*c*). Furthermore, we compared the learning speed of these four models across each target torque. As the bias of neuron MDs increased, the learning speed in the first and third quadrants on the torque plane decreased. However, the learning speed increased towards the directions of the second and fourth quadrants (150° and 330°, respectively; figure 6*b*), which corresponded to approximately 20° anticlockwise from the bimodal axis of the neuron PD distribution (*θ*^{bias} = 130.1° and 310.1°).

## 4. Discussion

### 4.1. Modularity maximizes the learning speed

Rearrangement of the neuronal circuitry is an important property of the nervous system that allows it to adapt to the external environment [45,46]. In the general scheme during sensorimotor adaptation, such as rotational perturbation simulated in this study, muscle activation is corrected based on an internal model encoded in the nervous system, which reflects the relationships between musculoskeletal properties and external environments [27,28,47]. In an isometric condition, the endpoint forces are directly produced by the combination of MDs of relevant muscles [48]. In the muscle synergy hypothesis, however, muscles are grouped into clusters whose MDs in the involved joints have a lower dimensionality than the muscles [49]. Therefore, the internal model would calculate the mechanical contribution of each modular controller [4] and the desired endpoint force would be determined by the combination of MDs in the modules [46,47]. The intrinsic mechanical property of the muscle MD was biased mainly due to the presence of bi-articular muscles. For example, there are no bi-articular muscles spanning elbow extension and shoulder flexion, and furthermore, the mechanical contributions of relevant muscles are different across each joint (figure 1*c*). Then, the state of M1 neurons fundamentally reflects this high degree of bias in the musculoskeletal system. However, our simulation demonstrated that the modularity reduced the influence due to the mechanical bias by constructing MDs of modules in all quadrants on the torque planes, which led to the reduction of the bias in the distribution of neuron MDs. Therefore, modularity would play an important role not only in reducing neuronal redundancy but also in reducing the bias of the muscle MDs. This reduction of the bias leads to an improvement in sensorimotor adaptation to multi-directional force generation (figure 4) because the internal model was modified based on the arranged MDs of modules rather than the biased muscle MDs. The bias reduction of the neuron and muscle MDs due to modularity also decreased the bias of the distribution in PDs, which were coordinated perpendicular to the distribution of the MDs. Consequently, the motor effort cost for the targets in the second and fourth quadrants on the torque plane (figure 1*c*) where no bi-articular muscles were distributed was reduced. This contribution of the modularity led to lower activation of neurons and muscles as observed in the modular model than in the non-modular model (figure 4*d,e*).

Using a closed-loop intracortical brain–computer interface learning paradigm in rhesus macaques, a previous study found that easy-to-learn and difficult-to-learn motor tasks arose from neuronal constraints in the cortical circuitry [50]. Moreover, an empirical study showed that the difficulty of motor learning also depended on muscle weights in muscle synergies [15], indicating that muscle synergies reflect constraints of the neuromuscular control. Modularity limits motor learning when the desired movements could not be reproduced by the combination of existing low-dimensional neuromuscular controllers. In addition, we demonstrated that the presence of modules in the neural network model overcame the inherent bias of easy-to-learn and difficult-to-learn movement directions. Therefore, the neuromuscular constraints contribute to reducing the mechanical bias due to the musculoskeletal geometry and speeding up sensorimotor learning in task space compatible with the existing modular controllers.

In multi-dimensional arm movements in extrinsic task space, the distribution of the neuron MDs was also found to be biased (electronic supplementary material, figures S5–S8), likely due to the bias in the musculoskeletal system [18,42]; however, this proposal is controversial [51]. Hence, the contribution of the modularity observed in the present study plays an important role during sensorimotor adaptation in these conditions.

As comparing the learning speed and errors after learning between models with and without module layers, the depth, i.e. three and two intermediate layers in modular and non-modular models, respectively, and dimensionality of the neural network could be important factors that gained the learning ability of the models. The deeper networks in the model with modules, however, did not make a significant improvement in the learning performance if the MDs of muscles were uniformly distributed (figure 5*a,c*). In addition, the results also demonstrated that the learning performance was higher in the modular model than in the non-modular model even if the dimensionality of the models was comparable to each other, i.e. four dimensions (figure 5*b,d*). Therefore, the modularity that reduces the bias in the musculoskeletal geometry was the crucial factor to enhance the learning capacity of the model rather than the depth and the dimensionality of the network.

### 4.2. Functional roles of modularity

Understanding the relationship between cortical neurons and extracted muscle synergies is important to clarify how cortical neurons modulate muscle activation in a low dimension [24,25]. We examined the relationship between neuron PDs and the synaptic weights (** W_{neu}**) from the neurons to each module (figure 7). Although the much larger number of neurons than modules appeared to lead to a redundancy in the motor execution, the neurons with similar PDs connected to the same modules. Therefore, the building blocks in the cortical layer may be related to the corresponding modules, which reduce a high degree-of-freedom in a high-level system. Neurons with similar PDs also showed a smaller variance of their MDs in the modular model than in the non-modular model (figure 8), indicating that the population of neurons with similar PDs was preferentially associated with the production of the similar joint torque due to the constraint of neuromuscular control. These results suggest that the modularity contributes to arranging the neuronal circuitry not only to reduce bias in the musculoskeletal system but also to explicitly define the building blocks in a cortical region.

The number of modules is an important parameter to determine the dimension of muscle activation that can be controlled. In the present study, we selected the model based on four modules that best fitted the target torques in all calculated models with the different number of modules according to the AIC and the BIC. However, even if the number of modules was more than 4, the contribution of the modularity to speeding up adaptation to the rotational perturbation was not affected (electronic supplementary material, figure S9). Therefore, in the present task space, four modules were sufficient to reduce the bias due to the distribution of the muscle MDs.

In the physiological system, minimal motor effort is often required to optimize the motor commands to precisely produce behaviour [52,53]. Previous computational studies using modular controllers required minimization of the motor effort cost to yield motor commands [54–56]. In our study, the motor effort cost was theoretically minimized by the decay term [18]. In the experimental case, however, several muscle synergies varied across individuals [12,48], implying that muscle synergies reflect not only an optimal solution but also a good-enough solution due to motor history, such as experience and training [57,58]. This study focused on the motor adaptation of cortical neurons after constructing the modular controller, whereas it is important to clarify based on what criteria the modular architectures were constructed and the structural roles of extracted muscle synergies.

According to previous simulations during balance control, the control effort was increased in the model with modules compared to individual muscle control, although the experimentally observed force was well reconstructed [56]. By contrast, our simulation demonstrated that the bias reduction of neuron MDs due to the modular controller contributed to a reduced total motor effort cost of neuron and muscle activity during multi-directional force generation (figures 4*d,e* and 6*c*). These results suggest that the contribution of the modularity to motor effort depends on the task contexts [56,59]. In the precise control of the limb such as the present tasks, reducing motor efforts would be needed, which leads to decreasing the endpoint variability [60], although strict minimal-effort solutions would not be necessary for discrete and instantaneous tasks such as postural responses against external perturbations [56].

### 4.3. Modularity in the neuronal circuit

Whether the modular controller is of neural origin is still an open question [29]. Indeed, muscle synergies extracted from the experimental electromyogram (EMG) data were also seemingly accounted for by the biomechanical constraints [61] or the product of optimization [62]. However, evidence that the neuronal representation of muscle synergies, i.e. the modular structure from a last-order spinal interneuron to motor neurons for multiple muscles, was provided in physiological studies [22,23,26]. Accordingly, in this study, we examined the functional roles of the modularity as the constraint of neuromuscular control. Although to clarify the relationship between the muscle synergies extracted from EMG data based on the decomposing technique and the physiologically observed neuronal structure of the modularity [22,23,26] was beyond the scope of this paper, we demonstrated that modular controllers that would be represented by muscle synergies play various important roles in the coordination in the nervous system to generate desired motor behaviours.

Previous findings have shown that muscle synergies were shared across behaviours with different biomechanical demands [11–14]. Shared modules were also demonstrated in the present study; the models that produced force in the intrinsic, two-dimensional extrinsic and three-dimensional extrinsic task spaces, were constructed based on the similar modules despite the separate training across each model (figure 3, electronic supplementary material, figures S6 and S8). Furthermore, the adaptation of the model to the novel rotational perturbation did not change the modules (electronic supplementary material, figure S3), which was observed in the previous experiment [40]. These results suggest that the modules in this study reflect the stable neuronal system rather than plastic. Robust modules might be constructed over evolution or long-term development, whereas adaptable modules would be necessary to learn unusual and highly skilled behaviours, such as novel tool use [63]. Future studies will test the stability and the plasticity of modules based on the model.

The existence of modular controllers as the interneuronal pathway in the spinal cord was previously identified in spinalized frogs using direct stimulation to the spinal circuitry [64,65] and the recording of muscle activity induced by perturbing proprioceptive feedback [2,66]. Hart and Giszter [20] demonstrated that the activity recorded from frog spinal interneurons was correlated with the activity of statistically estimated muscle synergies rather than individual muscles [20]. Other researchers found the neural basis of muscle synergies as last-order interneurons in the spinal cord using retrograde trans-synaptic rabies virus injected into mice muscles [22] and spike-triggered averaging of hand muscle activation in primates [23]. Furthermore, the correlation between extracted muscle synergies and cortical neurons was observed during both intracortical electrical stimulation of cortical regions and voluntary reaching and grasping movements [25,67], indicating that the modularity may originate in the cortex in some contexts. Based on the physiological findings, our neural network model, which included a modular controller at the spinal level [2,20,22,23], was able to reconstruct the neural state observed in the experiment [30]. This result supported the possibility that the modular controller exists in the descending pathway between cortical neurons and motoneurons. However, the existence of the modularity in the spinal cord was no more than an assumption. Many studies also reported monosynaptic connections with motoneurons from different regions, such as reticulospinal and vestibulospinal pathways [68] or the cortico-motoneuronal (CM) pathway [69] that do not include spinal interneurons, i.e. the spinal modules in the present model. Indeed, physiological and anatomical findings implied that the modular controllers might be located in the spinal cord and the medulla of a frog [70] and also in cortical circuits of a primate [71,72]. In addition to the layer of modules, these neuronal pathways were driven depending on the task context: the reticulospinal pathway contributes to coordinating trunk and proximal muscles and inducing locomotor activity [73], whereas CM neurons predominantly project to the distal hand muscles [74] and are important for generating and controlling highly skilled movements [75]. Although the simulated tasks in this model may arise primarily from the modular controller encoded in the spinal circuit, the actual biological system consists of various neuronal pathways that would include a hybrid of the modular and non-modular model architectures [76]. It is important for the identification of muscle synergies in the neuronal system to consider other important neuronal pathways and spinal characteristics [77] in future studies.

In summary, the neural network model that included modules showed an increased learning speed compared with the model without modules during multi-directional torque production tasks. This improvement in the learning performance in the model with modules was attributed to a reduction of the bias in the musculoskeletal geometry, which reflected the coordination of neuron MDs. These results suggest that the functional roles of modular controllers represented by muscle synergies are not only to simplify the redundancy in the biological system but also to arrange the network so that the CNS rapidly and efficiently adapts to various environments.

## Data accessibility

The Matlab codes for running the simulations are deposited in the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.3h70d78 [78].

## Authors' contributions

Conception and design of the experiments: S.H. and M.K. Collection, analysis and interpretation of data: S.H. Drafting of the article and revising it critically for important intellectual content: S.H. and M.K. Final approval of the version to be published: S.H. and M.K.

## Competing interests

The authors declare no competing interests.

## Funding

This work was supported by JSPS KAKENHI grant number 13J00233.

## Acknowledgements

We are grateful to Dr Masaya Hirashima (Center for Information and Neural Networks) for his technical comments. We also thank Dr Keisuke Fuji (Institute of Physical and Chemical Research [RIKEN]) for his helpful discussions.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4235432.

- Received April 12, 2018.
- Accepted September 5, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.