Neuromechanical Cost Functionals Governing Motor Control for Early Screening of Motor Disorders

1. Introduction

Primary motor cortex and supplementary motor areas govern the execution and planning of the voluntary motor actions of the skeletal muscles. The projections from primary motor cortex traverse through the corticobulbar and corticospinal tracts and synapse on to the lower motor neurons in the brain stem and spinal cord, respectively. The system of extrapyramidal tracts controls the movement through a non-pyramidal route composing of the nigrostriatal pathway, the basal ganglia, the cerebellum, the vestibular nuclei, and different sensory areas of the cerebral cortex (Hall and John, 2005). Reflexes on the other hand are the result of a hardwired mechanism where interneurons of spinal cord inhibit the stimulation of antagonistic muscles resulting in appropriate movement. Dysfunction of any of these areas can engender a motor deficiency (Kandel et al., 2000). Generally, damage to lower motor neurons produces local effects in the muscles such as weakness, muscle atrophy, and fasciculation. However, upper motor neuron damage can have diffused effects in the body-eliciting spasticity and babinski reflex in adults (Kandel et al., 2000). Quantification of the dynamics governed by the upper and lower motor neurons is a difficult problem.

Computationally, this problem is approached in two ways. First methodology is the use of detailed neural system models with varying complexity (integrate and fire neurons (Brunel and Hakim, 1999), spiking models (Izhikevich, 2003)) to learn the motor control and explain the motor behavior (Grillner et al., 2007; Eliasmith et al., 2012). These methodologies generally use reinforcement learning methods in conjunction, which derives an optimal policy under a set of constraints. Parameter estimation is difficult using these approaches due to (1) their large number and (2) difficulty in defining a tractable cost function. Even in the case where a parameter estimation is possible the models are still not in one to one correspondence with the physiology leading to abstractness in the result.

The second one involves defining an abstract functional, which is optimized to derive the control inputs for the mechanical system (Erdemir et al., 2007). This top-down approach to model the neural controller acts in tandem with the mechanical system that governs the trajectory of motion for a given control input. Therefore, quantifying the motor control involves quantification of two different systems - a mechanical system and the neural controller. These methods, even though abstract, can be solved perfectly to generate results that are optimal.

In this study, we use an approach closer to the second one. The neural controller is known to have characteristics of an optimizer where the cost functional minimized can be jerk, metabolic energy, etc. Therefore, in this study, we use inverse optimal control to find the cost functional that governs the motion. Inverse optimal control refers to a set of techniques that are used to learn the objective functional that governs the system using a data set (Mombaur et al., 2011). These techniques can identify underlying optimality criterion in human motions (Mombaur et al., 2011). In Berret et al. (2011), authors show evidence of composite cost functions governing the arm motion. But this work does not examine normal and patient population and uses numerical techniques that are computationally less efficient. Also, here the arm is modeled using two links rather than three links ignoring the interaction torques that is generated from the third link. The optimal control algorithm used in this model being non-linear also has an issue of non-uniqueness, which makes the applicability of the solution less likely. There is a need to develop a simpler but more efficient technique, so that, it can be used in clinical practice. The method also needs to be reproducible and solution of the optimal control algorithm unique making the results between different patients and at different time frames comparable. A typical reaching task in the horizontal plain neglects gravity. Therefore, the task needs to simulate a real life situation where the arm moves against gravity. This type of modeling although is not a replacement to a more detailed model of neural control, it is useful to obtain personalized neural models. This method overcomes the problems described in detailed models by abstracting away the parameters to be identified. The model-based parameter estimation techniques give insights into the governing motor dynamics and therefore can be used to replace or augment traditional scales used for the quantification. This model can be studied further to understand the behavior of the system against perturbations such as oscillations.

Oscillations are ubiquitous in the brain (Engel and Singer, 2001; Varela et al., 2001), which occur in different scales and at all levels of organization. The organizational levels can range from microscopic (Wang, 2010), mesoscopic (Nunez and Srinivasan, 2006; Cardin et al., 2009) to macroscopic (Llinas et al., 1998; Bollimunta et al., 2011) depending on their source of origin, that is, if oscillations arise from single neurons, a small groups of neurons, or different regions in the brain, respectively.

While intrinsic properties such as conductance govern the oscillations in cellular level (Llinas et al., 1991), the network connection strengths and excitatory-inhibitory nature of the connections control them at network level (Kilpatrick, 2013). A group of neurons engaging in oscillatory activity can be mathematically represented as a lumped system (Haken, 1996). Therefore, in this work, we study how oscillations can affect the motor control.

