Abstract
How does the motor cortex (MC) produce purposeful and generalizable movements from the complex musculoskeletal system in a dynamic environment? To elucidate the underlying neural dynamics, we use a goal-driven approach to model MC by considering its goal as a controller driving the musculoskeletal system through desired states to achieve movement. Specifically, we formulate the MC as a recurrent neural network (RNN) controller producing muscle commands while receiving sensory feedback from biologically accurate musculoskeletal models. Given this real-time simulated feedback implemented in advanced physics simulation engines, we use deep reinforcement learning to train the RNN to achieve desired movements under specified neural and musculoskeletal constraints. Activity of the trained model can accurately decode experimentally recorded neural population dynamics and single-unit MC activity, while generalizing well to testing conditions significantly different from training. Simultaneous goal- and data-driven modeling in which we use the recorded neural activity as observed states of the MC further enhances direct and generalizable single-unit decoding. Finally, we show that this framework elucidates computational principles of how neural dynamics enable flexible control of movement and make this framework easy-to-use for future experiments.
1 Introduction
Behavior arises from the interaction among functionally and anatomically distinct entities of the brain, body and physical environment. The computational principles underlying movement generation, such as optimal feedback control, have often provided robust behavioral insights, but failed to connect these to their corresponding neural implementations. The basic relationship between cortical activity and the corresponding movement thus remains poorly understood. This understanding can allow us to uncover the underlying neural computations and dynamics required to perform diverse movements in a dynamic environment.
The representational perspective posits that cortical activity relates to abstract movement representations, such as movement direction or muscle features [1–11]. In recent years, this representational perspective has been challenged by the dynamical systems perspective, which states that motor cortical activity mainly contributes to supporting intrinsic dynamical features, and exhibits representational tuning only incidentally [12–17]. Moreover, goal-driven models trained to perform a motor task considered analogous to that performed by the biological population of neurons have been successful at explaining population level dynamics [18–24]. However, these models lack biological and physical realism: they do not act on the real-time sensory feedback, do not interact with complex musculoskeletal dynamics, do not coordinate with intricate physical laws of the environment, are often trained on a very small subset of the features that the biological MC may consider as an input or output, and are not based on optimality principles. Therefore, these models often fail to generalize to unseen movements. For example, once a movement is learnt, we have the ability to produce this movement arbitrarily fast or slow, within a range, while maintaining accuracy. As a result, these models cannot be used for prediction of cortical activity during these unseen conditions. Building on previous studies [21–23, 25], we show that building in the biological realism of the musculoskeletal dynamics and incorporating the notion of optimality [26] in goal-driven models of MC can help better capture neural dynamics underlying movement generation while vastly improving generalization.
Here, we instantiate the point of view that the motor pathways in the brain act as a controller for the musculoskeletal system and drive optimal behavior under appropriate constraints [27–36]. Optimal feedback control theory has been successful at generating computational insights about the underlying control laws and providing behavioral-level predictions [37–52]. However, these models lack flexible neural network implementations of the controller, and fail to generate neural level predictions. A recent approach called MotorNet [53] includes training controllers with differentiable biomechanical models; however, the correspondence of the resulting controllers to the MC or generalization to unseen conditions has yet not been established. Moreover, neural constraints that arise from the evolutionary standpoint cannot be implemented using these models, such as minimization of neural firing rates.
Recently, deep reinforcement learning (DRL) has emerged as a promising way to train controllers for musculoskeletal models in highly complex tasks [54–58]. Using DRL, it is possible to train these models directly on high-dimensional sensory inputs to produce muscle activations that drive musculoskeletal models to produce complex movements, such as quick turning and walk-to-stand transitions [59]. The resulting behavioral-level dynamics have been shown to resemble experimentally recorded data [59].
In this research, we validate the neurophysiological plausibility of models obtained using DRL and propose a computational framework to generate behavioral- as well as population- and single-unit-level neural predictions for unobserved movements. We build in biological realism in this computational framework by basing it on the biological sensorimotor loop, incorporating anatomically accurate musculoskeletal models. We train these models using DRL as it shares the same notion of optimality (see Methods) as the optimal feedback control framework [60]. We also incorporate physical realism by simulating the resulting movements in highly advanced physics simulation engines. Moreover, for the first time, we show that the realistic neural-network-based implementation of a controller enables the implementation of neural constraints, such as minimization of neural firing rates, for movement generation.
This realism, along with the notion of optimality under suitable neural constraints, thus enables this framework to better capture and infer the behavior as well the underlying neural dynamics and single-unit firing rates. This framework also captures the computational principles and properties of motor control, such as generalization to novel movement conditions. Further, our results indicate that we can decode single-unit firing rates of MC using the activity of trained models even during unseen conditions. Thus, allowing these models to be used for hypothesis generation for prediction and analysis of neural activity during novel limb movements. We show that the developed framework is flexible enough to incorporate other known properties of MC, such as sensorimotor delays and forward models. This framework can thus significantly reduce the gap between the computational principles of movement generation and their corresponding population- and single-unit level neuronal implementation.
2 Results
2.1 Development of biologically and physically accurate goal-driven sensorimotor framework
We develop a goal-driven computational framework, µSim, to obtain a model of the MC with biological structural and functional properties (Fig. 1a). µSim consists of various independent, interactive and interchangeable modules (See Github page https://github.com/saxenalab-neuro/muSim). The modularity and flexibility allows the use of various musculoskeletal models, physics simulation engines and training algorithms with minimal modifications (Fig. 1b and Supplementary Fig. 1). µSim can thus be easily adapted for diverse use-cases, such as simulating realistic environmental conditions, building neural and behavioral constraints, and using musculoskeletal models of different animal species. Here, we explore the framework for a validated macaque limb model with 38 muscles [61], that was adapted from OpenSim to MuJoCo as in [62, 63] for computational efficiency.
µSim mimics the biological sensorimotor loop where the controller, representing the model of the MC, receives real-time sensory feedback and delivers muscle excitations to drive the musculoskeletal model producing diverse movements in the physics simulation environment (Fig. 1a). µSim consists of independent and interactive controller, environment, musculoskeletal model and goals specification modules (Fig. 1b). µSim supports the individual development and modification of these modules to incorporate biological and physical realism into the simulations. We base the controller on recurrent neural networks (RNNs) to mimic the recurrent connections of the biological MC. The controller consists of three layers, an RNN layer representing the model of MC, along with two feedforward layers representing the modular structure and the processing downstream and upstream of MC. The sensory feedback consists of visual and proprioceptive feedback, specifically muscle excitations, joint positions, joint velocities, hand and target positions and velocities and any high-level task input to the MC. We use anatomically accurate musculoskeletal models that receive muscle excitations to simulate the resulting movement. The physics simulation engine enables the implementation of diverse experimental conditions, such as contact / friction forces between bodies and ground reaction forces.
We use DRL to train the µSim controller to produce experimentally recorded movements under specified biological constraints (see Methods) (Fig. 1b). DRL trains the µSim controller to maximize a given reward function under specified neural- and kinematics-level constraints and is based on the same notion of optimality as optimal feedback control. In our approach, the reward function is designed such that the muscle excitations resulting in specified goals receive the maximum reward. The networks associated with reward function learning may be considered equivalent to other areas of the brain’s reward system, such as basal ganglia, and may reflect dopaminergic projections to MC [64]. Here, we used the soft actor-critic (SAC) algorithm in a maximum entropy framework to train the controller (see Methods).
2.2 Trained µSim controller achieves high kinematic accuracy
Using our computational framework, µSim, we detail out the procedure for one specific set of experiments where an adult male rhesus macaque monkey was trained to produce a cyclic movement of the forelimb at 8 different speeds, while kinematics and single-unit activity from premotor and motor cortex were recorded [22]. We used experimentally recorded kinematics during a subset of four different speeds as a task specification for µSim to train the controller, leaving out other speed conditions for testing purposes. It is important to note that the recorded neural data was not used at all during µSim training.
While several µSim controllers were trained with high kinematic accuracy, here we analyzed one example model for comparison with the recorded kinematics and neural data. The model achieved a high kinematic accuracy on all the four training speeds and we observed a significantly low mean squared error (MSE) between the experimental and simulated trajectories (Fig. 2a).
2.3 MC population dynamics and single-unit activity can be decoded from the trained µSim controller’s activity
We analyzed how well the trained µSim controller was able to capture the neural population dynamics of experimentally recorded MC neurons in the experimental conditions used for its training. We used canonical correlation analysis (CCA) for assessing correlations on the level of neural population dynamics and linear regression analysis (LRA) for assessing single-unit level decoding accuracy (see Methods).
CCA showed a high R2 between the µSim RNN activity and the experimental MC population activity for the conditions that the µSim was trained on (Fig. 2c). Next, we trained a linear regressor from the µSim RNN’s hidden activity to the recorded MC firing rates using data from all speed conditions except the held-out speed, and tested the linear model on the held-out condition using LRA (see Methods). The ‘seen’ speed is used for µSim training, but is held-out for LRA for finding the single-unit encoding fits below.
We compared the single neuron decoding accuracy (R2) of the µSim RNN to the existing goal-driven and representational MC models on the seen speed using LRA (Fig. 2e). We considered the following commonly used goal-driven- and representational-based models for comparison [1–11, 21, 22]: 1) Goal-driven open-loop RNN trained to transform a scalar speed signal to corresponding experimental recordings of muscle signals (electromyography; EMG) to provide a comparison against existing goal-driven models of MC ; 2) Empirically recorded EMG to test if the trained controller provides a comparison against the representational view that the cortical activity represents muscle commands; 3) Experimentally recorded kinematics to test if the MC model outperforms the representational view that signals mirroring the variables relevant to the task can capture cortical activity. The µSim RNN outperforms all alternative models in decoding accuracy for the speed condition tested (Fig. 2e). The firing rates of 4 example MC neurons reconstructed with LRA using different models for the seen condition are shown in Fig. 2h.
2.4 µSim generalizes to unseen conditions
We test the generalization ability of the µSim controller on ‘interpolated’ and ‘extrapolated’ speed conditions unseen during µSim training. The interpolated speed condition lies in between the range of training speeds and the extrapolated speed condition lies well outside the training speeds range. µSim, without further training, generalizes well to both the unseen conditions (Fig. 2b). We observed significantly low MSE values between the experimental and simulated trajectories averaged across x and y coordinates for the interpolated and extrapolated speed conditions, respectively.
After validating that µSim achieves a high kinematic accuracy on unseen speed conditions, we tested how well it captures the experimental neural data during those conditions. We first perform CCA separately on the interpolated and extrapolated speeds to recover a lower-dimensional subspace of maximum correlations between recorded MC and µSim RNN activity. We then perform inverse CCA to transform the µSim RNN activity from the maximum correlations subspace into the experimentally recorded MC activity subspace and reconstruct its top 3 PCs (Fig. 2d). CCA shows that the population activity of µSim RNN perfectly captures the MC population activity.
We then compared the decoding accuracy R2 of µSim against the traditional goal-driven and representational based models of MC for the unseen interpolated and extrapolated speed conditions. For single-unit decoding accuracy comparisons, we performed LRA separately for interpolated and extrapolated speeds with each speed treated as held-out speed for LRA. µSim outperforms the traditional models in explaining cortical neural activity (Fig. 2f and 2g). The firing rates of 4 example MC neurons reconstructed with LRA using different models for each held-out unseen condition are shown in Fig. 2i and 2j.
The difference in the relative decoding accuracy is even more pronounced for the unseen speeds as compared to the training speed conditions. These results support the dynamical systems view of the MC. Furthermore, µSim outperforms other simple goal-driven models confirming that the biological and physical realism plays an important role in explaining the cortical neural activity during limb movements.
MC population activity has been shown to exhibit strong rotational dynamics for various movements [20, 21]. We used jPCA (see Methods) to quantify the rotational dynamics in µSim population activity. jPCA revealed the presence of strong rotational dynamics in µSim population activity for all speed conditions with rotational dynamics explaining a significant variance of the µSim activity (Fig. 2k and 2l). The mean of distribution of angles between the µSim state and its derivative lies close to , indicating the presence of a strong rotational component. Moreover, for the monkey cycling task, MC population activity exhibits stacked elliptical neural trajectories having low trajectory tangling and separated along a speed dimension to enable flexible control of speed [22]. We found that µSim employs strikingly similar neural strategies for this task with elliptical trajectories separated along a speed dimension (Fig. 2m). Therefore, we conclude that various features of neural population dynamics, like rotational dynamics and low trajectory tangling, can be considered a consequence of optimally transforming the sensory feedback into muscle excitations under neural and musculoskeletal constraints to achieve desired movements. The emergence of biological neural population dynamics and structure thus allows µSim to capture cortical properties, such as generalization to unseen speed conditions, and the underlying computational principles.
2.5 Simultaneous goal- and data-driven modeling enables generalizable single-unit level decoding on unobserved conditions
We showed that under appropriate constraints, µSim is able to predict single-unit and population- level neural activity using post-hoc analyses such as linear regression and CCA. Here, we constrain the solution space of µSim further to develop a direct correspondence with the recorded single units. The resulting constrained µSim is denoted as ηµSim. We develop a dynamical systems model of MC that is both goal- and data-driven: it transforms sensory feedback into muscle excitations that are converted to kinematics, while directly encoding recorded MC single units. During training, we constrain a subset of the ηµSim RNN’s nodes to the experimentally observed MC single unit activity for training conditions (see Methods). We show that the ηµSim controller trained using DRL directly reproduces both the experimental MC single unit activity and the resulting kinematics (Fig. 3a). We then show that the trained model is generalizable to testing conditions that are significantly different from the training conditions (out-of-distribution generalization): it directly predicts the activity of the MC single units while simultaneously reproducing the behavior during those unseen conditions. The kinematics reproduced by ηµSim during these unseen conditions are shown in Fig. 3b.
To further interrogate the resulting structure of population dynamics in the trained networks, we use Procrustes Analysis (see Methods) to quantify how well the ηµSim RNN activity matches the experimentally recorded neural activity. Procrustes Analysis scales and rotates the network’s activity to align it optimally with the experimentally recorded single unit activity. Procrustes Analysis is a well-established technique for single-unit experimental and network activities comparison and is a more stringent test of alignment than CCA and LRA [23, 65, 66].
We compare the single unit encoding accuracy of the trained ηµSim controller against the existing data- and goal-driven models of the MC: 1). Data-Driven Encoding EMG-based MC model trained to simply reconstruct experimental MC single unit activity from the corresponding muscle excitations (EMG). 2). Goal-driven neural-unconstrained µSim trained without enforcing the direct neural activity reproduction constraint, as in the preceding Section.
The scatter plots (Fig. 3c and 3d) compare the R2 per neuron between its reconstruction produced by an alternative model (horizontal axis) against that produced by the ηµSim (vertical axis). We see that for all the seen and unseen speed conditions tested, most of the neurons lie above the dashed unity line, which means that they are being better predicted by the goal- and data-driven ηµSim model as compared to the alternative goal-driven neural-unconstrained µSim and data-driven EMG-based models. We also observe the lowest disparity D (L2 norm between the network reconstruction and recorded single-unit activity) for the ηµSim as compared to the alternative µSim and data-driven EMG-based models. Moreover, jPCA analysis reveals that the network dynamics of ηµSim mimic the MC neural population dynamics uncovering the underlying neural strategies for task representation (distinct initial conditions and trajectories) and movement generation (oscillatory dynamics) (Fig. 3e and 3f).
2.6 Modeling the contribution of sensory feedback and internal models on motor cortex activity dynamics and task execution
Here, we leverage the generalization ability of the ηµSim to explain and predict the computational principles underlying movement generation and their corresponding neural dynamics implementation. Little is known about how sensory feedback and internal forward models shape motor cortex (MC) dynamics during movement generation. We take a modeling perspective using the developed computational framework to probe this question. Goal-driven models of the MC based on recurrent neural networks (RNNs) when trained to transform high-level task-specific inputs into experimentally observed EMG exhibit rotational dynamics resembling the MC rotational activity patterns [20, 21]. Therefore, it is assumed that the intrinsic recurrent connections of MC give rise to its experimentally observed rotational dynamics with negligible contribution from sensory feedback. However, recent studies show that the models of MC based on feedforward networks without any recurrent connections also exhibit rotational activity patterns when trained to transform proprioceptive feedback into muscle activations, thus suggesting that sensory feedback contributes substantially to MC rotational dynamics [67]. We leverage the predictive ability of the ηµSim to evaluate these conflicting theories about the role of sensory feedback in MC dynamics.
To dissect the role of sensory inputs, we perform ablation studies by eliminating specific inputs to ηµSim and observe the effect on resulting movement kinematics and ηµSim’s RNN activity patterns (Fig. 4). We found that the recurrent connections’ ablation significantly disrupted the initial rotational dynamics (Fig. 4d) and the variance explained by jPCs significantly reduced from 47% to 31%. To compensate for this disruption ηµSim RNN employed angles greater than. Eventually, ηµSim is able to overcome this disruption and restore the rotational dynamics as the movement progresses (Supplementary Fig. 2a) and achieve the kinematics relatively well (Fig. 4b and 4c).
The ablation of the task-specific input significantly disrupts the rotational dynamics and separation of ηµSim trajectories across different conditions (Fig. 4e). However, the network is able to overcome the disruption of rotational dynamics as the movement progresses (Supplementary Fig. 2b) and achieve the task kinematics with negligible decrease in kinematics accuracy (Fig. 4b and 4c). The eventual emergence of condition-specific separation among trajectories (Supplementary Fig. 2b) indicates that ηµSim may be able to integrate the initially available sensory information to form task-specific internal forward models. The task-specific scalar input may thus reflect communication from the upstream brain regions representing the internal forward models that the network uses to separate the rotational trajectories across different conditions and execute the task relatively efficiently. Our results are consistent with previous studies suggesting the role of internal models in controlling the separation between neural trajectories [68].
Lastly, we examine the ablation of specific feedback to the model. If the proprioceptive feedback is ablated, ηµSim is not able to execute the task well (Fig. 4b and 4c). Surprisingly, the proprioceptive feedback ablation significantly increases the initial rotational dynamics (Fig. 4f). However, these rotational dynamics are significantly disrupted as the movement progresses (Supplementary Fig. 2c). If the proprioceptive and visual feedback is ablated, ηµSim is not able to execute the task at all (Fig. 4b and 4c). The ηµSim trajectories quickly settle at a fixed point (Fig. 4g) and the dynamics resemble pure translation/scaling as the mean of the distribution of angles approaches π. This correspondingly results in a linear movement of the musculoskeletal model’s hand (Fig. 4b). Similarly, the visual feedback ablation significantly disrupts the task kinematics and rotational dynamics (Fig. 4b and 4h). Our results are consistent with previous computational approaches to study the role of proprioception during motor control [69].
As the original rotational dynamics are significantly disrupted by sensory feedback, including both the visual and proprioceptive feedback, high-level task scalar and the recurrent connections ablation, we conclude that these rotational dynamics are modulated synergistically by sensory brain regions, internal models and task-specific modules to achieve the desired movement. Therefore, rotational dynamics alone do not guarantee efficient task execution.
Biologically relevant sensorimotor delays dictate that MC acts on delayed feedback information. If the musculoskeletal simulator is in the current state st, the µSim RNN receives an outdated version of the state feedback, st−Δt (Supplementary Fig. 3a). We observed that the µSim was able to achieve the delayed feedback task with high kinematic accuracy for Δt = 60ms [70]. We also observed a high R2 between the µSim RNN activations at time t and the non-delayed state st. This may stem from the µSim RNN developing a representation of the state st from the delayed state st−Δt in order to solve the task efficiently (Supplementary Fig. 3b).
2.7 µSim flexibly models diverse tasks and different animal species
Finally, we tested the modularity and flexibility of µSim for seamless modeling of diverse tasks and animal species. First, we tested the ability of µSim to capture neural dynamics across different tasks. For this purpose, we trained µSim on a center-out reaching task with 8 straight reach trajectories to the outer targets (see Methods) [71]. The controller achieves a high kinematic accuracy on all the experimentally observed conditions for the reaching task (Fig. 5a).
CCA revealed a high correlation between the trained µSim activity and the experimental MC population activity across all the reaching conditions (Fig. 5b and Supplementary Fig. 4a). LRA resulted in a high R2 between µSim activity and single-unit activity for the 8 conditions, indicating that the controller was able to capture the firing rates of experimentally recorded MC single units. We then compared the LRA-based single neuron decoding accuracy (R2) of µSim with the kinematics-based model of the MC (Fig. 5c and Supplementary Fig. 4b). We observed that µSim significantly outperformed the kinematics-based model of MC. Four example neurons reconstructed with LRA using different models are shown in Fig. 5d and Supplementary Fig. 4c.
Does the kinematic accuracy of the network determine the neural accuracy? We used Procrustes Analysis to quantify the role of the µSim’s kinematic accuracy in enabling neural decoding. We trained several µSim controllers with low to high MSE values between the simulated and experimental kinematic trajectories averaged across the reaching conditions. We observed that the controllers with relatively lower MSE (higher kinematic accuracy averaged across the 8 reaching conditions) had lower disparity between their activity and the experimental MC firing rates (Fig. 5e). This shows that the kinematic accuracy of the trained controller plays an important role in decoding the experimental MC neurons.
Next, we tested if a trained µSim is able to generalize across tasks (out-of-distribution generalization). For this purpose, we compared the kinematic accuracy of three different controllers on the reaching task: 1) The µSim trained and tested on the monkey reaching task, Reaching µSim 2) The µSim trained to perform the cycling task and tested on the reaching task, Cycling µSim. 3) A µSim controller initialized with random weights, Random µSim. The Cycling µSim significantly outperformed the Random µSim on the reaching task (Fig. 5f). We observed a relatively small difference between the kinematic accuracy of the Cycling µSim and Reaching µSim on the reaching task. This shows that the controllers obtained using this computational framework are able to generalize across diverse tasks.
Lastly, we trained µSim on a mouse alternation task with experimental single-unit MC activity, kinematics and EMG as in [72] (see Methods). We used an anatomically accurate musculoskeletal model of the mouse (see Methods) [73]. The µSim trained on the experimental kinematics achieved a high kinematic accuracy for the three training conditions (Fig. 5g). CCA revealed significantly high correlations between the top 5 PCs of experimentally recorded MC activity and the µSim RNN activity for the three conditions (Fig. 5h). We then used LRA to assess correlations on the level of single-unit activity. LRA resulted in high correlations between the µSim activity and experimental single-unit activity. Four example neurons reconstructed with LRA using different models are shown in Fig. 5i for the three conditions. LRA-based single neuron decoding accuracy (R2) comparison of µSim with the EMG- and kinematics-based model of the MC showed that µSim significantly outperformed the alternative models (Fig. 5j and 5k). Therefore, µSim can be used in a seamless way across well-researched animal species in neuroscience.
2.8 Discussion
The brain is evolved to optimally interact with the physical environment under neural and musculoskeletal constraints. Here, we develop an anatomically and physically accurate computational framework, µSim, based on the biological sensorimotor loop, to obtain a model of MC. Under specified neural constraints, we use DRL to train the models of MC to optimally achieve diverse movements. The trained controller is able to capture and predict the neural dynamics and single-unit activity underlying movement generation in a generalizable manner. Constraining a subset of the µSim RNN’s nodes to the experimental data during training conditions further enables it to develop one-to-one correspondence with the experimentally recorded single-units of MC and enables their direct prediction during unseen conditions. The resulting models capture the experimentally recorded neural population dynamics and structure surprisingly well. These models may capture higher-cognitive properties such as generalization to unseen conditions, even to those that are significantly different from training. We finally show that the developed framework accurately captures the computational principles governing the integration of sensory feedback, task specification, and recurrent neural dynamics underlying movement generation. Therefore, this computational framework can be used for prediction and hypothesis generation on both the behavioral and neural level. Future directions include the investigation of network remapping during the learning of new behaviors [74, 75].
µSim consists of various independent and interactive modules, such as the musculoskeletal model, physics simulation engine, controller and optimal training algorithm, experimental data and goals specification. These modules can be used and enhanced independently by machine learning engineers, computer scientists, control theorists, roboticists, biomedical engineers, medical professionals and neuroscientists to enable diverse applications. µSim can enhance our understanding of motor control on both behavioral and neural level, and uncover the role of sensory, cognitive and reward regions of the brain in movement generation. µSim will enable significant applications, such as enabling enhanced real-time decoding and encoding models for brain computer interfaces, discovering muscle excitations or designing stimulation signals to achieve various movements, building predictive simulations to improve rehabilitation and assistive devices research, discovering the effect of neural stimulation on resulting movement in deep brain stimulation, and building biologically plausibile control and machine learning models.
3 Methods
3.1 Experimental Data
3.1.1 Monkey Cycling Data
Head-restrained rhesus macaque monkeys were trained to perform a cycling task while sitting in a customized chair. During the controlled experiments, monkeys manipulated a pedal-like device with their right-arm while the left arm was loosely restrained. The real time horizontal and vertical hand positions were recorded. The wrist movement was restrained and the cycles were driven mainly by the movement of the elbow and shoulder. This resulted in a highly stereotyped arm movement across cycles.
The monkeys rotated the pedal to track a moving target with respect to their virtual first-person location on the monitor in front of them. Juice reward was dispensed so long as they maintained their virtual position close to the displayed target.
Conventional single electrodes driven by a hydraulic microdrive were used to make sequential neural recordings from a broader range of sites in primary motor cortex and adjacent aspect of dorsal premotor cortex. For all further analyses, these recordings were treated together as a single motor cortex population. Neural signals were amplified, filtered and manually sorted using Blackrock Microsystems Digital Hub and 128-channel Neural Signal Processor. During cycling, nearly all of the isolations made were responsive and those with low signal-to-noise ratios or insufficient trial counts were discarded. Gaussian filtering (20ms SD) was then applied to filter the spikes of each recorded neuron to produce estimated firing rate, which was then trial averaged.
The trials were divided into 8 speed bins. Kinematics, EMG and neural data were then averaged across the trials in each speed bin. In each speed bin, a high degree of training ensured stereotypical behavior across trials. Adaptive alignment procedure was applied to correct remaining slight misalignments. Aberrant trials deviating from stereotypical behavior were discarded. Kinematics, neural and EMG were then averaged across the resulting trials in each speed after applying this adaptive alignment procedure. Further details about data preprocessing are given in [22].
3.1.2 Monkey Reaching Data
A rhesus monkey was trained to perform a center-out eight-target reaching task while grasping a two-link manipulandum. The outer-targets were spaced 45° intervals around a circle. The monkey had to reach the outer target and hold for a liquid reward. The neural data was recorded from proximal arm area of primary motor cortex and dorsal premotor cortex contralateral to the arm used for performing the task. Kinematics and neural data were then averaged across stereotypical trials from the same session.
Further details about behavior and neural data recording and preprocessing are given in [71].
3.1.3 Mouse Alternation Data
Alternation task elicited MC dependent flexor-extensor alternation. Mice were head-fixed and trained to perform alternation task for a water reward while muscle activity from the forelimb was recorded. Neural activity in the left caudal forelimb area was recorded after mice had been fully trained. Kinematics, EMG and neural data were then averaged across stereotypical trials from the same session.
Further details about behavior and neural data recording and preprocessing are given in [72].
3.2 Musculoskeletal Models
3.2.1 Musculoskeletal Model of a Monkey Limb
We build on the 3D musculoskeletal model of a macaque monkey arm developed in [61]. The developed model consists of seven degrees of freedom (DoF) that include extension and flexion of the elbow, 3D rotation about the shoulder joint, supination and pronation of the lower forelimb, adduction/abduction and flexion/extension of the wrist. The five segments represented in the model are the hand, the torso, the radial side of the lower arm, the ulnar side of the lower arm and the upper arm. Shoulder abduction and adduction is its rotation around the x–axis and its rotation around y-axis is defined as internal and external rotation. Rotation of the shoulder about the z-axis is shoulder flexion and extension. Elbow flexion is rotation about the z-axis. An intermediate x rotation is used to define the off-axis pronation and supination axis. Rotation about the x-axis at the center of the wrist is wrist flexion and extension. It is followed by rotation around the z-axis which is wrist abduction and adduction.
The model also consists of 38 muscles. These muscles are based on the anatomical data obtained from cadaveric studies [76, 77]. The muscle properties such as muscle/tendon length and pennation angle were obtained from the literature [78].
The musculoskeletal dynamics used for developing the forward dynamic model can be represented by the following equation Where M represents the mass distribution of the system given the current joint angles Θ. and represent the joint velocities and accelerations, respectively. R is the moment arm matrix. F is a vector of muscle forces. G, V and E describe the moment contributions of gravitational, internal and external forces, respectively. The model is designed to be scalable to a generic monkey arm given the monkey mass. Further details of the model are given in [79].
This model is adapted for use in OpenSim [80]. We first replaced the existing muscle model with a more biologically accurate and computationally stable Millard muscle model [81]. However, the forward simulations based on biologically accurate musculoskeletal models in OpenSim are computationally expensive. This makes the downstream learning of the controller extremely inefficient. To overcome this challenge of slow forward simulations, we first convert the monkey arm model from OpenSim to an equivalent model in MuJoCo physics simulation engine [62, 63, 82]. MuJoCo is a state-of-the-art physics-based joint-constrained physics simulation engine that can reach forward simulation speeds of more than 600 times those of OpenSim [62]. Therefore, we optimized the OpenSim monkey arm model for use in MuJoCo by approximating the musculo-tendon units: we minimized the squared difference of joint positions between the OpenSim and MuJoCo forward simulations over all training trajectories. This resulted in an anatomically accurate musculoskeletal model in MuJoCo that can facilitate the downstream learning of the controller by producing fast forward simulations.
3.2.2 Musculoskeletal Model of a Mouse
We utilized an anatomically accurate mouse musculoskeletal model developed in [73], and built in the PyBullet framework [83]. PyBullet is a fast, stable, and open source physics simulator with a python application programming interface (API). The authors of [73] developed a simulator-agnostic muscle library to be integrated with PyBullet for the purposes of this model. µSim controller interacts with the right forelimb of this model, directly activating its 18 muscles while the rest of the skeletal model is fixed in place. Forelimb muscle attachment points were determined based on embryonic studies of mice; due to lack of experimental data, mainly distal muscles were included. The model does not include proximal muscles that originate from the spinal segment.
We now list the DoF for each joint in the forelimb. For shoulder joints, which were modeled as a spherical joint, there are three rotational DoF: retraction-protraction (rotation about the transversal axis), abduction-adduction (coronal axis), and external-internal rotation (sagittal axis). For elbow joints there are two DoF: extension-flexion (transversal axis) and supination-pronation (sagittal axis). Wrist joints additionally contain two DoF: extension-flexion (transversal axis), and abduction-adduction (coronal axis). More details as well as the joint range-of-motion and limits can be found in [73].
The muscles themselves are of the Hill-type [84], with activation dynamics given by where u(t) are the activation dynamics, a(t) is the muscle excitation, and τact is the time constant. The muscle excitation signal a(t) is determined by µSim controller, directly controlling the muscle activation in order to produce the necessary movements. We chose a biologically realistic starting position of the forelimb based on the joint range-of-motion given in [73].
3.3 Kinematics MSE
The MSE between the kinematics produced by µSim and experimental kinematics is calculated as follows: where x and y represent the hand’s x and y coordinates, respectively. h∗ and h represent experimental and µSim hand kinematics, respectively. T represents the total number of time points in the corresponding cycling speed trajectory.
3.4 Formulation of the Reinforcement Learning Problem
We formulate the task of controlling musculoskeletal models to perform different movements as a reinforcement learning problem. This formulation consists of a policy network, critic networks and an environment. Given the current state st ∈ S of the environment, the policy network outputs the probability distribution, represented by πθ(a|s), over the possible actions a ∈ A. θ denotes the parameters of the policy network implemented using RNNs. The environment can be implemented using different physics simulation engines. In this work, we use MuJoCo [82] and PyBullet [83] physics engines for implementing the environment. The environment consists of musculoskeletal models that can interact with other physical objects. Advanced physics engines allow the implementation of physically realistic simulations under different constraints such as contacts. At each timestep t, the environment receives an action at from the policy network. at represents the muscle excitation signal or motor command. Given at, the resulting movement is executed and the environment transitions from the current state st to next state st+1. st+1 represents the physical effects of the motor command executed in our biomechanical simulation. The state transitions can be stochastic in µSim and the conditional probability of reaching the next state st+1 is given by p(st+1|st, at). The state transitions are governed by the underlying dynamics of the musculoskeletal system and the physics environment. These dynamics are defined by differential equations represented by F : Given the probability of starting in some initial state s0, the probability of realizing a trajectory T = (s0, a0, ……, sN) under the policy πθ(a|s) is given by: In current work, we implement deterministic transitions of the environment. At each timestep t, a reward rt is also generated to quantify the performance of the policy network.
3.5 Environment Design
The environment consists of anatomically accurate musculoskeletal models of different subjects that can interact with other objects while obeying physical laws and constraints, such as collisions, friction, torque, contacts and ground reaction forces. This allows the simulation of freely moving subjects that can interact with their environment. This framework can thus be used to investigate motor control for freely moving subjects. Motor control for freely moving subjects has not been explored well previously due to limitations in training algorithms and simulation frameworks.
In current work, the objective is to reproduce experimentally observed behavior. The musculoskeletal models of different species in this work are fixed in space while the limbs/extremities and thus the end-effector can move freely. The experimentally recorded kinematics are represented in the simulation using abstract objects known as targets. At each timestep t, the position of the target is updated to reflect the experimentally recorded position to be tracked by the end-effector in a 3D plane. The performance of the policy network at each timestep t can thus be quantified by the distance between the end-effector position and the target position. The policy network is then trained to output muscle commands such that the end-effector tracks target position at each timestep t, thus reproducing the experimentally recorded behavior.
The same framework can be used to model the task or abstract behavior in the absence of experimentally recorded kinematics using sparse metrics for quantifying the performance of the policy network. For example, in the absence of experimental kinematics, the reaching task can be achieved by fixing the abstract target at the desired final position.
3.6 States/Actions
The environment state se(t) is filtered to form the state feedback st for the policy network. st may consist of complete or partial information required to solve the control problem. If partial information about the environment state se(t) is available such that the state feedback is imperfect or partially observable, the control problem in the DRL framework can be formulated as a partially observable Markov decision process (POMDP). Various formulations of DRL, such as deep hierarchical RL, can be used to solve the control problem for POMDPs [85]. In previous work, a DRL formulation with RNN implementation of the policy network has also been shown to solve control problems with a partially observed state [86]. This also justifies the use of RNN implementation of the policy network in current work, in addition to the existence of strong recurrent connections in the motor pathways. This may also point to the fact that the motor pathways may be able to construct complete state from partially observed information through these recurrent connections. In a dynamic environment with freely moving subjects, often only partial state information is available. Therefore, this computational framework can be used to solve such complex tasks.
In this work, st represents the sensory feedback and consists of both the visual and proprioceptive feedback. Specifically, state feedback st to the policy network consists of muscle kinematics (muscle activations at the last timestep), joint angle and velocity for each joint, positions and velocities of the end-effector in 3D space, positions and velocities of the target in the 3D space, the difference vector representing the distance between the end-effector and the target, and a scalar representing different conditions within a task (‘task information’).
The action vector a ∈ A represents the muscle excitations or the motor command that are then transformed into muscle activations. These activations actuate some DoF by applying motor torque. Physics simulation engines, such as MuJoCo, can also accurately model additional active and passive forces, including gravitational and contact forces, in addition to these actuated forces. Such additional forces are not modeled accurately using biomechanics engines, such as OpenSim. Therefore, such biomechanics engines can not be used to model neural dynamics accurately in the context of freely moving subjects interacting with each other or the dynamic environment. The developed computational framework is therefore expected to advance our understanding of motor control in such complex tasks.
3.7 Reward Function Design
The reward function r(st, at, nt) specifies the behavior of the policy network. Motor pathways drive optimal behavior under specific neural and behavioral constraints. Therefore, the reward function consists of two parts: task specification T (st, at) and constraints specification C(st, at, nt). where, nt represents neural parameters, such as synaptic weights or firing rates. In the developed computational framework, these neural parameters nt correspond to the parameters θ and activity of the policy network. To keep the notation concise, we assume nt ⊂ st.
The generated rewards can be immediate, such as provided at each timestep t, or they can sparse, such as provided only at the last timestep. Freely moving subjects in dynamic environments often receive sparse rewards that can be modeled using the proposed computational framework.
It has been proposed that the central nervous system (CNS) achieves a specific task under various kinematics- and behavioral-level constraints. Previously, constraints based on the minimization of muscle effort [26, 87] or maximizing the smoothness of the end-effector trajectory and of the torque commands have been proposed [88, 89]. The models based on these constraints have been successful at reproducing empirical kinematic data. However, it is not clear how these are measured by the CNS. Other models assume noise in the motor command that scales with its magnitude [90]. The CNS aims to then minimize the end-point variance given this motor noise. However, is also not clear if such scaling is a consequence of neural constraints, such as neuronal noise. Therefore, it remains an open question how we approach this constraint for novel, unrehearsed movements [87].
In this work, we instead assume and validate the existence of neural constraints, such as minimization of neural firing rates, that are suitable from an energy minimization and evolutionary standpoint. These neural constraints can also be thought of giving rise to other previously proposed kinematics- and behavioral-level constraints. We also validate that such constraints generalize across novel and unrehearsed movements. To test specific hypotheses, this framework can be used to implement various behavioral- and neural-level constraints, such as the ones described above. After training, the resulting behavior and network activity under specified constraints can then be analyzed or compared with experimental data to validate their existence.
For the purposes of this work, these neural constraints can be implemented using regularizations on the policy network and are described below. Here, we design the reward function as consisting purely of the task specification T (st, at). In this work, the task specification is designed to make the end-effector reproduce the empirical kinematics by tracking the target at each timestep t: where and h z,t represent the x, y and z coordinates of the end-effector’s position and and represent the x, y and z coordinates of the target’s position at timestep t. This reward function r(st, at) produces the maximum reward rt for the motor command that results in the minimum distance between the end-effector and the target for the given timestep t. rt is designed to decay exponentially with the increasing distance between the hand and the target to accelerate the learning of the policy network parameters. We use wd = 5.0.
3.8 Notion of Optimality
Optimality principles have been used in motor control to specify the control laws giving rise to the generated behavior and to theoretically explain why the motor system behaves as it does [91]. Given the dynamics of the motor system, musculoskeletal models and environment as in (4), these control laws are usually specified by the cost functions, c(st, at), that can be considered as analogous to the reward functions r(st, at). For example, the control laws that achieve the task as accurately as possible while minimizing the energy consumption are considered more suitable than those that do not follow either of these constraints. The constraints that result in controls laws under which the experimental behavior can be inferred provide insights into quantities that the motor system is trying to optimize. Such constraints, such as minimization of energy consumption, should thus help us generalize across different novel and unrehearsed movements. Usually such constraints are specified in kinematics space, such as maximizing the smoothness of end-effector trajectory. The realsitic neural network implementation of control laws provides a way to instead specify these constraints in the neural space, such as minimization of neural firing rates. We show that the experimental behavior and neural dynamics across diverse and novel movements can be inferred under such neural constraints together with the biological accuracy of the developed sensorimotor loop.
3.9 Soft Actor-Critic (SAC) Design
We adapt the SAC algorithm to train the policy and critic networks to achieve the desired task. SAC is an off-policy RL algorithm based on the maximum entropy framework [92]. This algorithm is chosen because of the complex environmental dynamics and computational complexity of the biomechanical simulations.
The reward rt generated at each timestep t determines the behavior of the policy trained using this algorithm under the neural constraints. The sum of the rewards discounted by the factor γ defines the return R. Where T is the last timestep. The discount factor γ < 1 determines the relative importance of the future rewards relative to the earlier rewards. We set γ = 0.99 for training µSim.
The high dimensional action space and the complex environmental dynamics in biomechanical simulations make exploration quite inefficient using standard RL algorithms. Here, an entropy term H(πθ(.|s)) is incorporated in the return to achieve efficient exploration required for training. The entropy term determines the stochasticity of the trained policy and is used to achieve efficient exploration and robust convergence towards the global optimum given the high dimensional action space and complex environmental dynamics. Moreover, it has been shown to considerably improve learning speed and stability over the training algorithms that maximize the standard RL objective of maximizing only the return R. Probabilistic matching that has been used to explain human decision making can be considered the biological equivalent of the maximum entropy RL framework [93].
The objective is to learn a policy that maximizes the following soft return: The temperature coefficient κ controls the stochasticity of the trained policy by determining the relative importance of reward against the entropy term. κ is adjusted automatically during training using the dual gradient-descent implemented in the SAC algorithm [92].
The SAC algorithm uses policy iteration to train the parameters of the policy network to maximize the objective function J(θ). To reduce the variance of the sampled trajectories during training, the SAC algorithm makes use of the critic network in addition to the policy network.
3.10 Critic Network
The critic network represents the transformation from the state-action pair (st, at) to its q-value q(st, at) which in turn represents the expected soft-return of motor command at given the sensory feedback st in the maximum entropy RL framework [92]. The critic network parameterizes the q-value qϕ through parameters ϕ. In this work, the critic network consists of two RNN-layers with 256 nodes each.
Dopaminergic projections from the ventral tegmental area to the MC may reflect the neural correlates of the reward signals [94]. The reward regions of the brain, such as basal ganglia along with these dopaminergic projections, can be considered equivalent to the critic network. Disruption of these reward signals can inhibit further learning during forelimb reaching movements in rats, a prediction that can be obtained using the actor-critic architecture used [95].
3.11 Policy Network
The actor/policy network parameterizes the policy πθ through the parameters θ. The policy network represents the mapping from state feedback st to the probability distribution over the muscle excitation space. In the brain itself, various regions are involved in the transformation from the sensory feedback to the motor command. The sensory feedback is first processed in the sensory and visual processing regions of the brain that share strong reciprocal connections with the premotor areas. The premotor cortex processes the feedback further and has strong reciprocal connections to the MC. The MC then transforms the processed feedback into motor commands through subcortical and spinal cord projections.
We base the architecture of the policy network to mimic the modular structure and reciprocal connections of the motor pathways involved in the movement generation. Therefore, the policy network consists of three layers representing the sensory and premotor/motor regions and the final subcortical and spinal cord projections. The layers representing the premotor and motor regions are based on RNNs to mimic the recurrent connections. For comparison with the recorded cortical data, we use the activity of the RNN layers representing the premotor and motor cortex.
The policy network consists of three layers. The first layer is a feedforward forward layer with the following input-output transformation: where sk(t) is the input sensory feedback with dimensionality I1. σ1 is the non-linearity for the first layer. WU and bU represent the weights and biases for the first layer.
The second layer consists of an RNN. RNN can be considered a discrete dynamical system with the following dynamics: with Where r represents the inputs to the non-linear activation function σ2 and x represents the corresponding output of RNN hidden layer. N1 is the number of units in the first feedforward layer and N2 is the dimensionality of the hidden layer of the RNN. σ2 is the non-linearity for the RNN layer. WI and bI represent the input weights and biases, respectively. WH and bH represent the recurrent weights and biases for the RNN layer, respectively.
The final layer is a feedforward layer with the following input-output transformation: The output of the third layer, z, represents the muscle command with dimensionality N3. WZ and bZ represent the weights and biases for the final feedforward layer, respectively.
3.12 Neural Constraints
Here, we hypothesize that the brain is evolved to produce optimal behavior under neural constraints. These neural constraints can be implemented as regularizations on the networks. Here, we implement three neural constraints as regularizations on the policy network activity.
The first regularization term is an L2 penalty on the input and output weights of the three layers which encourages the minimization of synaptic weights or sparse connections between the network nodes. The second regularization term encourages the minimization of the neural firing rates. It also prevents the network activity from saturating: Where T consists of cumulative timepoints for all the conditions the network is trained on.
Finally, we implement a third regularization term that encourages the network to achieve the task while making simple trajectories in the low-dimensional state-space as proposed in [21]. In this work, we observed that RSD does not play a significant role in increasing the correlation or similarity between the network and neural activities.
Therefore, the final loss function consists of the following terms: Where θ = WU, WH, WI, WZ, bU, bH, bI, bZ represents all the parameters of the policy network. The policy network parameters are then trained using gradient descent to minimize the loss L(θ). We use Adam optimizer to update the parameters of the loss function. This loss function is used for µSim training.
We used α = 0.001, β = 0.01 and ζ = 0.1. We used N1 = N2 = 256. Additionally, we used σ1 = σ2 = tanh.
3.13 Simultaneous goal- and data-driven modeling
To implement the simultaneous goal- and data-driven modeling, we enforce an additional constraint on a subset of the RNN units to follow recorded neural activity for the training conditions. Specifically, we minimize the follow loss between the network and recorded single unit activity: Where Tc represents the number of timesteps per condition, C is the total number of training conditions, n represents recorded neural activity, and NR is the total number of recorded neurons. At each timestep t, the simultaneous goal- and data-driven modeling loss minimizes the difference between the activities of the subset of RNN units and the corresponding recorded neurons. The final loss is: This loss function is used for ηµSim training.
We used τ = 104 to account for relatively less number of experimentally recorded neurons as compared to total number of units in the policy network and relatively smaller magnitude of LGDM as compared to other terms in (20).
3.14 CCA
We used CCA to find the correlations R2 between the trained network and experimentally recorded neural population responses [96]. CCA finds the weightings for the units in network and experimental datasets such that the reweighted datasets are maximally correlated. Both the network and recorded neural activities were first reduced to ten dimensions using principal components analysis (PCA). CCA was then applied to recover the subspace of maximum correlations between the recorded and network activities during one complete cycle of the movement. Inverse CCA was applied to transform the network activities from the maximum correlations subspace back into the recorded activities subspace. We reported the reconstruction comparison and R2 between the network and experimentally recorded activities in this subspace.
3.15 jPCA
We used jPCA [20] to quantify the oscillatory dynamics in neural state x(t, c) across times, t, and conditions, c. jPCA provides summary features, such as quality of fit and variance explained, relevant to the hypothesis that the neural state evolves according to the oscillatory dynamics. It also allows the visualization of the two-dimensional projection of the neural data containing the oscillatory dynamics. As a preprocessing step, we first applied PCA to the neural data having dimensionality equal to the number of recorded neurons to reduce it to the top 4 PCs capturing ≥ 90% variance. Using jPCA, we analyzed two different time periods of the neural activity: 1) 150ms after movement onset; 2) 620ms after movement onset.
3.16 Linear Regression Analysis
We used linear regression analysis (LRA) to compare the network and experimentally recorded neural single-unit level activities. For the monkey cycling task, we first fit a linear model with ridge regression on all the conditions except the held-out condition. The network and neural activities for training conditions excluding the held-out condition are first concatenated along the time dimension separately. Ridge regression is then used to fit a linear model in which the concatenated activity for each recorded neuron is determined by the concatenated network activity separately. For testing, we used this trained linear model to transform the network activity for each held out condition into the recorded neural activity subspace. The transformed network activity is then compared with the corresponding actual single-unit level neural activity for the held-out condition during the movement period to find the correlations. We used a similar procedure for comparing the kinematics and EMG models.
For the monkey reaching and the mouse alternation task, we fit a linear regressor using ridge regression separately on each condition. Ridge regression is used to fit a linear model in which the activity of each recorded neuron is determined by the network activity separately for the movement period. We used a similar procedure for comparing the kinematics and EMG models. We used a regularization coefficient of 5 × 10−2 for LRA.
3.17 Procrustes
We used Procrustes Analysis to compare the network and neural single-unit level activities [65, 66]. Procrustes applies linear transformations, such as scaling and rotation, to the network activity to align it with the neural activity for the given condition. Procrustes minimizes the following loss: Where datanet represents the network activity and dataneu represents the recorded neural activity. The rows of the data matrix correspond to the timepoints and the columns correspond to the units/neurons. If the number of units in the neural data is less than the network activity, we first apply PCA to the network activity to make the number of units equal. Conversely, if the number of units in the network activity is less than the neural data, we first append the columns with zeros to make them equal. Here, D is also known as the disparity.
4 Supplementary Figures
2.9 Acknowledgements
This work was supported by NIH Brain Initiative grant 1RF1DA056377-01. We are very grateful to the following researchers for making their experimental data available: Abigail Russo and Mark Churchland for the monkey cycling dataset, the Slutzky laboratory for the monkey reaching dataset, and Claire Warriner and Andrew Miri for the mouse alternation dataset.
Footnotes
Acknowledgements for datasets used are updated.
References
- [1].↵
- [2].
- [3].
- [4].
- [5].
- [6].
- [7].
- [8].
- [9].
- [10].
- [11].↵
- [12].↵
- [13].
- [14].
- [15].
- [16].
- [17].↵
- [18].↵
- [19].
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].
- [29].
- [30].
- [31].
- [32].
- [33].
- [34].
- [35].
- [36].↵
- [37].↵
- [38].
- [39].
- [40].
- [41].
- [42].
- [43].
- [44].
- [45].
- [46].
- [47].
- [48].
- [49].
- [50].
- [51].
- [52].↵
- [53].↵
- [54].↵
- [55].
- [56].
- [57].
- [58].↵
- [59].↵
- [60].↵
- [61].↵
- [62].↵
- [63].↵
- [64].↵
- [65].↵
- [66].↵
- [67].↵
- [68].↵
- [69].↵
- [70].↵
- [71].↵
- [72].↵
- [73].↵
- [74].↵
- [75].↵
- [76].↵
- [77].↵
- [78].↵
- [79].↵
- [80].↵
- [81].↵
- [82].↵
- [83].↵
- [84].↵
- [85].↵
- [86].↵
- [87].↵
- [88].↵
- [89].↵
- [90].↵
- [91].↵
- [92].↵
- [93].↵
- [94].↵
- [95].↵
- [96].↵