Neuromuscular Control With Forward Dynamic Approximation In Human Arm

The central nervous system (CNS) directs a large number of muscles to produce complex motor behavior’s. Moreover, human movement control is significantly compromised in neuromuscular diseases; many of which result from the imbalance in the sensorimotor control system. Forward models are also known to capture the casual relationship between the inputs to the system and its output. The aim of this study was to approximate a forward dynamic simulation by using feed forward neural networks in human elbow arm movement. Motion capture data was used to generate c3d data from Qualisys software via Qualisys track manager. The c3d data was then converted into marker data which contains the markers location in the form of x-y-z coordinate. OpenSimbiomechanical software was used to process the marker data for scaling, inverse kinematics and the subsequent forward dynamics simulation in the right human upper arm model. The results for the training error of the approximated forward dynamics simulation for 3 degrees of freedom (DOF) were 2.804429, 1.468017 and 2.475500 with a constant validation error of 0.000000. The proposed control algorithm served as robust method for approximating a forward dynamic simulation, and can broaden the understanding the representations of neuromuscular control in the central nervous system as representations of movement plans that are eventually executed by the spinal cord and muscles in the periphery. .


INTRODUCTION
Human movement control is significantly compromised in neuromuscular diseases; many of which result from the imbalance in the sensorimotor control system. Muscle excitation that follows the neural impulses to generate movements common in daily activities like reaching, grasping holding plays a key role in the mechanism of motion from the sensory input 1,2 . Neural Networks have also gained popularity in their usage in fields such as Natural Language Processing, Computer Vision, reinforcement learning, and even animation 3,4 . The primary motivation of this study was similar previous studies which discuss the control of multijoint movements by the CNS to perform complicated transformations from the desired behavioral goals into appropriate neural commands to muscles 5 . A model thus formulated was proposed by defining an objective function, which is considered as a performance measure for any possible movement: square of the rate of change of torque integrated over the entire movement [6][7][8] . Furthermore, many studies have also discussed the use of neural network model was developed to understand the internal organization of movement generation [9][10][11] . The use of recurrent neural networks is discussed in detail where recurrent neural network is trained to produce the muscle activity of the reaching monkeys resulting in a dynamics where simple inputs were transformed into temporally and spatially complex patterns of muscle activity 12 .Past researches in the field of motor control have also discussed the criterion that is adopted for determination of the trajectory wherein the common invariant features were found in several researches that measured the hand trajectories of skilled movements 13 . The hand movement for instance, between a pair of targets, subjects tended to generate roughly straight hand paths with bellshaped speed profiles. These observations along with dynamic optimization theory, led to the design of a mathematical model which accounts for formation of hand trajectories 13,5 .The use of neural networks which contains internal model of inverse dynamics as essential learning parts of the musculoskeletal system was the main feature in the learning of voluntary motor control 14 . Other studies emphasize the forward model as representation of motor system that predicts the next state using the current state of motor system and motor commands 15 .
This study will discuss its application to approximating and controlling muscle activation with feed forward neural networks. The only goal of this study is to see if we could achieve values using the above strategy in order to have a broader understanding as to how neural signals of the motor cortex relate to movement. We have selected a musculoskeletal model that we will use for our neural network approximation. In order to achieve our goal, we generate forearm flexion movements in an allowable range in order to train and test our neural network.

MATERIALS OpenSim
OpenSim is a software framework designed to build and share musculoskeletal models, simulate movements, and analyze and visualize those movements using specialized tools. This framework contains, amongst other things, a graphical user interface (GUI) scripted in java language, various built-in musculoskeletal models developed and published by an open source community of researchers and users, a software development kit, APIs and other various tools that are used by researchers or hobbyists for analysis 16 . OpenSim has many open source musculoskeletal models available in their repository. These include human and animal models that can vary in complexity and scope. In this study, we are only looking at the right arm of a human which consists of 2 degrees of freedom, with 6 muscles and muscle elements such as bones, ligaments, etc. and provides a dynamic graphical simulation to visualize and analyze moments about the joints of a human arm 16 . There are multiple arm models that are available and we have the liberty of choosing one that we believe can achieve our objective simply and effectively. There are many tools available on the platform that can be used to analyze muscle actuation and joint position. We use an arm model whose muscles are based on Hill's muscle model 16,17 , which is shown in fig 1.

Motion Capture System
The motion capture device used for this study was Qualisys motion capture system (Qualisys Inc) with Qualisys Track Manager (QTM) allows users to perform 2D, 3D and 6DOF capture of data in realtime by the way of C3D data which is then converted into in a .trc file which provides the distances between markers on the model and experimental marker positions which are compared and determined by the scaling factors 6 . For this study the scaling factor was set to 1.0 which is the default scaling factor used in OpenSim software. The experimental markers set for 15 healthy subjects were right acromium; right humerus epicondyle and right styloid process of radius to match with the generic marker in model. The unscaled model should have a set of virtual markers placed in the same anatomical locations as the experimental markers. The modeling of the neural network architecture was carried out in Matlab 2013 platform. Matlab is a high level language and interactive environment for numerical computation, visualization and programming. MATLAB provides a range of numerical computation methods for analyzing data, developing algorithms, and creating models .The inverse and forward simulation along with the conversion of motion capture data were performed via Matlab-OpenSim interface 18,19 .