The primary objective of this study is the development of a methodology to determine a functional governing the motor commands in such a way that it quantifies the neuromotor control abstractly. This method shall give a unique solution for the optimal control problem analytically. This will help in developing an efficient solution that is comparable within and between patient groups at different time frames. The second objective is to validate this model in normal and patient population with mild neuromotor disabilities. (We chose the population with mild neuromotor disabilities as mild neuromotor disability is hard to quantify. The subtle changes in those patients can only be detected by an expert in the field.) The secondary objective is to test the model developed against perturbations and study its effects. We hypothesize that there will be a significant statistical difference between the normal and patient population when the parameters governing the optimal criterion are compared. We also hypothesize that the perturbations of at frequencies around beta frequencies will have higher energy as motor problems manifest around these frequencies.

2. Methodology

We present the detail methodology of deriving the cost functional governing the arm movement using inverse optimal control. The block diagram for the overall algorithm is given in Figure 1 The transformation matrix is analyzed to understand the principal direction of resistance and a time frequency analysis is done to study the effect of perturbations on the dynamical system. A flowchart is given in Figure 2 with various functional blocks relating with the equations presented in this work.

FIGURE 1

Figure 1. The overall algorithm: the Kinect data after resampling was used as the input to the algorithm. These data were then scaled to correct the model with the joint coordinates of the subject. This was then feed to an inverse optimal control algorithm where the functional J is optimized to match the Kinect data. The optimal control in itself involved an optimization of a functional using variational calculus formulation where the search is done in the space of smooth functions matching the boundary conditions.

FIGURE 2

Figure 2. Flowchart shows how the set of equations (1)–(30) are used. The equation numbers are given in brackets. The optimal control and model equations cannot be fully separated to have one causally influencing the other, therefore, it is shown as coupled. The last two equations come from the results section, which depicts the derived functional.

2.1. Optimal Control Formalism

Optimal control formalism relies on the variational formulation (Laub, 1979; Bhattacharyya et al., 2009) of the problem where a functional is minimized under differential constraints. This is particularly useful in situations where the cost has to be described as an integral of a function over a trajectory and the variables of the integrand are governed by system dynamics. A special case of this problem with linear differential constraints is described below.

The problem of minimizing a functional equation (1) with differential constraints equation (2) can be solved in the following way. In equation (2) only c(t) is unknown and is the control input, this need to be found in such a way that the cost functional J is minimized.

Let the functional to be minimized be

$\mathbf{\text{J}}={\int }_0^{\infty }\left({\mathbf{\text{x}}}^{\mathbf{\text{T}}}Z\mathbf{\text{x}}+{\mathbf{\text{c}}}^{\mathbf{\text{T}}}U\mathbf{\text{c}}\right)\mathit{dt},$ (1)

where Z and U are known and let the differential constraint be

$\frac{d\mathbf{\text{x}}}{\mathit{dt}}=\mathbf{\text{Ax}}\left(\mathbf{\text{t}}\right)+\mathbf{\text{Bc}}\left(\mathbf{\text{t}}\right)$ (2)

and x(t₀) = x₀ then the optimal

$\mathbf{\text{c}}\left(t\right)=-\mathbf{\text{Kx}}\left(t\right),$ (3)

where

$\mathbf{\text{K}}={\mathbf{\text{U}}}^{-1}{\mathbf{\text{B}}}^T\mathbf{\text{S}}$ (4)

and S is the solution of the following equation

$-\mathbf{\text{SA}}-{\mathbf{\text{A}}}^T\mathbf{\text{S}}+\mathbf{\text{SB}}{\mathbf{\text{U}}}^{-1}{\mathbf{\text{B}}}^T\mathbf{\text{S}}-\mathbf{\text{Z}}=0.$ (5)

This equation is known as algebraic Riccatti equation, which is easy to solve using standard techniques. For each J, there is a solution x that is optimal.

2.2. Mathematical Modeling and Cost Functional

A 3 linked arm model with 3-revolute joints was developed with a sigmoid non-linear constraint function (k ∈ C^∞)¹ imposed on the joints to avoid elbow hyper extension. This constraint is imposed so to restrict the motion in the physiological range. These were rigid links without friction and with shoulder joint as inertial frame of reference.

Control parameter c(t) are the torques at different joints; gravity is assumed to be acting at center of mass. Kane’s method (Kane et al., 1987) was used to derive the equations of motion (EOM) in the following form:

$\mathbf{\text{f}}+{\mathbf{\text{f}}}^{\ast }=0,$ (6)

where f is the generalized active force vector and f* is the generalized inertial force vector. These set of equations were then converted and assembled to the following form and solved using Hindmarsh (1983):

$\mathbf{\text{M}}\frac{d\mathbf{\text{x}}}{\mathit{dt}}=\mathbf{\text{r}},$ (7)

where M is the mass matrix, x is the state vector, and r is the forcing vector. Here, the generalized coordinates and speeds were chosen to be the angles and angular velocities, respectively. The system of equation (7) is non-linear. This was linearized at the unstable equilibrium point to convert to the following form. This equation was only used to derive the necessary controller c(t). c(t) and c are used to refer these control inputs interchangeably.

$\frac{d\mathbf{\text{x}}}{\mathit{dt}}=\mathbf{\text{Ax}}+\mathbf{\text{Bc}},$ (8)

where A ∈ ℝ^{6 × 6} and B ∈ ℝ^{6 × 3} represent² the linearized coefficients of the differential equation. The model was scaled for each subject to match the experimentally obtained lengths of the different joints. Equation (9) was then used as the cost functional. First term inside the integral represents cost of change in state x, and the second term represents the cost of controls c. The final time of the integral is assumed ∞ as the subjects were not given an end time to complete the task.

$\mathbf{\text{J}}={\int }_0^{\infty }\left({\mathbf{\text{x}}}^{\mathbf{\text{T}}}\mathbf{\text{Z}}x+{\mathbf{\text{c}}}^{\mathbf{\text{T}}}\mathbf{\text{Uc}}\right)\mathit{dt}={\int }_0^{\infty }{\mathbf{\text{I}}}_v+{\mathbf{\text{I}}}_c,$ (9)

where

${\mathbf{\text{I}}}_v={\mathbf{\text{x}}}^{\mathbf{\text{T}}}\mathbf{\text{Zx}},$ (10)

${\mathbf{\text{I}}}_v={\mathbf{\text{c}}}^{\mathbf{\text{T}}}\mathbf{\text{Uc}},$ (11)

$\mathbf{\text{x}}={\left[{\theta }_1,{\theta }_2,{\theta }_3,{\omega }_1,{\omega }_2,{\omega }_3\right]}^{\mathbf{\text{T}}},$ (12)

where θ₁–θ₃ are the joint angles in shoulder, elbow, and wrist, respectively, and ω₁–ω₃ are the angular velocities. x is the state variable, c is the control torques and Z ∈ ℝ^{6 × 6} and U ∈ ℝ^{3 × 3} are the weights of independent and controllable variables of the differential equation, respectively. These were constrained to be positive definite. As the function was linearized at the final state x*, and minimization of the cost function J also results in minimization of the Euclidean distance between this equilibrium point and the current state. Here, the matrices Z and U were assumed to be diagonal and arbitrary with the following structure for Z.

$\mathbf{\text{Z}}=\left[\begin{array}{cccccc}\hfill {\theta }_p\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\theta }_p\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\theta }_p\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\omega }_p\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\omega }_p\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\omega }_p\hfill \end{array}\right],$ (13)

where θ_p and ω_p represent penalty to angular displacements and velocities, respectively. We assumed U = 𝕀^{3 × 3} is the identity matrix making I_c fixed with respect to which I_v can be thought to be varied. In the inverse optimal control problem (Figure 1), a requirement for increased penalty to torques applied would translate to an appropriate change in the penalties θ_p and ω_p and therefore determination of the function is given by equation (14)

${\mathbf{\text{I}}}_v={\omega }_p\left({\omega }_1^2+{\omega }_2^2+{\omega }_3^2\right)+{\theta }_p\left({\theta }_1^2+{\theta }_2^2+{\theta }_3^2\right).$ (14)

Let Ω and Θ be mean squared values of ω₁₋₃ and θ₁₋₃, respectively.

${\mathrm{\Omega }}^2={\omega }_1^2+{\omega }_2^2+{\omega }_3^2,$ (15)

${\mathrm{\Theta }}^2={\theta }_1^2+{\theta }_2^2+{\theta }_3^2,$ (16)

$\Rightarrow {\mathbf{\text{I}}}_v={\omega }_p{\mathrm{\Omega }}^2+{\theta }_p{\mathrm{\Theta }}^2.$ (17)

Making a set of weights constant in relation to another also helps in reducing the search space for the algorithm thus reducing the computational cost considerably. The problem was formulated as an infinite horizon to account the fact that final time was not given to the patients as a criterion while performing the experiments. In summary, the cost functional J was solved with differential constraints in equation (8). The inverse optimal control problem of determining the matrix Z was done using Nelder and Mead (1965) and Wright (1996). The overall algorithm is depicted in Figure 1.

Mathematically, the optimization problem of finding the personalized J* can be described as in the following way. Each J maps to only one trajectory x; This makes x a function of J or f x(J)

${\mathbf{\text{J}}}^{\ast }=\underset{\mathbf{\text{J}}\in {\mathbf{\text{C}}}^2}{\text{argmin}}\mathbf{\text{\hspace{0.17em}Er}}\left(\mathbf{\text{x}}\left(\mathbf{\text{J}}\right),\widetilde{\mathbf{\text{x}}}\right),$ (18)