Scaling
S c a l i n g i s p e r fo r m e d b a s e d o n a combination of measured distances between x-y-z marker locations which were obtained from the data specified in the track row column (trc) files and manually-specified scale factors. The marker locations are usually obtained using motion capture equipment. For pair 1 on the model (m1) the distance is computed by placing the model as all joint angles get assigned their default values. The experimental distance between pair 1 (e1) is computed by looking at each frame of experimental marker data in the given .trc file, computing the distance between the pair for that frame, and taking the average across all frames in a user-specified time range. The scale factor due to pair is s1=e1/m1.

METHODS
The schematic of the proposed approach to perform neuromuscular control with forward dynamic approximation is shown in Fig.2

Inverse kinematics
The Inverse Kinematics (IK) tool computes generalized coordinate values which positions the model in a pose. This pose "best matches" experimental marker and coordinate values for that time step or the frame of motion 1,15 .Mathematically, the "best match" is expressed as a weighted least squares problem, whose solution minimizes both marker and coordinate errors. The marker error is the distance between an experimental marker and the corresponding generic marker on the model; the IK solver computes the generalized coordinates and is used to position the markers. Each marker has a weight associated with it, which specifies the extent to which that marker's error term should be minimized. The IK solver is given by the following equation as the weighted least squares problem solved by IK.
... (1) where q is the vector of generalized coordinates being solved for, x i exp is the experimental position of marker i, x i (q) is the position of the corresponding marker on the model (which depends on the coordinate values), q j exp is the experimental value for coordinate j. The marker weights (w i 's) and coordinate weights (É j 's). This least squares problem is solved using a general quadratic programming solver, with a convergence criterion of 0.0001 and a limit of 1000 iterations.
This generates motion file with the times histories of joint angles [1] shown in table 1.
Each of the positions are in units of radians. For our model, we have 2 Degrees of Freedom (joint angles) that correspond to 2 joints, namely glenohumeral (GH) and elbow (EL). Each joint has a minimum and maximum value for position expressed as radians.

Forward dynamics
The Forward Dynamics Tool can drive a forward dynamic simulation. A forward dynamics simulation is the solution (integration) of the differential equations that define the dynamics of a musculoskeletal model 1,18 . By focusing on specific time intervals of interest, and by using different analyses, more detailed biomechanical data for the trial in question can be collected. From Newton's second law, accelerations (rate of change of velocities) of the coordinates can be described in terms of the inertia and application of forces on the skeleton as a set of rigid-bodies: This is given by the following equation Where q is the coordinate accelerations due to joint torques T Coriolis and centrifugal forces, C(q,q) as a function of coordinates q ,and their velocities q, gravity G(q), and other forces applied to

Neural Network Architecture
Neural Network architecture is imbibed from the computational phenomenon of the brain.
In the brain neurons are connected to each other and fire when the sufficiently appropriate stimuli are presented. Inputs gathered by the Neuron and then passed through an activation function which produces an output from the neuron. These outputs are then propagated throughout the neural network and in the final stage, combined into an output value used to interpret the relevant state of the application shown in fig 3. Neurons are arranged in layers labeled as input, hidden, or output. Input layers represent the inputs to the neural network. The hidden layers are between the inputs and outputs. The output layer is the observable. Each neuron in layer l receives weighted inputs from layer l-1and sums them together along with a bias before passing them through activation function. The bias shifts the input sums along the sigmoid function and can be beneficial to biasing the output towards a particular value. The activation function may change depending on the situation of the neural network, but typically it is the sigmoid function expressed in eq (3) and shown in fig 4. g(z) = 1/(1+e -z ) ... (3) Feed forward network with 1 input layer, 3 hidden layers, and a single output layer was used. It was fully connected such that each neuron in the previous layer is connected to all neurons in the next layer. Each hidden layer size is the same as the number of inputs created at the beginning of the network shown in fig 5.

Sampling
The sampling required motion kinematics data (.mot file) and state data in storage (Controls. sto) file. The states file can contain the time histories of model states, including joint angles, joint speeds, muscle activations, muscle fiber lengths, and more. The forward dynamics Tool to set the initial states of the model for forward integration. Muscle states can be estimated by solving for tendon and muscle fiber force equilibrium when the solution for equilibrium for actuator states is checked. The control data contains the time histories of the model controls (e.g., muscle excitations) to the muscles and/or joint torques, for this study biceps long head and short head was considered. It is possible to specify the controls as .sto files instead, with columns corresponding to desired excitations.

Fig. 4: Sigmoid function
Ci is the i th degrees of freedom. The max and min are the maximum values and minimum values of the degree of freedom, ej and is the j th element.

Initialisation
Initialization of the network was done with 13 random inputs each with 12 dimensions which correspond to the down sample rate of 3 for 100 time steps. The network has 5 layers where the first and the last are input and output layers, leaving 3 hidden layers. There are 50 training rounds in each of the 50 epochs. The neural network was tested by down sampling the motion files and controls file by a factor k = 3. This decreases the amount of neurons used for our network and significantly lowers the amount of time needed for training and testing. The training consists of each input which is a 1X3 vector which represents the state at some time step. Each weight is multiplied elemental wise with the input and subsequently summed together column-wise. The bias is then added to the result and then each element of the new vector is passed through the sigmoid function.