where x is the solution of the optimal control problem for a given J and $\widetilde{\mathbf{\text{x}}}$ is the measured state trajectory. The error function $\mathbf{\text{Er}}\left(\mathbf{\text{x}}\left(\mathbf{\text{J}}\right),\widetilde{\mathbf{\text{x}}}\right)$ maps the trajectories of the states $\widetilde{\mathbf{\text{x}}}$ and x to the l₂ norm of the difference between the corresponding displacements in the Cartesian y-coordinate of the wrist. The projection of the data in the way described helps to reduce the complexity of the optimization problem significantly in a way that is physically meaningful. The data collection procedure is described in Section 2.5.

2.3. Stiffness Ellipsoid or Principal Directions of Resistance in the State Space

The control torques τ depends on the states x in the following way:

$\mathbf{\tau }=\mathbf{\text{Kx}}.$ (19)

Let

$\mathbf{\text{K}}=\left[\begin{array}{cc}\hfill {\mathit{\text{K}}}_{\mathbf{\text{1}}}\hfill & \hfill {\mathit{\text{K}}}_{\mathbf{\text{2}}}\hfill \end{array}\right],$ (20)

$\mathbf{\text{x}}=\left[\begin{array}{c}\hfill {\mathit{\text{x}}}_{\mathbf{\text{1}}}\hfill \\ \hfill {\mathit{\text{x}}}_{\mathbf{\text{2}}}\hfill \end{array}\right],$ (21)

$\Rightarrow \mathbf{\tau }={\mathit{\text{K}}}_{\mathbf{\text{1}}}{\mathit{\text{x}}}_{\mathbf{\text{1}}}+{\mathit{\text{K}}}_{\mathbf{\text{2}}}{\mathit{\text{x}}}_{\mathbf{\text{2}}},$ (22)

where K ∈ ℝ^{3 × 6} is the transformation matrix, x₁ and x₁ are the angles and angular velocities, respectively. For constants K₁ and K₂ and x₁, the value of x₂ determines the value of applied torque. We define

${\mathbf{\tau }}^{\mathbf{\ast }}=\underset{{\mathit{\text{x}}}_{\mathbf{\text{2}}}\in {ℝ}^{1\times 3}}{\text{max}}\parallel \mathbf{\text{Kx}}\parallel ,$ (23)

where x₂ is the vector of angular velocities. The resistance torque applied by the neuromuscular system depends on two characteristics of the angular velocity vector x. First one, the magnitude of the angular velocity vector, which will directly increase the torque applied. The second one is the direction of the angular velocity vector. We denote this direction by ${\mathit{\text{x}}}_{\mathbf{\text{2}}}^{\mathbf{\ast }}$ and this is the principal eigenvector of the matrix K₂. If angular velocities are perturbed in this direction the resistance offered by the neuromechanical system peaks. This torque vector acts solely against the angular velocities resisting a higher velocity of motion. The largest component of this vector is indicative of the joint velocity that is resisted the most and the difference these eigen values among subjects represent intra-subject variability in resistance to motion.

2.4. Frequency Response and Time Frequency Analysis

To understand the response of the dynamical system behavior against perturbations, we studied the response of patient and normal systems to spectral perturbations. The new set of equations governing the perturbed system is same as that is described previously except equation (3), which is replaced by equation (24)

$\mathbf{\text{c}}\left(t\right)=-\mathbf{\text{K}}{\mathbf{\text{x}}}_n\left(t\right),$ (24)

${\mathbf{\text{x}}}_n=\mathbf{\text{x}}+\mathbf{\lambda }\mathit{cos}\left(\mathrm{\gamma }t\right)),$ (25)

$\mathbf{\lambda }=\mathrm{\kappa }\ast {\mathbb{I}}_{6\times 6},$ (26)

where x_n is the “perturbed” state signal, which represents different oscillations that arise in the brain and 𝕀_{6 × 6} is an identity matrix. The resulting solution of the optimal dynamical system for different γ values was computed, and a time frequency analysis of the solutions was performed using equation (27). κ (κ = 0.2) was chosen to be of a value that produces physiological relevant ranges for the results and easy to manage numerically.

$W\left(t,f\right)=\int \mathbf{\text{x}}\left(t+\mathrm{\tau }∕2\right){\mathbf{\text{x}}}^{\prime }\left(t-\mathrm{\tau }∕2\right)e^{-j2\pi \mathrm{\tau }f}d\mathrm{\tau }.$ (27)

The function W (t, f) corresponding to each solution x is analyzed the following way by using equation (28) to understand the time frequency pattern against different frequency of inputs. Here, τ is an arbitrary delay variable.

$E_{t,f}={\int }_0^f{\int }_0^{t}W\left(t,f\right).$ (28)

The E_t,f was computed for a time window of 60 s for all frequencies estimated. This limit was chosen to understand the effect of the frequencies during the motion rather than at the trailing end of the bell curve of velocities.

2.5. Data Collection Using Kinect and Signal Processing

The range of motion for left shoulder flexion is collected using Microsoft Kinect XBox One.³ A total of 13 normal subjects and 19 patients with neuromuscular disorders were selected for the study. Each session was repeated four times for all normal subjects resulting in total of 52 samples whereas a single trial was used for patients generating total 19 samples.

The clearance on ethical issues corresponding to the patient data collection has been obtained from the Institute Ethics Committee (IEC) in Advanced Medical Research Institute (AMRI) Hospitals. Informed consents are also taken from the subjects in written form. The data are annonymized by representatives of AMRI in a community setup and provided to TCS Research Lab for purpose of the analysis presented in this paper. The clearance on ethical issues corresponding to the handling and analysis of the data collection has been obtained from relevant body in TCS. The patients were of the age group range 56–74 with a mean of 65.6 and a SD of 4.62. The normal subjects were of the age group 22–45 with an SD of 10.3 and mean of 30.1. For the healthy group male, the average height was 172 (163–183) cm, and for the females, it was 158 (150–161) cm. For male patients, the average height was 171 (164–180) cm, and for the female patients, it was 160 (156–170) cm. The average weight of the healthy male group was 78.25 (70–84) kg, and for the female group, it was 63.5 (56–68) kg. The average weight of the patient male group was 74.83 (62–90) kg and that of the female patients was 67.2 (60–74) kg. The patient group had mild neuromotor abnormalities due to diabetes and old age.

These data were then resampled from 30 to 60 fps for the analysis using built-in python Fourier method. This resampling was done to make the numerically generated trajectory vector, and the data collected have the same dimension. The lengths of different bone segments were then computed for each subject using the 3D coordinates of the skeleton data of the initial 5 s of static window before the start of the exercise. This was used for scaling, and the scaled model was given as input to the inverse optimal control algorithm to optimize the cost functional.

3. Results

We observed that the displacements obtained from the model matched with the experimental results from the Kinect with negligible error, see Figure 3. In Figure 3, each subject is represented using a single color. The subfigure of Figure 3 titled Model vs Data shows how the modeled displacement is matched to the measured displacements (noisy due to the measurement noise) using the optimization procedure given in Subsection 2.2 using equation (18). The computed velocity is also shown to match the measured values as shown in Figure 4. This verifies the ability of the model to approximately reproduce the kinematics of the subjects.

FIGURE 3

Figure 3. Model predicted torques for wrist, shoulder, and elbow joints for 4 normal subject samples are shown (each of the subjects is assigned a single color in the graph). The model is personalized to match the experimental observations of displacement. The range of motion is very well matched by the model used.

FIGURE 4

Figure 4. The modeled velocity profiles match qualitatively with the Kinect data. The deviation of the Kinect data from the model is due to the fact that the process of finding the instantaneous velocity involves taking numerical derivatives, which amplify the high-frequency noise. Modeled velocity is in blue, and the data are shown in green.

Figure 5 shows the comparison of the θ_p and ω_p in normal and patient categories. To analyze the significance, a Welch’s t-test (Welch, 1947) was performed without assuming equal variance or equal number of samples (Ruxton, 2006). We found these cost-function parameters vary significantly between normal subjects and patients (see Figure 5) with a test statistic value of −2.28, p-value of 0.027 and test statistic value of 2.3, p-value of 0.023 for ω_p and θ_p, respectively.

FIGURE 5

Figure 5. Variation in θ_p and ω_p for normal and patient categories. The θ_p (black) is significantly higher in patients than in normals whereas ω_p (red) is significantly lower in patients than in normals.

The test used did not assume the equal variance, which if assumed the p-value will be lower. Therefore, these parameters can be considered as a good measure to understand the underlying optimality criterion for the motor control mechanism of a patient. A lower ω_p was observed in patients indicating jerky and fast movements. On the other hand, a higher θ_p change was observed in patients indicating a lower ability to hold the hand at intermediate position. There are two objectives for this study, first one being the detection of optimal functions governing the neuromotor control and the second one finding the response of the model against perturbations.

3.1. Optimal Functions—Patient and Normal

Here, we detail the optimal functions of the normals and patient groups. Using the optimal control equations described in the methods sections, we estimated the functions I_v for the normals (I_vn) and the patients (I_vp), which take the forms described in equations (29) and (30).

${\mathbf{\text{I}}}_{\mathit{vn}}=9,315,830{\mathrm{\Omega }}^2+12,926,292{\mathrm{\Theta }}^2,$ (29)

${\mathbf{\text{I}}}_{\mathit{vp}}=5,245,197{\mathrm{\Omega }}^2+13,950,617{\mathrm{\Theta }}^2.$ (30)