Training
The training involved modified back propagation algorithm along with stochastic gradient descent 10 . The error function is expressed in equation 5 n Arg min =1/2 (y i -w i ) 2 . where y i is the associated true output and w i is the predicted output from the neural network. The weights and biases of the network are represented by w and b. We sum across all the outputs of our network to gather all the contributing errors. By finding the gradients with respect to each parameter, we will then update each parameter by taking a step in the direction of the gradient to decrease this error function. The weights and biases are updated in the direction of the negative gradient of the performance function.
The algorithm for training is as follows: A random training example from training pool is chosen. Pass training input through the network and determine objective loss for Calculate each nodes in output layer by 1) Computing gradients with respect to biases 2) Computing gradients with respect to weights 3) Computing gradients with respect to previous layer end for (Calculate gradients with respect to outputs backwards through the layers) Repeat using output gradients from the next layer for: 1) Computing gradients with respect to biases 2) Computing gradients with respect to weights 3) Computing gradients with respect to previous layer until input layer is reached for Update all weights and biases based on the computed gradients by a step size end for For each neuron output: g is the sigmoid activation function. O l j is the output of node j in layer l where O i l-1 is the output from the previous layer and W and b are the weights and biases of the current layer respectively. The gradients for weight, biases and previous output are as follows: O i l-1 contributes to all of the outputs in the next layer and so in order to back propagate the contribution, we have to average all of its contributions as follows: The chain rule is applied to find the gradients with respect to the weights and biases of previous layers using the above equations. Once the gradients have been calculated, the weights were updated update each weight element by a step size into the gradient direction. The step size in the current study is 0.5 since the gradients were small compared to actual weights. This is probably due to the fact that normalization required each weight and input to be between 0 and 1 and thus these values will be multiplied together during training to result in a very small value. Therefore a step size large enough to influence the objective error function in each epoch was used.

Validation
In this study, 15 samples operating on 15 training examples for each of the 20 epochs were presented to the network. The network selects 95% of the data for training leaving 5% for validation The data produced by forward dynamics simulation from the 15 samples were validated. The network updates the weights and biases and uses those updated parameters for the validation set through training and is shown in fig 6.  In the past several models were developed to understand motor control in human that is incorporated in designing the robotic system. At present there are various models that use recurrent neural network along with supervised learning. A study by Kawato et al proposed, proposed a mathematical model which accounts for formation of hand trajectories. This model is formulated by defining an objective function, a measure of performance for any possible movement. This objective function is determined complex nonlinear dynamics of the musculoskeletal system square of the rate of change of torque integrated over the entire movement. This study was done in an effort to understand the trajectory planning and control of voluntary human movements 6,9 .

Results of Inverse
A more recent study on reaching tasks in monkeys was carried out by Susillo et al concluded the solution generated by recurrent neural network which was a low dimensional oscillator resulting in multiphasic commands; in order to understand the dynamics that transformed simple inputs into both spatially and temporally complex muscle activity 11 . The study by Jordan et al, demonstrated that certain classical problems associated with the notion of the "teacher" in supervised learning con be solved by judicious use of learned internal models as components of the adaptive system 10 . In particular it focused on how supervised learning algorithms can be utilized in cases in which an unknown dynamical system intervenes between actions and desired outcomes. Thus in the light of previous studies, the present study proposes a control algorithm based on feed forward that uses stochastic gradient descent to compute the errors for learning strategy. The study executes a forward dynamics simulation approximated by the feed forward neural network, in an attempt to understand the inputs from the CNS that are transformed into movements by the way of synergistic muscle activities. The network only updates the weights and biases through training and uses those updated parameters for the validation set. The validation graphs reflect the training error graphs in their behavior. The degree of freedom which represents the time index shows graph where the output times and input times are similar and that this error is the lowest error of all the degrees of freedom as would be expected.Further in the results it is seen that the validation error follows a similar pattern as the training error, suggesting the possibility approximating the forward simulation. This would indicate that forward dynamics simulation can be approximated by a neural network. It also important to considered the relatively small sample size used for the study as a limitation which resulted in similar patterns of the cumulative errors throughout the epochs. A larger sample size could have shown more pronounced results with regard to the changes in the errors values for the corresponding training errors with greater computational costs. Moreover, the study was circumscribed by the constraints in the degrees of freedom which limited the movements to elbow flexion in allowable range. The activities of the secondary muscles that aid in elbow flexion were consequently beyond the scope of this study.

CONCLUSION
A common cumulative error was achieved that corresponded to the training errors in order to approximate the simulation of forward dynamics. Hence it can be concluded that control algorithm that uses feed forward neural network along with back propogation by gradient descent is a robust method in an effort to mimic the execution and motor planning of the brain when performing simple tasks. Furthermore, the proposed control algorithm can be judiciously used with both supervised and unsupervised learning algorithms.