The contour plots of I_vp and I_vn with respect to different values of Θ and Ω are shown in Figure 6. The contours in the case of normals are more circular than that of the patients. This is due to the fact that ratio of ω_p to θ_p is lower in the case of patients than in the case of normals. This increases the tendency to move faster in the case of the patients leading to a more jerkier motion and inability of the patients to slowdown in an intermediate position. Figure 7 shows how in different subject groups the θ_p and ω_p varies in detail. Here, each point represents a subject. The value of θ_p is higher in the case of patients (shown in red) and lower in the case of normals (shown in black). Figure 7 also shows the pattern of θ_p and ω_p with respect to the average velocity of motion. An increased θ_p and a lowered ω_p show similar effects on the average velocity. The optimal torques for both patient and normal cases are shown in Figure 8. As expected from the cost functions, the normals use lesser torques than that of the patients for reaching the same targets. The shoulder experiences maximum range of torque inputs compared with the elbow and wrist. This higher range of torques in the shoulder helps it resist the higher moments against which it moves.

FIGURE 6

Figure 6. Shows the contour plot of I_v(Θ, Ω) of normals and patients with abscissa show change in Ω and ordinate change in Θ. The normal I_v is more circular compared with the patient I_v. This owes to the fact that the motion is more smoother and of lower velocity in normals compared with the patient group.

FIGURE 7

Figure 7. The variation of parameters θ_p and ω_p with respect to the average velocity is shown. The blue curve is a fit that predicts the value of parameters for a given value of average velocity. The values of θ_p are higher for the patients and lower for the normals; ω_p shows an opposite trend.

FIGURE 8

Figure 8. Torques of the normal (black) and patient (red) are shown here for shoulder, elbow, and wrist. The patients tend to use a higher torque than normals, and shoulder torque has a higher range in compared with elbow and wrist.

3.2. Response to Oscillations and the Eigen Values

We then studied the response of the system. The time frequency analysis performed using equation (27) showed and increased energy in earlier in time frame than later as shown in Figure 9. The energy levels, computed using equation (28), are higher at lower frequencies and around the frequencies that correspond to the beta oscillations. To illustrate this, we have computed E_t,f for t = 60 and f = max (as there is a bound to the maximum frequency content that can be estimated) was computed, normalized by dividing with the maximum value of the response and plotted in Figure 10.

FIGURE 9

Figure 9. The WVD of the velocity waveforms is shown. The patient group tends have a higher energy at the start of the motion compared with the normals.

FIGURE 10

Figure 10. The value of ${\int }_{\text{0}}^{\text{30}}{\int }_{\text{0}}^{\text{60}}W\left(t,f\right)$ for different frequencies is shown. At low frequencies and close to beta oscillation frequencies, the response of the system is observed to be higher.

Even though, all the principle directions remained the same for patient and normal groups, normals had lower eigen values compared with the patients. This shows why normal population apply lesser torques compared with the patient. The eigen directions ${\mathit{\text{x}}}_{\mathbf{\text{2}}}^{\mathbf{\ast }}$ remained the same for the normals and patients. The wrist velocities corresponding to the input oscillation at frequencies 0.5–45.5 Hz are shown in Figure 11.

FIGURE 11

Figure 11. Velocity profiles of motion corresponding to different frequency of inputs (0.5, 15.5, 30.5, and 45.5 Hz left to right) are shown. The higher frequencies do not affect the velocity profiles as much as the lower frequency values. Also there is a qualitative change in the velocity profiles at 30.5 Hz. The normal subjects are shown in black, and the patients are shown in red.

4. Discussion and Future Research

The objective of this work was to develop a methodology to quantify and understand the neuromechanical pathologies and validate it using the methodology on two subject groups one normal and another with mild neuromotor abnormalities. We found that the inverse optimal control methodology can be used to personalize model of the arm dynamics. Thus, the cost functional J parameterized by θ_p and ω_p can be used to understand subject specific variations. We also have showed how the eigen vectors of the matrix K₂ controls the torque applied. The different eigen values controlling the resistant torques are indicative of the higher resistance offered by the normals against a drastic change in velocity. This protects the population against natural wear and tear. In Figure 6, for the patient case moving through the ordinate will reach the maximum values (shown in red) compared with moving through the abscissa. We also note that in the case of normals as the cost function is more steeper ease of reaching the minima is more faster. It can be speculated that in a noisy case, this will govern the ease with which motor control decisions are taken making the control in case of the patients more difficult.

The lower frequencies and the frequencies around 30 Hz show an elevated response eliciting the differential response of the controller against different frequencies. This also indicates why some oscillations in the brain can result in motor symptoms.

In future, we intend to collect data from patients with neuromuscular disorders and estimate the parameters as the disease progresses. We argue that reversal of the parameters from diseased states to the normal (see Figure 5) is a necessary condition to validate treatment efficacy (from elliptical to normal). The cost functional can be further generalized, and a detailed set of parameters can be estimated to better understand the physiology. Careful considerations are needed from optimization front while generalizing the cost functional as this will increase the dimensionality of the problem.

Even though we have used a simple planar 3-link arm that is not a major limitation of the study as the task was constrained to a planar motion. But constraining the task to a planar motion may not be able to reveal all aspects of proprioception and motor control but this helps in making the task easily reproducible in a clinical setting. Another limitation is the linearization that was carried out to make the solution unique for the optimal control problem. Linearizing at different points, solving the optimal control problem and comparing the results will resolve this limitation. Extending the model to a realistic motion involves adding more dimensionalities to shoulder, elbow, and wrist joints and constraining the motion realistically. Although, this not necessary to simulate a planar motion could be used for other tasks such as pick and place. But this will render the problem of optimization computationally expensive. But simulating and doing the optimal control on such a task may provide very valuable information to the clinician such as how different muscles are affected by the neuromotor disorders. This will help in designing a task to strengthen the appropriate muscles. The shoulder flexion task involves movement of the arm against gravity. Most of the practical real life scenarios such as combing hair involve movement against gravity. θ_p and ω_p are abstract quantities. Supplementary motor area and other areas of the brain involved in planning of the motion directly affect the magnitude of these quantities. Velocity-dependent conditions such as hyper reflexia and spasticity can directly affect the eigen values of the matrix K₂ and thus resistance against velocities.

More precisely, an ischemic stroke in the middle cerebral artery (MCA) affecting the motor cortex or sensory cortex will have an impact on the parameters mentioned. The parameters can also be affected by stroke in subcortical regions which in turn result in motor deficits. In addition, the parameters could be influenced by aging-related changes such as changes in spindle innervation, increased co-contraction of agonist and antagonist muscles, and decreased reaction times due to decreased motor conduction. Precisely, mapping different disease conditions to parameters is a future work.

Noise in smaller levels will not affect the computation of these parameters as the differential equation can only give smooth solutions, the method will act as a filter. But in cases where the noise is very high and skewed will contribute to the study results. A careful study may be conducted later to understand the effects of noise in the analysis.

Ethics Statement

The clearance on ethical issues corresponding to the patient data collection has been obtained from the Institute Ethics Committee (IEC) in Advanced Medical Research Institute (AMRI) Hospitals. Informed consents are also taken from the subjects. The data are annonymized by representatives of AMRI in a community setup and provided to TCS Research Lab for purpose of the analysis presented in this article. The clearance on ethical issues corresponding to the handling and analysis of the data collection has been obtained from relevant body in TCS.

Author Contributions

MU conceived the model, developed it, analyzed the data, and prepared the manuscript. AS conceived the model and prepared the manuscript. KC, DC, and AD designed the experiment and collected the data.

Conflict of Interest Statement

MU, AS, KC, and DC are employed as researchers by TCS Research and Innovation, Tata Consultancy Services Ltd. AD is the Director of Neurorehabilitation Senior Consultant Neurologist, AMRI Institute of Neurosciences Kolkata. He received honorarium as a research advisor from TCS Research and Innovation, Tata Consultancy Services Ltd. All the authors declare no other competing interest.

Funding

This work was done in TCS Research & Innovation, funded by Tata Consultancy Services Limited.

Footnotes

1. ^Cⁿ is the set of all n times differentiable functions.
2. ^ℝ^m×n is the set of all m × n matrices.
3. ^http://www.xbox.com/en-IN/xbox-one/accessories/kinect.

References

Berret, B., Chiovetto, E., Nori, F., and Pozzo, T. (2011). Evidence for composite cost functions in arm movement planning: an inverse optimal control approach. PLoS Comput. Biol. 7:e1002183. doi: 10.1371/journal.pcbi.1002183

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharyya, S. P., Datta, A., and Keel, L. H. (2009). Linear Control Theory: Structure, Robustness, and Optimization. Boca Raton, FL: CRC Press, 33.

Google Scholar

Bollimunta, A., Mo, J., Schroeder, C. E., and Ding, M. (2011). Neuronal mechanisms and attentional modulation of corticothalamic alpha oscillations. J. Neurosci. 31, 4935–4943. doi:10.1523/JNEUROSCI.5580-10.2011

CrossRef Full Text | Google Scholar

Brunel, N., and Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11, 1621–1671. doi:10.1162/089976699300016179

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardin, J. A., Carlén, M., Meletis, K., Knoblich, U., Zhang, F., Deisseroth, K., et al. (2009). Driving fast-spiking cells induces gamma rhythm and controls sensory responses. Nature 459, 663–667. doi:10.1038/nature08002

PubMed Abstract | CrossRef Full Text | Google Scholar

Eliasmith, C., Stewart, T. C., Choo, X., Bekolay, T., DeWolf, T., Tang, Y., et al. (2012). A large-scale model of the functioning brain. Science 338, 1202–1205. doi:10.1126/science.1225266

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., and Singer, W. (2001). Temporal binding and the neural correlates of sensory awareness. Trends Cogn. Sci. 5, 16–25. doi:10.1016/S1364-6613(00)01568-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Erdemir, A., McLean, S., Herzog, W., and van den Bogert, A. J. (2007). Model-based estimation of muscle forces exerted during movements. Clin. Biomech. 22, 131–154. doi:10.1016/j.clinbiomech.2006.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S., Kozlov, A., Dario, P., Stefanini, C., Menciassi, A., Lansner, A., et al. (2007). Modeling a vertebrate motor system: pattern generation, steering and control of body orientation. Prog. Brain Res. 165, 221–234. doi:10.1016/S0079-6123(06)65014-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Haken, H. (1996). Principles of Brain Functioning: A Synergetic Approach to Brain Activity, Behaviour and Cognition. Berlin, Heidelberg: Springer.

Google Scholar

Hall, A. C. G., and John, E. (2005). Textbook of Medical Physiology. Technical Report, ISBN 978-0-7216-0240-0. Philadelphia: WB Saunders.

Google Scholar

Hindmarsh, A. C. (1983). “Odepack, a systematized collection of ode solvers,” in Scientific Computing, Vol. 1, eds R. S. Stepleman, et al. (North-Holland, Amsterdam: IMACS Transactions on Scientific Computation), 55–64.

Google Scholar

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural. Netw. 14, 1569–1572. doi:10.1109/TNN.2003.820440

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandel, E. R., Schwartz, J. H., Jessell, T. M., Siegelbaum, S. A., and Hudspeth, A. J. (2000). Principles of Neural Science. New York: McGraw-Hill, 4.

Google Scholar

Kane, T., Ryan, R., and Banerjee, A. (1987). Dynamics of a cantilever beam attached to a moving base. J. Guid. Control Dyn. 10, 139–151. doi:10.2514/3.20195

CrossRef Full Text | Google Scholar

Laub, A. (1979). A schur method for solving algebraic Riccati equations. IEEE Trans. Automat. Contr. 24, 913–921. doi:10.1109/TAC.1979.1102178

CrossRef Full Text | Google Scholar

Llinas, R., Ribary, U., Contreras, D., and Pedroarena, C. (1998). The neuronal basis for consciousness. Philos. Trans. R. Soc. Lond. B Biol. Sci. 353, 1841–1849. doi:10.1098/rstb.1998.0336

CrossRef Full Text | Google Scholar

Llinas, R. R., Grace, A. A., and Yarom, Y. (1991). In vitro neurons in mammalian cortical layer 4 exhibit intrinsic oscillatory activity in the 10-to 50-hz frequency range. Proc. Natl. Acad. Sci. U.S.A. 88, 897–901. doi:10.1073/pnas.88.3.897

PubMed Abstract | CrossRef Full Text | Google Scholar

Kilpatrick, Z. P. (2013). “Wilson-cowan model,” in Encyclopedia of Computational Neuroscience, eds Jaeger, Dieter, Jung and Ranu (New York, NY: Springer New York), 1–5. doi:10.1007/978-1-4614-7320-6_80-1

CrossRef Full Text | Google Scholar

Mombaur, K., Laumond, J.-P., and Truong, A. (2011). “An inverse optimal control approach to human motion modeling,” in Robotics Research, eds P. Cédric and H. SiegwartGerhard (Berlin, Heidelberg: Springer), 451–468.

Google Scholar

Nelder, J. A., and Mead, R. (1965). A simplex method for function minimization. Comput. J. 7, 308–313. doi:10.1093/comjnl/7.4.308

CrossRef Full Text | Google Scholar

Nunez, P. L., and Srinivasan, R. (2006). Electric Fields of the Brain: The Neurophysics of EEG. USA: Oxford University Press.

Google Scholar

Ruxton, G. D. (2006). The unequal variance t-test is an underused alternative to student’s t-test and the Mann–Whitney u test. Behav. Ecol. 17, 688–690. doi:10.1093/beheco/ark016

CrossRef Full Text | Google Scholar

Varela, F., Lachaux, J.-P., Rodriguez, E., and Martinerie, J. (2001). The brainweb: phase synchronization and large-scale integration. Nat. Rev. Neurosci. 2, 229–239. doi:10.1038/35067550

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J. (2010). Neurophysiological and computational principles of cortical rhythms in cognition. Physiol. Rev. 90, 1195–1268. doi:10.1152/physrev.00035.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Welch, B. L. (1947). The generalization of student’s’ problem when several different population variances are involved. Biometrika 34, 28–35. doi:10.2307/2332510

CrossRef Full Text | Google Scholar

Wright, M. H. (1996). Direct search methods: once scorned, now respectable. Pitman Res. Notes Math. Ser. 344, 191–208.

Google Scholar