University of Southern California Dissertations and Theses

Model-based studies of control strategies for noisy, redundant musculoskeletal systems

(USC Thesis Other)

Model-based studies of control strategies for noisy, redundant musculoskeletal systems

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MODEL-BASED STUDIES OF CONTROL STRATEGIES FOR NOISY, REDUNDANT MUSCULOSKELETAL SYSTEMS by Dan Song A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (BIOMEDICAL ENGINEERING) August 2008 Copyright 2008 Dan Song ii DEDICATION To my dearest parents – JINGLING XIE & XUEFENG SONG iii ACKNOWLEDGEMENTS I want thank all the people who made this work possible and helped me go through this long, intensive, sometimes even tearful, but most of time happy, challenging and very fruitful process. Among all the people, I would like to especially thank the special members of my Ph.D. Committee – Professor Gerald E. Loeb, Professor Francisco J. Valero-Cuevas, Professor Jill McNitt-Gray, Professor Ning Lan, and Professor James Gordon – for their professional guidance and technical support. Particularly, I would like to greatly thank Professor Loeb, my primary advisor, for providing your constant support, generous guidance, intelligent advices and academic challenge to help me grow into an independent thinker and researcher. I learned from you a way of critical thinking not only as an engineer, but as a life scientist that is responsible to discover the truth of “motor nature”. Your engagement into my thesis work is really appreciated, and working together with you has taught me a lot more than I can ever obtain from textbook. You would be my lifelong mentor who can lead me and guide me in every step of my future development. I also thank Dr. Lan, my secondary advisor, for your steady-fast academic and financial support in the first three years of my research work, and for giving me many chances to attend conferences and publicize this precious work. I am very grateful to Professor Gordon, in that, you not only have introduced me into the promising project, and provided me with generous instruction and advices, but whenever I was in a conflicting mind or downturn state, you have always been there to encourage me and provide timely help and support. Many thanks to Professor Valero-Cuevas for being in my committee and guiding me in this last year of the Ph.D. program. Although only for a very short period, working with you has greatly enlightened my mind, and brought many refreshing and inspiring ideas into my research. Dr. McNitt-Gray is the first professor I talked to at USC. I really thank you for bringing me into this wonderful school and this excellent program, and for offering your intelligent advice and guidance whenever I need them. I also want to express my gratitude to Professor Rahman Davoodi for the technical support, and Professor Michael Khoo for providing financial support in the last year of my Ph.D. study. The VA model development was based upon work supported by the National Science Foundation under Grant IOB 0352117, and VM work and its on-line release was supported by Alfred E. Mann Institute for Biomedical Engineering at the University of Southern California. iv I also want to thank my previous advisor, Professor Ning Pan from the department of Biological System Engineering at University of California, Davis (UCD), for introducing me to the field of biomechanics and movement control. Although I did not finish my Ph.D. at UCD, you are still my life long mentor. I also thank Professor Maury L. Hull from the department of Mechanical and Aeronautical Engineering at UCD. The research I conducted under your guidance has helped me make my mind in pursuing the field of human motor control. Many thanks would go to Dr. Veronica J. Santos for your generous help and sharing your research ideas and thoughts with me both as a friend, a colleague and a mentor. And Dr. Weiwei Li, thank you for helping me in understanding the optimal feedback control theory, which have broadened my view in robotics and control theory. My sincere gratitude to many good friends, Yi Fan, Zhou Lei, Ming Suvimol, Limei Cheng, Christina Lin, Li Wen, Jiongjiong Xu, Xin Chen, Hongjuan Jiao, Mingyuan Lv, Juha Kaarle Pietari Heiskala, Gunther Krausz, Manish Kurse, Giby Raphael, Markus Hauschild, Younggeun Choi, Feng Qi, Choel Han, Jay Mung, Jeong-Yoon Lee, Djordje Popovic and his wife Karolina, Nicholas Wettels, Rahul Kaliki … I can always count on your support, even sometimes far away from another side of the ocean; just a brief phone call or a warm email have greatly lifted my spirit. As I mentioned, there are too many names to be thanked in this short page. But I can never forget my ever-lasting, unconditionally supportive, dearest parents. Thanks! mom and dad! For your love, patience, guidance, support; for giving me life, for letting me grow safely under your wings, for opening the entire world in front me, and for encouraging me to fly myself to explore it! Father, you have sacrificed yourself for me, quietly, lovingly, without any complaints. Without you, I would not have come to the Ph.D. program at the first place, and now you can not even see me in this beautiful gown and hood. I dedicate this dissertation to you for your happy viewing in another world. Mother, you provided me with not only love, but also precious guidance and unquestionable trust that has been the most incentive for every step forward in my life. Now you are here, helping me, caring me, encouraging me at the last few crazy months. You ensure me, once again, that I can achieve this goal of my life. ☺ - Dan Song v TABLE OF CONTENTS DEDICATION ........................................................................................................................................ ii ACKNOWLEDGEMENTS ................................................................................................................... iii LIST OF TABLES ................................................................................................................................ vii LIST OF FIGURES.............................................................................................................................. viii ABSTRACT......................................................................................................................................... xiii CHAPTER I: INTRODUCTION .............................................................................................................1 i. ORGANIZATION OF THE THESIS........................................................................................2 ii. MOTOR REDUNDANCY ........................................................................................................4 a) What is motor redundancy?...................................................................................................4 b) How to resolve motor redundancy?.......................................................................................5 iii. MOTOR VARIABILITY AND MOTOR NOISE.....................................................................8 a) Where does the movement variability come from?...............................................................9 b) What determines the nature of signal dependent motor noise? ...........................................12 c) What does the noise and redundancy imply about movement control?...............................15 iv. OPTIMALITY PRINCIPLE IN MOTOR CONTROL............................................................21 a) What is and why need optimization?...................................................................................22 b) What constitute performance criteria?.................................................................................23 c) Open-loop and closed-loop control .....................................................................................28 d) Further considerations in optimal control............................................................................30 v. IMPEDANCE CONTROL IN BIOLOGICAL MOTOR SYSTEM ........................................31 a) What is and why need impedance control? .........................................................................33 b) Has impedance control scheme been seen in human movements?......................................36 c) How do impedance control and signal-dependent noise interplay? ....................................38 d) In view of optimal impedance control.................................................................................40 vi. REALISIM IN SYSTEM-LEVEL MODEL AND SIMULATED TASK...............................41 vii. SIGNIFICANCE .................................................................................................................44 CHAPTER II: VIRTUAL MUSCLE MODEL ......................................................................................45 i. ABSTRACT.............................................................................................................................45 ii. INTRODUCTION ...................................................................................................................45 iii. MATERIAL AND METHODS ...............................................................................................47 a) Musculotendon dynamics....................................................................................................47 b) Recruitment types................................................................................................................48 c) Contractile dynamics...........................................................................................................51 d) Validation............................................................................................................................56 iv. RESULTS ................................................................................................................................56 a) Stability ...............................................................................................................................56 b) Accuracy for static conditions.............................................................................................57 c) Accuracy for dynamic conditions........................................................................................58 d) Simulation speed .................................................................................................................60 v. DISCUSSION..........................................................................................................................60 vi. CONCLUSION........................................................................................................................62 CHAPTER III: SENSORIMOTOR SYSTEMS MODEL OF HUMAN ARM......................................63 i. ABSTRACT.............................................................................................................................63 ii. INTRODUCTION ...................................................................................................................63 iii. MATERIAL AND METHODS ...............................................................................................65 vi a) Musculoskeletal model........................................................................................................66 b) Virtual Muscle model parameters .......................................................................................70 c) Implementation of proprioceptor models ............................................................................73 d) Model integration in Matlab/Simulink ................................................................................75 e) Dynamic simulation using the VA model ...........................................................................76 iv. RESULTS ................................................................................................................................77 a) Parameterizing the VA model .............................................................................................77 b) Dynamic validation of the VA model .................................................................................83 v. DISCUSSION..........................................................................................................................85 vi. CONCLUSION........................................................................................................................88 CHAPTER IV: CONTROL STRATEGIES FOR SIGNAL DEPENDENT MOTOR NOISE ..............89 i. ABSTRACT.............................................................................................................................89 ii. INTRODUCTION ...................................................................................................................90 iii. MATERIAL AND METHODS ...............................................................................................94 a) Methods overview...............................................................................................................94 b) Musculoskeletal model and posture tasks ...........................................................................96 c) Signal-dependent noise and force variability in Virtual Muscle .........................................97 d) Co-contraction in a single-joint system for postural accuracy ...........................................98 e) Monte Carlo dynamic simulations ......................................................................................98 f) Postural variability versus energy consumption................................................................100 g) Control of postural variability along desired axis of accuracy ..........................................101 iv. RESULTS ..............................................................................................................................102 a) Natural Continuous Virtual Muscle simulates realistic force variability...........................102 b) Co-contraction can suppress SDN for single-joint postural accuracy ...............................104 c) CNS coordinate redundant muscles for multi-joint postural accuracy ..............................105 d) Muscle redundancy allows task-relevant control of SDN .................................................108 v. DISCUSSION........................................................................................................................115 vi. CONCLUSION......................................................................................................................122 CHAPTER V: SUMMARY AND FUTURE WORK..........................................................................123 i. SUMMARY...........................................................................................................................123 ii. FUTURE WORK...................................................................................................................124 a) Generalization to 3-D musculoskeletal systems................................................................124 b) Experimental testing of optimal control hypotheses .........................................................125 BIBLIOGROPHY ................................................................................................................................131 vii LIST OF TABLES Table II.1: Summary of model equations and best-fit constants ............................................................53 Table II.2: Symbols and definition of VM parameters...........................................................................55 Table III.1: Muscle origin (O) and insertion (I) points...........................................................................69 Table III.2: Segment mass and inertia parameters..................................................................................70 Table III.3: Muscle architectural parameters..........................................................................................71 Table III.4: Muscle fiber type fractional distribution and number of motor units..................................72 Table III.5: Comparison of muscle peak forces with literature values...................................................80 Table III.6: Comparison of optimal fascicle lengths with literature values............................................81 Table III.7: Comparison of tendon slack length with literature values...................................................82 viii LIST OF FIGURES Figure II.1: Mechanical structure of a single and lumped unit in the VM model. The contractile element (CE) operates in parallel with the passive element (PE), which consists of stretching (PE1) and compressing components (PE2), to represent the fascicles. The fascicle element is in series with muscle mass element and series-elastic element (SE), which represent the combined tendon and aponeurosis. F pe1 results from the well recognized non-linear spring K pe1 with a viscosity element ( η pe1 ) that resists quick stretch and compression in the passive muscle, while F pe2 results from a non-linear spring K pe2 resisting compression at the thick myofilaments during active contraction at short lengths, thus is proportional to the activity of the muscle; F se results from a non-linear spring K se with a low stiffness toe region. The F ce’ is the force produced by combined contractile and passive components in the fascicle or contractile element; the difference of F ce’ and F se operate on the muscle mass to drive the muscle contraction dynamics........................................................................................................47 Figure II.2: Continuous and discrete recruitment algorithms in Virtual Muscle model. (a) Natural Discrete recruitment algorithm as applied to a muscle consisting of 3 simulated slow-twitch and 3 fast-twitch motor units, respectively. U th i is the recruitment threshold of i th motor unit; U r is the activation level at which all the motor units are recruited. Once a motor unit is recruited, the firing frequency of the unit will rise linearly with U between f min and f max . This recruitment scheme mimics biologic recruitment of motor neurons. (b) Natural Continuous recruitment algorithm as applied to two fiber types, slow and fast. U th i is the recruitment threshold of i th fiber type or unit; the first fiber type U th 1 =U th slow is fixed at 0.001. Similar to the discrete algorithm, once a fiber type is recruited, the firing frequency of the fiber will rise linearly with U between f min and f max . .....................................................................................48 Figure II.3: Effects of L eff on musculoskeletal dynamics. Simulations in (a) and (b) are actuated by discretely recruited muscle models with and without L eff . (a) A demonstration of the instability caused by L eff in Virtual Muscle through an open-loop dynamic simulation using a 2-DOF-6-Muscle arm model. The muscle activation commands of the antagonist muscles across shoulder (sh), elbow (el) and both joints (bi): U sh =0.9, U el =0.4, U bi =0.0. (b) A demonstration of stable open-loop response at the shoulder and elbow joints after removal of the effective length state. ..............................................................................................................57 Figure II.4: Comparison of steady-state force-activation relationship of continuous and discrete recruitment schemes. Deltoid posterior was at L mt =16 cm between the Natural Continuous algorithm and Natural Discrete models with a small motor neuron pool (3 slow and 3 fast motor units) and a large motor neuron pool (30 slow and 30 fast motor units). All the models have 60% slow fibers and 40% fast fibers. (a) F-U curves spanning entire activation rage (U=0-1); (b) F-U curves zoomed in at the low activation levels (U=0-0.3). ................................58 Figure II.5: Comparison of dynamic responses between the Natural Continuous muscle model and Natural Discrete muscle model with large motor neuron pool (30 slow and 30 fast motor units). Both models have 60% slow fibers and 40% fast fibers. (c) shows the dynamic force responses given the muscle inputs depicted in (a), sinusoidal activation (U: 0-1 at 3 Hz) and constant length (L mt = 16 cm); (d) shows the dynamic force responses given the muscle inputs depicted in (b), sinusoidal muscle length (L mt : 14-18 at 3 Hz) and constant activation (U=0.5). (e) shows the ratios of mean forces and force modulation ranges between the continuous and discrete models over the physiological range of modulation of activation U (1-10 Hz).......................................................................................................................................59 ix Figure II.6: Comparison of dynamic force responses between the Natural Discrete muscle models with and without effective lengths (L eff ). The models include a large motor neuron pool (30 slow and 30 fast motor units). Both models have 60% slow fibers and 40% fast fibers. (a) shows the dynamic force responses to the sinusoidal activation (U: 0-1 at 3 Hz) and constant length (L mt = 16 cm) inputs; (b) shows the dynamic force responses to sinusoidal muscle stretching (L mt : 14-18 at 3 Hz) and constant activation (U=0.5) inputs. .......................................59 Figure II.7: The comparison of computation time of different recruitment models and implementation. The original Natural Discrete model by Brown and Cheng includes L eff and is implemented by interconnected Simulink basic-blocks; the Natural Discrete and Natural Continuous models do not include L eff and are implemented using Simulink CMEX S- functions. ......................................................................................................................................60 Figure III.1: Diagram for the sensorimotor systems model integration. The Virtual Muscle receives α-command (u) and produces force outputs to drive SIMM musculoskeletal model for joint kinematics. The joint dynamics in turn passes the changes in musculotendon length (L mt ) back to the Virtual Muscle models. GTO model receives muscle force output from Virtual Muscle blocks and produce I b afferent firings. Muscle Spindle model is driven by dynamics, static γ-commands ( γ stat , γ dyn ) and muscle fascicle lengths (L ce ) from Virtual Muscle blocks, and outputs primary (I a ) and secondary (II) firings..............................................65 Figure III.2: The virtual arm musculo-skeletal model in SIMM. (a). The planar view of right arm skeletal model demonstrating the segments (ribs, sternum, clavicle, scapular, humerus, radius, ulna and bones of the wrist and hand), the joint coordinate systems, the degrees of freedom (DOF), shoulder flexion/extension (F/E), elbow F/E and forearm pronation/supination (P/S). (b). The origin, insertion points and paths of the 15 major muscle elements across the three joints (shoulder, elbow and forearm) of skeletal model. The 15 muscle elements include: deltoid anterior (DA), deltoid posterior (DP), clavicle portion of pectoralis major (PC), supraspinatus (SS), infraspinatus (IS), biceps long (Blh), biceps short (Bsh), triceps long (Tlh), triceps lateral (Tlt), triceps medial (Tmd), brachialis (BS), brachioradialis (BR), pronator teres (PT), pronator quadratus (PQ), and supinator (SP). The thorax functions as a ground and was included as the reference of the clavicle and scapular, which provided the attachment sites for shoulder muscles. The 2 DOF and 6 muscles used in later simulation are bold labeled in (a) and (b) .............................................................................67 Figure III.3: S-function implementation of the spindle model. (a) Matlab/Simulink S-function implementation of each intrafusal fibers in the spindle model. The spindle model which consists of three intrafusal fiber models (bag1, bag2, chain); two afferent firing summation nodes (I a /II afferent firing models); and the partial occlusion effect in primary afferent firing. Each intrafusal fiber model responds to two inputs: fascicle length and the relevant fusimotor drive. The spindle model generates two outputs: I a and II afferent activity. (b) Frequency-to-activation conversion and saturation effect for bag2 and chain fibers. The S- function spindle model skips the frequency-to-activation conversion in the original model by directly scaling the gamma static commands ( γ stat ) for bag2 and chain fibers according to their different saturation frequencies. Receiving the same γ stat drive (0-1), bag2 saturate at 0.5 corresponding to its saturation frequency of 100 Hz, whereas chain fiber saturate at 1, corresponding to its saturation frequency of 200 Hz. ...................................................................74 Figure III.4: Simulink implementation of a GTO model representing piece-wise linear static relation between the total muscle force and I b afferent firing to account for the ensemble response of the GTO in the muscle. The higher slope at the lower levels of muscle forces is due to the recruitment of tendon organs, which makes the predominant contribution to the total response at low forces...........................................................................................................75 x Figure III.5: Moment arm matching of shoulder flexion and extension (F/E) (a), elbow F/E (b), and forearm pronation and supination (P/S) (c). The dotted lines are the experiment measurements, the solid lines are from SIMM model. The muscle paths are validated through matching the model moment arm profiles of shoulder, elbow and forearm joints with those obtained from experimental measurements. Since the shoulder F/E moment arms of two heads of biceps and the long heads of triceps muscles are only measured at one single joint configuration, the model moment arm profiles are compared to the single data points in (a)..............................................................................................................................................78 Figure III.6: Diagram showing the operating range of fascicle length (L ce ) predicted for each muscle in the virtual arm model. The left and right edges of each dark bar indicate the values of min ce L and max ce L corresponding to the minimum and maximum physiological lengths of each musculotendon elements respectively, and the dark bars illustrate the portions of the length – tension curve on which muscle develops active force. ...................................................79 Figure III.7: Dynamic responses of one simulation at hand position C: (a) hand kinematics in the Cartesian coordinate originated at shoulder, (b) shoulder and elbow joint kinematics, (c) α command of the six muscles, (d) L ce , and (e) L ce , (f) muscle force F m , (g) spindle I a and (h) II firings, and (i) GTO I b firings of 6 muscles. .................................................................................83 Figure III.8: (a) Three hand positions (A, B, C) in space achieved with three sets of open-loop activations. (b) With the constant fusimotor drives ( γ stat = γ dyn =0.3), the steady state values of I a and II firings (left axis, red) and muscle fascicle lengths (right axis, black) are shown at the three hand positions (A, B and C) in the workspace. Each plot is for one of the six muscles. There is a monotonic relation between afferent firings with the arm positions and muscle fascicle lengths. ................................................................................................................84 Figure III.9: The relation between equilibrium I a (1 st row) and II (2 nd row) firing with γ ( γ stat = γ dyn ) activations for the six muscles at the three workspace locations (A, B and C). The figure demonstrates that once the primary and secondary afferents were activated, the firing frequency shows a monotonic relationship with the fusimotor drives at all of the three hand positions........................................................................................................................................85 Figure IV.1: (a) Two-joint, 6-muscle virtual arm model. The model is kinematically constrained n the horizontal plane at shoulder level; the 2 DOFs includes shoulder horizontal flexion/extension, and elbow flexion extension, each of which is modeled as hinge joint. q 1 , q 2 are the generalized coordinates representing shoulder and elbow flexion angles respectively. The six muscles include two pairs of mono-articular antagonists and one pair of biarticular antagonists; the muscles are color-coded that is consistent with the muscle control patterns shown as Fm bar plots in Fig IV 4.7-4.9. (b) A schematic figure showing the coordinate frame and the three task configurations.................................................................95 Figure IV.2: A schematic diagram showing the two sets of desired axis of accuracy (DAA): (1) four axes (a1-a4) are evenly distributed on the extra-personal space, separated by 45 o ; (2) two orthogonal axes are defined in the intro-personal space according to the arm manipulability ellipsoid: m1 is the major axis and m2 is the minor axis of the ellipsoid. ..........102 Figure IV.3: Force variability expressed as the SD as a function of mean force levels at full range of Deltoid Posterior (DP) muscle operating lengths (L mt =14~18cm). Given Gaussian distributed SDN with the noise amplitude proportional to the mean level of activation (CV u =0.2), the natural continuous recruitment scheme of VM can generate signal-dependent force variability whose SD is linearly scaled with force magnitude with a slope, CV=2%, within the range of experimental measurement. .........................................................................103 xi Figure IV.4: Linear spectrum analysis of the natural continuous Virtual Muscle model of deltoid posterior at operating point of u=0.8, L mt =16cm and Fm=140N. (a) Time series of broad- band white noise in muscle activation input u, and its amplitude spectrum. (b) Magnitude and phase frequency transfer function of the Natural Continuous Virtual Muscle. (c) Time series of noise in muscle force output, and its amplitude spectrum. ...........................................104 Figure IV.5: Elbow joint kinematics variability (SD) as a function of muscle co-activation level at joint angle of 80 o , and isometric joint torque at 0.5Nm. .............................................................104 Figure IV.6: Three types of dynamics responses of shoulder and elbow joint kinematics from the activation patterns that produces same equilibrium posture and static end-point force. The there types are classified by the level of their stability: (1) unstable – joints can not be stabilized during the first 15 second of deterministic simulation; (2) meta-stable – joints can be stabilized in the deterministic simulation, but once the SDN is introduced, non-stationary oscillation occurs; and (3) stable – joints can be stabilized in both deterministic and stochastic simulations. ................................................................................................................105 Figure IV.7: Multi-objective optimization analysis of task configuration 1 (T1) when the error is measured by the size of hand random distribution (Size Hxy ). (a) Log-log plot of Size Hxy versus Energy consumption of all the control points ( ▪), the four identified Pareto frontiers ( □-1,2,3,4) , and one point ( ■-1*) that is the result of proportional scaled control pattern from Pareto frontier 1. The same sets of control points are plotted against same log scale of Energy for (b) shape, (c) orientation of the Hxy, and (e) shape, (f) orientation of the hand stiffness ellipses. (d) shows the hand random distribution (Hxy) traces and the contours of the five points. The four PFs are plotted on the same scale so relative size can be compared, whereas point 1* is in a much larger scale. The task configuration is also illustrated schematically on the first Hxy distribution. (g) shows the end-point stiffness ellipses of the five control points. The stiffness along the major axis of the ellipses is labeled at the center of each ellipsis. (h) shows the muscle control patterns of points 1,2,3,4 and 1* as stacked bar plots of normalized muscle force vectors. The length of each segment represents the relative contribution of each muscle in terms of muscle forces. From bottom up, the colored segments represent the relative contribution of shoulder flexor (PC-yellow), shoulder extensor (DP-green), biarticular flexor (Bsh-black), biarticular extensor (Tlh-blue), elbow flexor (BS-red), and elbow extensor (Tlt-magenta). (i) shows the absolute values of muscle forces. (j) shows the muscle activation levels.............................................................................107 Figure IV.8: Multi-objective optimization analysis of task configuration 1 (T1) when the error is measured by the standard deviation (SD) along a1 (see column 1) and a3 (see column 2) axes. (a1) and (a2) are the muscle control patterns of the Pareto frontiers as stacked bar plots of normalized muscle force vectors. The length of each segment represents the relative contribution of each muscle in terms of muscle forces. From bottom up, the colored segments represent the relative contribution of shoulder flexor (PC-yellow), shoulder extensor (DP-green), biarticular flexor (Bsh-black), biarticular extensor (Tlh-blue), elbow flexor (BS-red), and elbow extensor (Tlt-magenta). (b1) and (b2) are the absolute values of muscle forces. (c1) and (c2) are the muscle activation levels. (d1) and (d2) are the Log-log plot of SD versus Energy consumption of all the control points ( ▪) and the four identified Pareto frontiers ( □). The same sets of control points are plotted against same log scale of Energy for (e) size, (f) shape, (g) orientation of the Hxy, and (i) shape, (j) orientation of the hand stiffness ellipses. (h1) and (h2) are the hand distribution (Hxy) traces and the contours of the Pareto frontiers. They are plotted on the same scale so relative size can be compared. The task configuration is also illustrated schematically on the first Hxy distribution. (k1) and (k2) are the end-point stiffness ellipses of the Pareto frontiers. The stiffness along the major axis of the ellipses is labeled at the center of each ellipsis. ........................................................111 xii Figure IV.9: Multi-objective optimization analysis of task configuration 1 (T1) when the error is measured by the standard deviation (SD) along m1 (see column 1) and m2 (see column 2) axes. (a1) and (a2) are the muscle control patterns of the Pareto frontiers as stacked bar plots of normalized muscle force vectors. The length of each segment represents the relative contribution of each muscle in terms of muscle forces. From bottom up, the colored segments represent the relative contribution of shoulder flexor (PC-yellow), shoulder extensor (DP-green), biarticular flexor (Bsh-black), biarticular extensor (Tlh-blue), elbow flexor (BS-red), and elbow extensor (Tlt-magenta). (b1) and (b2) are the absolute values of muscle forces. (c1) and (c2) are the muscle activation levels. (d1) and (d2) are the Log-log plot of SD versus Energy consumption of all the control points ( ▪) and the four identified Pareto frontiers ( □). The same sets of control points are plotted against same log scale of Energy for (e) size, (f) shape, (g) orientation of the Hxy, and (i) shape, (j) orientation of the hand stiffness ellipses. (h1) and (h2) are the hand distribution (Hxy) traces and the contours of the Pareto frontiers. They are plotted on the same scale so relative size can be compared. The task configuration is also illustrated schematically on the first Hxy distribution. (k1) and (k2) are the end-point stiffness ellipses of the Pareto frontiers. The stiffness along the major axis of the ellipses is labeled at the center of each ellipsis. ........................................................113 Figure IV.10: Summary of directional control of the three task configuration: T1, T2, T3. First column shows the schematic diagram of the three task configurations. Their directional controls along the two sets of DAAs are presented in the same row. For each DAA set, the Pareto frontiers identified along each of the two orthogonal axes are presented as the Hxy size versus SD curves. The marker type represents the axis along which the SD is measured: ○-a1, □-a3, +-m1 and x-m2. The control patterns that correspond to the minimal SD (denoted as red marker in Hxy-size vs SD curves) along the two orthogonal axes are shown as the stacked bar plot.................................................................................................................114 Figure IV.11: Plots of the shape (column 1) and orientation (column2) of Hxy distribution versus end-point stiffness ellipses for the three task configurations. .....................................................115 xiii ABSTRACT Precise control of posture is necessary for the success of any motor actions. The strategies whereby precise movement is achieved reflect many complexly interacting factors, including the tendency of active muscles to generate substantial motor noise as well as to influence the mechanical impedance of a limb with multiple degrees of freedom. A realistic model of the human sensorimotor system would be instrumental in understanding the underlying principles whereby the CNS controls the kinematic effects of signal dependent noise (SDN) in the actuators of a redundant musculoskeletal system. In this work we first constructed a comprehensive sensorimotor system model of the human arm that accurately represents the physiological properties of human muscles and proprioceptors. A simplified two-joint six-muscle arm model was then used for the simulation study of hand postural control in the face of SDN. In single-joint systems, co-contraction has been shown to be effective to reduce the kinematic effect of SDN. But the strategy is energetically inefficient, and extension to a multi-joint system is not obvious because of the intersegmental dynamics and redundant musculature. In this study, we used Monte Carlo simulation to explore how various patterns of muscle activation affect both the magnitude and direction of such perturbations and the energy consumption under various limb positions and loading conditions. Multi- objective optimization analyses identified a set of control patterns that satisfied both precision and energy criteria. We found that lower postural variability could generally be achieved with higher energy consumption, but the relative levels of activation of each muscle changed drastically at each control level rather than simply scaling proportionately. The end-point variability could be measured and minimized for a task-specific axis, which also affected which strategy provided the best performance, which depended also on limb configuration and external loads. Interestingly, no correlation between end-point stiffness and kinematic variability was found. We conclude that in the presence of SDN, the CNS can utilize muscle redundancy to maximize task performance but that the set of optimal strategies does not follow simple rules such as stiffness modulation and co-contraction, and probably needs be learned rather than computed by the CNS. 1 CHAPTER I: INTRODUCTION The research presented in this thesis is based on computational biology, an interdisciplinary approach to integrating biological motor systems and modern control theories. It uses a realistic sensorimotor system-level model of the human arm to address two main questions, namely: 1. What are the potential kinematic effects of signal-dependent noise in the actuators of a redundant musculo-skeletal system? 2. What strategies might the central nervous system (CNS) use to control these effects? In particular, can strategies based on optimal control be used to predict emergent behavior? The question of how humans control their motor behaviors, such as reaching out the hand and holding a cup of tea, may seem a simple one: we carry the hand from the start position to a tea cup and pick up the cup; the rest is Newton’s third and second laws. A little reflection reveals, however, that skilled motor behaviors, from the graceful leap of a ballerina to a precise use of screwdriver, appears effortless but reflects an intimate interaction among neural, sensory, and motor systems, their muscle-body dynamics, and their environments. We are not aware of the inherent difficulties of motor control until we observe the clumsy movements of a child, the motor challenges faced by patients with neurological disorders, or the yet unsolved control problems in robotic systems. Among these difficulties, how the CNS coordinates the redundant and noisy neuro-musculo-skeletal apparatus to achieve desired task performances is a central question that has puzzled neuroscientists for nearly a century. This has led to three broad approaches to motor coordination: 1. Kinematic analyses in human motor behavior (Higuchi et al. 2002; Vereijken et al. 1992), mammalian locomotion (Lacquaniti et al. 1999) and robotic skill acquisitions (Berthouze and Lungarella 2004) have revealed a “freezing” strategy in resolving the degrees-of-freedom (DOFs) problem. 2. Another related approach to the redundancy problem at the levels of forces, muscles and motor units is “synergy”: multiple muscles are bound together such that a central control signal jointly and proportionally activates all muscles in the synergy. 2 3. Finally, CNS can select a performance criteria based on certain mechanical, engineering, psychological costs, where a minimum or maximum optimization principle is applied to find a unique solution. Although these approaches were often considered as logical alternatives, they do not contradict each other. Indeed, they may need to be combined complementarily, for instance, learning a new “synergy” through optimizing a performance criterion, or, selected according to specific circumstances, for instance, searching for the optimal control strategy that minimizes the effect of inherent motor noise. Many motor control studies dedicated into understanding control strategies of motor coordination have amassed vast amounts of data, but due to various reasons, no single theories have encompassed the whole problem in resolving motor redundancy and motor noise. Mathematical models, at various levels and complexities, can play a critical role in synthesizing parts of these data by developing unified neuro-musculo-mechanical descriptions of complex multi-joint multi- muscle coordinated motor behaviors. And in this exercise, they can help to identify the convincing testing methods, and guide the understanding of biological motor systems as well as bio-inspired robots. This dissertation work introduced a general problem of motor coordination of redundant, variable motor systems. We took the specific case of multi-joint muscle-redundant arm posture control tasks, and synthesized the redundancy problem with inherent property of motor noise through dynamical simulations using a realistic musculoskeletal arm model. We examined how the internal motor fluctuations (SDN) propagate through the complex neural, muscle and joint dynamics. The findings helped us to outline the plausible control principles that the CNS might adopt when controlling redundant and noisy motor plant. This work provided revolutionary implications in noise control in human motor systems, and would encourage further efforts either in developing the computational control theories or experimental examination. i. ORGANIZATION OF THE THESIS This thesis is presented in the order of the research process: first is definition of the problem with broad literature review – demonstrate why and what part of motor control problem we are addressing, and what method and tools we need to address them (chapter I); then is the development of necessary tools – 3 present how the model of sensorimotor system was developed (chapter II, III); finally the use of the model to address the formed research questions – explain, demonstrate and answer the focused questions. Furthermore, along the main text of the thesis, inserted boxes are used to present the representative biological facts from the literature, comparison with current modeling results, and important engineering and mathematical concepts. The remainder of Chapter I includes the review of the literature. The widely known background materials and the evolvement of rationale for this thesis work are elaborated here. They encompass the areas in properties of biological and artificial motor system (sections I.ii redundancy and I.iii noise), theoretical framework of sensorimotor control (sections I.iv optimality and I.v impedance control) and mathematical models of neuromusculoskeletal systems (section I.vi). Chapter II presents Virtual Muscle model – a sophisticated mathematical description of biological skeletal muscles based on a series of muscle characterization experiments. The model captures many important muscle physiology properties that have been absent in previous muscle modeling efforts, and provide a generalized tool (in a form of freely released software VM4.0) in studies of sensorimotor system control. This part of work also is building the foundation of the dissertation – computation investigation of motor redundancy and motor noise. Chapter III describes the building process and validation of a sensorimotor system model of human right arm called Virtual Arm model. The arm is kinematically constrained in horizontal plane at shoulder level but includes three degrees of freedom (shoulder, elbow flexion and forearm pronation); Virtual Muscle is used for modeling 15 major muscles according to realistic musculo-skeletal geometry (including muscle moment arm and series- and parallel-elasticity). Each muscle is also equipped with proprioceptive sensors including a realistic model of muscle spindle model and Golgi tendon organ. The model is developed for a general purpose of sensorimotor study. A subset model of Virtual Arm (2-DOF, 6-muscle) without muscle sensors is also used in this research. Chapter IV contains the rationale, methods, results and discussions on the research of motor redundancy and motor noise. A realistic musculo-skeletal model of the human arm developed in Chapter III 4 was used for the dynamic simulation of multijoint posture control tasks with redundant, noisy neuromuscular structure. Motor commands were corrupted with signal-dependent noise, and the resulting hand/endpoint kinematic variability was characterized. Various system states – muscle force, coordination patterns, limb stiffness, and energy costs – were quantified and examined to reveal the underlying control strategy. Finally, summary and future works are contained in Chapter V. ii. MOTOR REDUNDANCY a) What is motor redundancy? Voluntary actions of humans and other animals consist of coordinated movements directed toward specific motor goals. To drink a cup of tea on a table, person would first locomote, transporting the body through space to the table while remaining upright; then reach out the hand and grasp the cup; finally transport the cup stably to mouth without spilling the tea. At every level of analysis of the motor system, there are more elements contributing to performance than are absolutely necessary to perform these motor tasks, at least when they are defined according to these simplistic goals. This problem of redundancy has been recognized as a central one from the earliest days of the scientific study of motor control. About half a century ago, A.N. Bernstein (Bernstein 1967, 1996) identified the “degrees-of-freedom (DOFs) problem” in neuromuscular control, which may be exemplified as follows. “Typical limb movements, such as reaching to pick up a small object from a table, require precise fingertip placement, but leave intermediate hand, arm, wrist, elbow and shoulder joint angles and positions undetermined” (Holmes et al. 2006). This kinematic redundancy addresses the redundant modes of movement in the multi-link structure of the human body. How the central nervous system plans and coordinates these DOFs while producing the desired performance in a motor task? Are these DOFs independently programmed or they are coordinated following invariant patterns? When professional blacksmiths stroke the chisel with the hammer, Bernstein noticed that variability of the hammer tip trajectories across a series of strikes was smaller than variability of the individual joint trajectories of the subject’s arm. This observation suggested to him that the CNS does not try to find a unique solution for the problem of kinematic redundancy by 5 eliminating redundant DOFs but rather uses the apparently redundant set of joints to ensure more accurate performance of the task, and to compensate for environmental disturbances in the course of the movements. Analogous problems arise at the level of muscles. Typically, limbs have fewer DOFs than the number of muscles actuating them. For example, in human upper extremity, 15 DOFs representing the shoulder, elbow, forearm, wrist, thumb, and index finger are actuated by 50 muscle compartments (Holzbaur et al. 2005); and 4-DOF, 3-link index finger are controlled by seven muscles (Valero-Cuevas et al. 1998). There is the potential for muscles acting on the joint to join in varying combinations of activation and still produce the same directional force output or movement trajectories. Which combination of muscle forces does the CNS select to achieve a particular amount of torque or endpoint force? Are the particular patterns of activation associated with the environmental demands in performing a motor task? Changing muscle activation patterns in a multi-joint arm movement were shown to adapt the limb stiffness to the instability in a spring force field (Franklin et al. 2007; Franklin et al. 2004), and suggest that the CNS is exploiting the seemingly redundant muscle sets to compensate the environmental instability to successfully perform the tasks (see Section I.v). b) How to resolve motor redundancy? A large body of literature has addressed the various resolutions to the redundancy problem; they can be categorized as “elimination”, “synergy” and “optimization” approaches. All of these imply a concept – collapse of dimension: the emergence of a low-dimensional attractive submanifold in a much larger state space. The coordination models of central concern in this dissertation, to be introduced later in this section, at lease suffice to explain the observed mechanical patterns associated with this concept. Elimination A principle of elimination indicates that the CNS solves the redundancy problem by reducing the number of DOFs to the ones necessary to perform the task. This is a “freezing” approach that had been postulated by Bernstein, and has been commonly invoked in contemporary studies of human motor behavior (Higuchi et al. 2002; Vereijken et al. 1992), mammalian locomotion (Lacquaniti et al. 1999) and robotic skill acquisition tasks (Berthouze and Lungarella 2004). 6 This approach appears to be associated with a posture principle observed in some circumstances, where the motion is restricted to a low-dimensional subspace within a high-dimensional joint space. Evidence for such constraints comes from Donder’s law which, when applied to eye movements, states that the angular gaze positions do not routinely make use of all three DOFs but are constrained to a two- dimensional surface when the head is stationary (Donders 1847). The work by (Gielen et al. 1997) explores the validity of Donder’s law for arm movements. When subjects point with their hand to targets, the end- configuration of their arm is reproducible, irrespective of the initial posture of the arm. This reproducibility can be accounted for by identifying particular rotation axes defined by the initial and the final pointing directions, such that rotation about these axes completely characterizes the movement, while rotations about remaining axes reflecting the redundant DOFs are very small. This implies a reduction in the number of DOFs (collapse of dimensions) because number of variables that need to be assigned new values during the movement is smaller than the number of mechanical DOFs. Similar constraints were also verified in arm throwing movements (Hore et al. 1994). The nature of these constraints may provide information about the properties of the underlying coordinate system used by the CNS. But the notion of elimination of DOFs has been invoked primarily at the kinematic level, not account for how the redundant muscle patterns are selected among almost infinite solution spaces. Moreover, freezing the joint at a subset of dimensions places a nontrivial problem to the CNS: how the muscles work together to constrain the joint at its frozen value, while producing the desired kinematics in non-constrained planes? For instance, (Latash et al. 1999) showed complex control signals to a “frozen” wrist joint when a pure elbow movement was generated. Synergy The second class of solutions is synergy. In fact, “synergistic phenomena have been prodigious source of evolutionary novelty” (Corning 1995); it has been proposed that the functional advantages associated with various forms of synergy have been an important cause of the “progressive” evolution of complex systems over time. “Though it often travels in disguise, synergy can be found in the subject-matter of most, if not all of the academic disciplines” (Corning 1995). The notion of motor synergies, or high-level “control knobs” that have distributed action over sets of low-level actuators, arose in the context of motor 7 coordination (Bernstein 1967), and can be traced back to early classical works by (Hughlings Jackson 1889) and (Sherrington 1910). While synergies mean different things to different people (Latash 1999), dimensionality reduction is the signature of synergistic control. It can be most clearly illustrated in the context of muscle synergies: “Multiple muscles are bound together such that a central control signal jointly and proportionally activates the muscles in the synergy”(Latash et al. 2007). This coupling (synergy) of the independent elements in the motor effectors can render an effective control space of the CNS with a much smaller dimensionality than the actual mechanical system. The idea is to simplify the control problem faced by CNS on redundant biomechanical system by reducing the number of the variables it has to manipulate. However, “such a view implies that the number of synergies is small, and their structure is fixed in advance of learning a given task, and is therefore task-independent” (Todorov and Ghahramani 2004). This is problematic in that the control is only easy if the synergies are known beforehand, and if a novel task demands a new set of synergies, this learning process will present to the CNS a problem with same dimensionality as the original control problem. Optimization Optimality principles in motor control offer a very different perspective in solving the redundancy problem. These approaches essentially are seeking the unique solutions by minimizing or maximizing cost functions that obey certain mechanical, engineering, psychological, or composite criteria [reviews see (Seif-Naraghi and Winters 1990)]. The optimality yields computational-level theories, which try to explain why the system behaves as it does, and to specify the control laws that evolve the synergistic actions observed in these stereotypical behaviors. For example, in arm movements, it has been suggested that trajectories and muscle activation patterns are selected to optimize a cost that is integrated over the movements, such as jerk (time derivative of acceleration) (Hogan 1984), or torque change (Kawato et al. 1990; Nakano et al. 1999; Uno et al. 1989). However, there has been no principled explanation why the central nervous system should have evolved to optimize such quantities, other than that these models predict smooth trajectories. Indeed, since these qualitative behavior manners are directly transformed into the cost functions such that the minimization of them will self-evidently reproduce the observed quality. This circular reasoning process is highly sensitive to the physical plant controlled, the task studied, the 8 performance index optimized and the control constraints imposed, but ignores the fundamental physical reasons inherent in the generic biological motor effectors. One such reason, as we shall argue in the following section, is the instantaneous signal fluctuations and motor noise in the sensorimotor system. More sophisticated optimal control models concerning the noisy motor plant include the TOPS model by (Hamilton and Wolpert 2002; Harris and Wolpert 1998), and optimal feedback control by (Todorov 2004; Todorov and Jordan 2002). These models have predicted the invariant features of human movements or “motor equivalence” (see Section I.iii) that evolves from a combined feature of motor redundancy and motor noise. iii. MOTOR VARIABILITY AND MOTOR NOISE Bernstein’s subjects were all professional blacksmiths: they had performed the same movement hundreds of times a day for years. For these highly skilled tasks, we can easily assume that the CNS can just apply certain rules (for example optimization) to find a best solution and program its control laws adhering to this unique solution to ensure 100% accuracy hitting the target. But there was variability both in the tip trajectory of the hammer and in the joint configurations of the subject’s arm. As humans, we can never conduct any single motor task in the exact same ways either for the dynamic process or the task goals. When throwing the basketball we may hit the basket in the first two shots but miss it in others. We take for granted that our behavior is variable, and that repeated attempts will have variable results, but what is the source of this variability? When we err, the reasons can be something went wrong during the movements, either a perturbation from the environment or fluctuations (noise) of our motor commands; it is also possible that variability arises during motor preparation, well before the first muscle contracts. Identifying the origin of the variability is critical to the study of motor control. Not only is variability a part of the behavior that must be explained, but hypotheses about the motor-control strategies are fundamentally linked to hypotheses regarding the noise those strategies combat. Here, a linguistic problem has to be cleared, as each field has developed its own interpretation of terms such as variability, fluctuation and noise. The term variability is generally used to refer to changes in some measurable quantity, such as finger force, hand positions or movement duration. They are mostly 9 associated with trial-to-trial variation. Importantly, the term variability does not indicate that a particular mechanism has generated the variability, and does not suggest whether the variability is beneficial or detrimental (Faisal et al. 2008). Whereas noise, existing in various levels or compartments of the motor system, is often used as opposition to the word deterministic and acts as one of the major sources of the observed variability. It is, as defined in the Oxford English Dictionary, ‘random or irregular fluctuations or disturbances which are not part of a signal, or which interfere with or obscure a signal or more generally any distortions or additions which interfere with the transfer of information.’ In the context of this dissertation, one conceptual grasp of the difference between the two terms is that, noise is more towards an inherent property of the motor signals resulting from various random processes at biochemical, electrical or molecular, cellular levels [reviews see (Faisal et al. 2008)], whereas variability is more “controllable” , for example by some neural mechanisms (Latash et al. 2005), so that important performance variables are stabilized (see I.iii.c). Thus we will refer to the inherent fluctuations of muscle force output as motor noise, and changes in various task parameters such as hand force and kinematics as variability. In this section, we will first review various opinions existing in the current literature concerning the sources of motor variability – central or peripheral. Among the peripheral sources, we will then focus our attention to motor noise (random fluctuations of muscle force), as it is one of the major concern in this work. Here the physiological origins of muscle motor noise and its signal-dependent features will be presented. In the last part of this section, the implication of noise on motor control theories will be reviewed. a) Where does the movement variability come from? The origin of the movement variability, if categorized by its anatomical structures, could be: (1) the central motor preparation which may or may not be related to the sensorimotor transformations, and/or (2) the peripheral process of the movement execution that takes place after the movement commands have been planned. In this context, the movement variability is often quantified from trial-by-trial measures (Churchland et al. 2006a; Gordon et al. 1994a), which can arise from: (a) noises that exist across the entire loop of signal-processing in sensorimotor system, and (b) the deterministic properties of the system, for example, the initial state of the neural circuitry will vary at the start of each trial, leading to different 10 neuronal and behavioral responses. The variability in the response will be exacerbated if the system’s dynamics are highly sensitive to the initial conditions (Faisal et al. 2008). Cautions should be taken when considering various opinions about the locus of the origin of behavioral variability. The discussions about the central or motor-preparatory origin very often do not clearly differentiate the noise and deterministic process, due to integrative nature of signal processing (e.g. sensorimotor integration and state estimation), whereas the peripheral sources are mostly associated with the inevitable noises (motor noise or the noise of the muscle forces) propagated through the descending or forward segment of the sensorimotor loops. The nature of the motor noise (the signal-dependence and Gaussian distribution) as we will review in next section, has been determined by the physiological properties of the signal relays (assembles of motor neurons and muscle fibers), thus is in a sense uncontrollable, at least at the level of force outputs. Central preparatory origins The idea of central origins of the motor variability is firstly derived from behavioral observation. Gordon et al (Gordon et al. 1994b) pointed that motor variability can be parsed into different soucese, and inferred the central source through kinematic analysis of planar reaching movements. They observed that when healthy subjects perform two-joint, point-to-point reaching movements along different directions without visual feedback, the variability patterns of movement direction and extent were consistent across different movement speeds and the movements involving different joints. These results ruled out the possibility that the variability may come from biomechanical factors, such as multi-segment limb inertia (Hogan 1985), and suggested that they primarily originate from errors of movement plan in the central processing units. But absence of a direct link to the neural activities can not prove that the exclusive source resides in motor planning. A recent study by (Churchland et al. 2006a) directly measured single cortical neuron activities (dorsal premotor and primary motor cortex) while a monkey was performing highly practiced reaching movements; the strong correlation between the fluctuations of neural activity and the peak velocity variability identified during a preparatory period (100’s ms before movement onset) suggested certain motor planning or preparation process is an important source (50% or more) of the subsequent movement variability. The idea that planning noise is the major source of movement variability has been very influential (McIntyre et al. 1997; van den Dobbelsteen et al. 2001; Vindras and Viviani 1998). 11 However, these studies failed to provide evidence that all the observed variability arises in the planning process and not at other stages; neither could they address the proportions of the central and peripheral contributions. Peripheral execution origins Putting aside the debate of which factor plays the dominant role in changing the motor behavior, the movement execution inevitably gives rise to motor variability because motor output (motor commands, muscle contractions, or muscle torques) is inherently variable. This is evident in motor noises (random fluctuations of muscle force) produced during isometric muscle contraction (Jones et al. 2002; Schmidt et al. 1979). Several authors (Harris and Wolpert 1998; Meyer et al. 1988; Schmidt et al. 1979; Vangalen and Dejong 1995) have proposed that noise arising in the motor output stage gives rise to substantial movement variability, which could explain the speed-accuracy trade-off known as Fitts’ law (Fitts 1954). A series of studies has correspondingly suggested that a significant amount of movement variability is, indeed, coming from such a peripheral source. For instance, when human makes repeated center-out reaching movements at various instructed speeds, the variability of the speeds and movement extents can not be explained by planning noise but is well explained by noise in movement execution (van Beers et al. 2004). Computational modeling of various types of noise (signal-dependent and –independent) inferred the possible contribution of each type; and the results suggested that execution noise accounts for at least a large proportion of movement variability. A later work by the author (van Beers 2007) studying saccadic eye movements with modeling inferences extends the findings by quantifying their relative contributions among sensory (~57%), signal-dependent (21%) and constant noises (~21%) . The immense interests in identifying the sources of movement variability obviously come from the importance of them in inferring the control strategies that could suppress the noises and limit their harm. And surprisingly different conclusions reached by the reviewed studies imply considerable challenges to answer the question. Firstly, we can not rule out the possibility that the major source could differs for different body parts; and even if the origin is universal, the precise identification would require direct measuring of signals across entire loops of sensorimotor chains, or advanced signal processing technique to separate the composite effects. These challenges are universal in experimental motor control studies. By 12 contrast, one attractive feature of modeling is its convenience in segmenting individual components of the sensorimotor chain and noise types to examine their behavioral consequences [e.g. (Jones et al. 2002; van Beers 2007; van Beers et al. 2004)]. The results presented in this thesis (Chapter IV) derive from such a modeling approach, and focus our attention on signal-dependent execution noise and its effects on motor performance. To do so, we will first need to review the nature and sources of the signal dependent motor noise in the following section. b) What determines the nature of signal dependent motor noise? Recent studies (Hamilton et al. 2004; Jones et al. 2002; Osu et al. 2004; Sosnoff et al. 2006) have attempted to characterize online neuromuscular noise precisely because it is proposed to be the key factor limiting performance. Unlike the general neural noise (e.g. cortical) that has spiking statistics that approximate a Poisson process (Churchland et al. 2006b; Faisal et al. 2008; Tolhurst et al. 1983), the noise in motor neuron and muscle force can be described by a distribution close to Gaussian (Gomez et al. 1986). In addition, direct recordings of motor neuron signals in monkeys during saccades (Hu et al. 2007), and in cat (Gomez et al. 1986) and goldfish (Pastor et al. 1991) during fixation have shown that their firing exhibits “signal-dependent noise” (SDN). This is a white noise with zero mean and a standard deviation (SD) proportional to the absolute value of the signal, such that the noise level can be quantified by a constant coefficient of variation (CV=SD/mean). Such a relation has been observed for both the firing of motor neurons (Pastor et al. 1991) and for force production during isometric contraction (Galganski et al. 1993; Jones et al. 2002; Laidlaw et al. 2000; Schmidt et al. 1979; Slifkin and Newell 1999). In standard model, the CV of the ISI variability in the motoneurons are reported to be 20% (Fuglevand et al. 1993; Jones et al. 2002; Selen et al. 2005) and CV of isometric force is 2-5% over the large range of force magnitude (Jones et al. 2002; Taylor et al. 2003). Furthermore, while the power spectrum of neural discharge rates spans a large range, both experimental and simulated forces have demonstrated significant power in the frequency spectrum below 2 Hz (Taylor et al. 2003). This low frequency spectrum is likely a result of the mechanical characteristics of force production (low-pass filtering). 13 Where is this signal-dependent noise in force output generated? Is it in the planning or execution of motor commands, or is it coming from the peripheral machinery of the skeletal muscles? Examining sources of the signal-dependent motor noise is crucial for constructing neuromuscular system models with appropriate variables that are akin to the underlying control strategies on movement variability. Obviously, part, if not all of the variability in continuous isometric force production arises from the statistical variability and synchrony in the discharge of motor neurons supplying the muscle (Semmler et al. 2001; Semmler and Nordstrom 1998). Any process contributing SDN at any stage in the evolution of a motor command will affect the outflow of action potentials to the muscle. However, besides SDN, there are many different types of noise that exist in the motor command (motor neuron signals). For example, the variability of the inter-spike intervals (ISI) could be independent of the signal levels; this constant noise BOX.1| Sources of signal dependent noise of muscle forces (Jones et al, 2001) A systematic study was conducted by Jones et al (2001) to investigate the origin of the signal-dependent noise present in the continuous muscle isometric force production. The authors firstly distinguished the central (motor neuronal) and peripheral (muscle contractile machinery) sources through experimentally comparing the voluntary isometric force with those elicited by neuromuscular electrical stimulation (NMES) of the extensor pollicis longus muscle. Then a model of motor unit pool was constructed to determine factors of motor-unit physiology – including size of the pool, twitch force range, and recruitment range – that may contribute to SDN. The model was also used to investigate the effect of recruitment order as well as different types of ISI variability (constant, signal-dependent, or CV increases with mean of ISI) on the force variability. The main findings are summarized as follows. As shown in following graphs, a | During voluntary contraction, SDN was evident with the standard deviation (STD) of the force variability linearly scales with mean force level (O). However, the noise during NMES remained almost constant over the same range of mean forces ( ∆). When voluntary contraction is combined with NMES, SDN could be restored in force outputs ( □). b | The graph shows the number of motor-units in the pool was not a significant factor in determining the scaling of SDN. That is the slope of the log-log relation was not affected by the number of motor units and remained relatively constant and close to 1 over a large range. c | The range of the twitch forces in the motor-unit pool was a critical factor in determining the scaling of SDN. That is, the largest motor unit has to be at least 100 times the amplitude of the smallest. The slope dropped sharply once this range was below about 50. d | Similarly, the recruitment range, which determines the percentage of the MVC at which the last motor neuron is recruited, was also an important determinant of SDN scaling. As this range was narrowed below about 10 (i.e. the last motor neuron activated at 30% MVC), the slope of the log-log relationship began to fall sharply. In other simulation results that were not shown here, the authors demonstrated that the orderly recruitment in voluntary force production is another necessary condition in addition to the twitch force range and recruitment range in determining the linear scale of SDN. Further more, the scaling of SDN was independent of the presynaptic noise in the motor commands. That is, the linear scaling property of SDN is preserved no matter if the std of ISI is constant, signal-dependent or its CV increases with the mean magnitudes. 0 50 100 150 0.5 0.6 0.7 0.8 0.9 1 Motor unit pool size Slope b 0 200 400 600 0.5 0.6 0.7 0.8 0.9 1 Twitch force range 0 20 40 60 0.5 0.6 0.7 0.8 0.9 1 Recruitment range 0 50 100 0 1 2 3 Mean Force (%MVC) STD Force (%MVC) (a) (b) (c) (d) 14 could result from background activity of the motor neurons or many complex factors (preparatory or executional, synaptic or electrical) (Faisal et al. 2008). Or the CV of ISIs distribution can increase with the mean of ISIs (Clamann 1969; Fuglevand et al. 1993) within the subprimary range of motoneuron firing (Kudina 1999; Person and Kudina 1972; Piotrkiewicz 1999), where the ISI histogram is better fit to a Rayleigh as opposed to a Gaussian distribution (Jones et al. 2002). In such cases, there must be some other mechanisms at other levels, maybe the recruitment mechanism in the motoneuron pool (spinal cord) and/or the conversion of neural spikes into forces at each motor units, that contribute to the signal dependent feature of force variability. A systematic study of these questions was conducted by Jones et al (Jones et al. 2002). They used neuromuscular electrical stimulation to distinguish the muscle mechanical machinery from neural mechanisms as the possible source of the force variability; in addition to the experiments, a model of motor neuron pool was used to infer the effects of individual aspect of motor unit organization. The results (see Box 1) led to three important conclusions 1) the contractile machinery of skeletal muscle does not account for the signal dependent feature of the isometric muscle force variability; 2) the nature of the voluntary force production (orderly recruitment) and the physiological organization of the motor-unit pool, such as range of twitch amplitudes and range of recruitment thresholds, biases force output to exhibit linearly scaled SDN; 3) linear scaling feature of the SDN of motor output does not depend on SDN in the motor commands (ISI), but the variability in the motor drive does affect the magnitude of motor noise reflected as the values of CV. The relatively consistent measurement on the level of force noise [CV=2-5% (Jones et al. 2002; Taylor et al. 2003)] might results from relatively consistent distribution of the sizes across multiple motor units and the dynamical filtering properties of each unit. Several other studies (Selen et al. 2005; Taylor et al. 2003; Vangalen and Dejong 1995) also confirmed the importance of these properties of motor unit pool organization. Taylor et al. (Taylor et al. 2003) continued this line of research by introducing a model of motor neuron synchronization, and examined the force variability in both time and frequency domain. The major findings in this study suggest that the effect of synchronization is independent of discharge-rate variability; in fact the presence of motor-unit synchronization may act in some instances to dampen the effect of discharge-rate variability on force fluctuations. The pattern of force fluctuations does 15 appear to rely on the composite activity of the motor-unit population, but a single mechanism can not account for the entire operating range of the muscle. There is still uncertainty in the picture regarding where and what mechanisms determine the various types and features of the noises in the motor system, but there is certainty that noise, particularly signal- dependent motor noise, places substantial limitations in the information processing, and constraints CNS’s selection of control strategies. c) What does the noise and redundancy imply about movement control? We now combine the notions of motor noise and motor redundancy to account for the various emergent properties of biological motor functions as well as their impacts on the proposed control theories. As we shall recall from Section I.ii, redundant design of human motor system places a fundamental question to the neurophysiologist and control engineers: How does CNS select a particular solution from infinite number of possibilities? Whatever solution it selects, both the inputs and the outputs will In redundant motor systems where varying combinations of elementary variables (or DOF) can result in a common goal of a motor task, motor variability naturally arise either from the inherent noises within the sensorimotor system, or planning phase of a motor task. This variability is distributed among redundant DOFs with a certain structure, where those affect the task performances are smaller than those do not. This is so called “structured variability” in the motor control literature and the phenomenon has been mathematically described using a “uncontrolled manifold (UCM)” method (Scholz et al. 2003; Scholz and Schoner 1999; Scholz et al. 2000). F 1 F 2 Task goal: F TOT = F 1 +F 2 = 20 N V good =V bad V good V good V bad V bad UCM Task-irrelevant dir V good >V bad CM Task-relevant dir This concept can be exemplified using a simple two-finger force regulation task (Li et al. 1998; Scholz et al. 2003). The goal of the task is to produce a constant stable force at 20 N with two fingers pressing simultaneously on a force gauge. Variability of individual finger forces will cause the total force to deviate from the desired force level (20N), especially when the two finger forces co-vary positively with each other (positive co-variation). This is a direction that is most task-relevant, and can be treated as a “controlled manifold (CM)”; the variability along this direction will be bad variability (V bad ) because it severely hurts the task perfromance. In the orthogonal direction, the two finger forces co-vary negatively with each other, thus has the errors introduced by one finger is compensated by another finger, hence stabilizing the total force output. This direction of covariation is task- irrelevant or “uncontrolled manifold (UCM)”; the variability along this direction can be seen as good (V good ). To stabilize the task variable (F ToT ), the nervous system has been shown to be able to optimally coordinate the elementary variables (F 1 and F 2 ) in such a way that the variability distribution along the CM is kept minimal, whereas that along the UCM is uncorrected. This control structure has been proved to be energy efficient in that the CNS does not allocate control energy unnecessarily into constraining those variances that do not interfere with performance. This has been formalized as a “minimal intervention principle” by (Todorov and Jordan 2002). BOX.2| Structured variability and uncontrolled manifold 16 BOX.3| Fitts’ law Fitts’ law is a model of human movement which predicts the time required to rapidly move to a target area, as a function of the distance to the target and the size of the target. One common form of its mathematical formulations is the Shannon formulation ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + = 1 W D blog a T 2 where, T is the average time taken to complete the movement; a represents the start/stop time of the hand and b stands for the inherent speed; D is the distance from the starting point to the center of the target; and W is the width of the target measured along the axis of motion. W can also be thought of as the allowed error tolerance in the final position. From the equation, we see a speed-accuracy tradeoff associated with pointing, whereby targets that are smaller and/or further away require more time to acquire. necessarily suffer some degradation as a result of noise, so presumably the CNS should select strategies, at least in part, on the basis of their sensitivity to such noise. Noise exists in the entire chain of the sensorimotor loop and is inescapable (as reviewed in previous section, I.iii), so the motor system has to be aware of its property, origin and mechanism in order to construct effective strategies to minimize its effects on those criteria by which performance will be judged. Of course, we shall not rule out the beneficiary aspect of noise, such as information processing (Faisal et al. 2008). But it is the detrimental effects that place the principal constraints in the motor functions. Here we will firstly review some invariant characteristics of human movements despite a large pool of redundant solutions and noisy implementations. The relevance of noise property, sources and redundant biomechanical plant to various existing motor control models will then be described. Optimization, as a most popular theory in movement control, will be briefly revisited here in the context of motor noise (detailed review see next section, I.iv). Alternatives such as impedance control will be introduced as a preparation for further discussion in the following section (Section I.v). Invariant features of human movement In humans, stereotypical behaviors arise in various levels of motor functions. Voluntary arm movements appear invariant, with a characteristic near-straight hand path and a bell-shaped velocity profile in Cartesian space (Abend et al. 1982; Atkeson and Hollerbach 1985; Brown et al. 2003; Flash and Hogan 1985; Hogan 1988). These kinematic characteristics remain unchanged even when the motor effectors were subjected to external force disturbances (Burdet et al. 2006; Morasso 1981; Shadmehr and Mussaivaldi 1994; Thoroughman and Shadmehr 1999). To produce the desired kinematics, extensive experimental studies have been carried out to identify the biological control signals to the redundant muscle actuators. An invariant triphasic muscle activation pattern has been observed for both single-joint and multi-joint 17 movements (Flanders et al. 1994; Ghez and Martin 1982; Hannaford and Stark 1985; Karst and Hasan 1991a, 1991b; Lestienne 1979; Marsden et al. 1983). Though the invariance seems prevalent in many functional variables in the motor system, there is no unequivocal relationship between these “consistent” goals and motor patterns involving infinite combinations of redundant elements. This is a property first noted to be central to the functioning of motor systems by (Lashley 1933), which he called motor equivalence. For instance, the nervous system can preserve a common kinematic pattern for executing a movement (e.g. straight trajectory of hand) while varying the moving effector [e.g. the angular patterns of motion, or underlying muscular activations (Bernstein 1967; Burdet et al. 2001; Gribble et al. 2003; Levin et al. 2003)]. However, “solving the DOF problem should not only help us understand the origin of these flexibility, but also explain how the flexibility coexists with constraints inherent to the functioning of the motor system” (Guigon et al. 2007). Here it comes with the notion of motor variability and motor noise. For instance, when humans are free to move, they automatically scale movement duration with movement amplitude and choose a trade-off between speed and accuracy, known as Fitts’ law (Fitts 1954). For drawing movements, it has been found that hand velocity (V) decreases as the radius of curvature of the path (R) decreases – this has been formulated as the two-thirds power law: V=KR (1-β) , where K is a constant and β ≈ 2/3 (Lacquaniti et al. 1983; Viviani and Schneider 1991). Though it is still unclear which components of the noise train in the motor system contribute to this typical behavior (as we reviewed early in this section), the subsequent motor variances raised an important question: are there coherent patterns of motor coordination exist across individuals? Such patterns certainly exist: Bernstein, when observing the skilled blacksmiths hitting the chisel with the hammer, noticed that the trial to trial variability of individual joint trajectories of the subject’s arm are on average larger than fluctuations on the trajectories of the hammer, i.e. fluctuations in the redundant motor effectors are not equally constrained; rather, they are coordinated in a way such that only the variability that hurt the performance is limited. This so-called structured variability has been observed in wide range of motor actions such as posture (Balasubramaniam et al. 2000), locomotion (Winter 1989), 18 writing (Wright 1990), shooting (Scholz et al. 2000), pointing (Tseng et al. 2002), reaching (Haggard et al. 1995), grasping (Cole and Abbs 1986), sit-to-stand (Scholz and Schoner 1999), bimanual tasks (Domkin et al. 2002), multi-finger force production (Li et al. 1998; Scholz et al. 2003), as well as multi-muscle body sway (Danna-dos-Santos et al. 2007) and finger force production (Valero-Cuevas et al. submitted). Further more, perturbations, either from unexpected environmental changes (Cole and Abbs 1987; Robertson and Miall 1997)or inherent noise (Li et al. 1998; Valero-Cuevas et al. submitted) within the motor system are always compensated in a way that task performance rather than a specific pattern of elementary variables are maintained. The subspace where the motor variability is constrained to is often quantified as the uncontrolled manifold (Scholz et al. 2003; Scholz and Schoner 1999; Scholz et al. 2000), a mathematical concept that is consistent with the minimum intervention principle whereby the deviations from the average behavior are corrected only when they interfere with the task goals (Todorov 2004; Todorov and Jordan 2002). Noise and redundancy on control theories Breakthroughs into the understanding of human’s motor functions usually come from developed principles that successfully explain their motor regularities. The existing theories attempting to resolve problems of motor redundancy and variability can be broadly categorized into two lines: the principle of optimality and the principle of abundance. These two categories evolved from relatively separate views initially then converge as common constraints such as noise has to be considered. Principle of abundance The principle of abundance is a hypothesis proposed by (Gelfand and Latash 2002, 1998), where they stated that all the DOFs at all levels of a redundant system always participate in all the tasks ensuring both stability and flexibility of the performance. Contrarily to treating the extra DOF as a computational problems for the CNS, they view the apparently redundant design of the system as “a luxury that allows the controller to ensure both stability of important performance variables and flexibility of patterns to deal with other task components and possible perturbations” (Latash et al. 2007). This line of research is evoked 19 directly from the observation of structured variability as described before in this section. The essential goal is to form the synergistic actions of elementary variables and proper co-variation patterns so that the variability or noise in the redundant elements is largely projected into the uncontrolled manifold of the high-dimension state space. This subspace, in the linear algebra analysis, is essentially the null space of the Jacobian matrix (J(u)) that relates the small deviation of the high-dimension elemental variables (u=[u i ], i=1,2,…,m) to the error in the task performance space with lower dimension (p(u)=[p j (u), j=1,2.,…,n]) . To get co-variation that is beneficial for stabilization of the task variables, the CNS may be aware (Goodman and Latash 2006) or quickly learn [through a neuronal network (Bullock et al. 1993)] this coordinate transformation, in effect developing a forward mechanism that selects one solution out of a set of “motor- equivalent” solutions. Another physiologically based model, “central back-coupling model” (Latash et al. 2005) is based on a plausible neural mechanism of self- and lateral inhibition (similar to the system of Renshaw cells) among the elemental neurons. The tunable back-coupling loop in the neural model is a means of assuring error compensation in a task-specific manner. Both the forward and short-latency back- coupling models led to the idea of functional synergies, or high-level ‘control knobs’, that have invariant and predictable effects on the task-relevant movement parameters despite variability in individual DOFs. But how the synergies appropriate for a novel task, varied plant or a changing environment can be constructed, or adjusted/learned, and why the motor system should prefer such control mechanisms remain unclear. We need a computational framework that underpins these rather “descriptive” models with a “prescriptive” control scheme. A successful and influential work directed towards this goal is the theory of optimal feedback control (Todorov 2004; Todorov and Jordan 2002) that we will introduce in the following text and detailed in next section (I.iv). Principle of optimality As we have observed in the rich collections of motor functions, the existence of the stereotypical behaviors within and between species suggests that the motor evolution or motor learning must adopt a common rule of thumb – optimization. These approaches to the problem of motor redundancy involve application of optimization principles based on minimizing certain mechanical, engineering, psychological 20 or composite cost functions (detailed review see next section I.iv). It has been a common framework employed in engineering (Bryson and Ho 1975) as a foundation for the prescription of natural or synthetic motion control. Essentially the method is transferring the locus of parameter tuning from plant loop parameters to the chosen performance criteria or the cost functions; these choices hence largely determine the quality of the resulting solution (Holmes et al. 2006). Inversely, researchers tend to evaluate the “truthfulness” of a performance criterion for being what our body really wants to optimize by how well it predicts the average behavior. Use the example of multi-joint reaching movement, such cost function can be “jerk to reproduce the smoothed trajectory of human hand movements (Flash and Hogan 1985). This smoothness criterion has been successful in predicting straight paths, bell-shaped speed profiles of reaching movements, as well as a number of trajectory features in “via-point” tasks, thus could apparently be concluded as a “true” performance index that governs the sensorimotor control. However, the obvious ‘self-coding’ of the average behavior property (smoothness) into the parameter tuning function (cost) leaves the validity of this conclusion dubious. Furthermore, there has been no principled explanation why the nervous system should have evolved to optimize smoothness – are there fundamental physical reasons of why the CNS chooses to penalize jerk? The smoothness optimization has also been formulated on the level of dynamics, such as minimum torque change (Kawato et al. 1990; Uno et al. 1989), minimum command change. And similar to the minimum jerk model, the links of the inherent property of the neuro- musculo-skeletal system to the underlying determinant of movement control are lacking. One candidate for a fundamental (as opposed to arbitrary) determinant is signal dependent motor noise (SDN), a distinct feature of biological motor system as we reviewed early in this section. In presence of SDN, emergent properties of average behaviors in saccade eye and arm movements can be produced through minimizing the variance of position (‘endpoint variance’) (Harris and Wolpert 1998). Note that SDN in muscle force can be also be termed as “control-dependent” since the CNS “controls” the skeletal plant through directly commanding muscles; so the choice of control signals affects movement variability. This form of optimization explains the smoothness criteria including minimum jerk, torque change and command change, since any abrupt change in the trajectory will require larger muscle forces that induce larger noise level. Thus, SDN inherently imposes a “physical reason” behind the optimality principle. 21 However, the minimum variance model was critically dependent on the assumption of SDN arising during the movement execution. As we have mentioned, empirical studies have identified considerable amount of noises from central preparatory phase contributing the movement variability (van Beers 2007). This noise will presumably be larger for more difficult and/or less-practiced task (Churchland et al. 2006a). This may force us to reconsider those conclusions based on SDN assumption. For instance, “reaches may exhibit smooth velocity profiles to minimize the impact of signal-dependent noise. Alternately, smooth movements may be “easier” to plan and may incur less variability at the motor preparation stage” (Churchland et al. 2006a). An advance in this line is the theory of optimal feedback control by (Todorov 2004; Todorov and Jordan 2002) which accounts for sensory or other noise sources (Faisal et al. 2008), and incorporated variability in motor preparation. Most aforementioned theoretical frameworks along the line of principle of abundance and principle of optimality have separately emphasized the goal achievement and the richness of motor variability, but failed to reconcile the two. The one exception is the optimal feedback control theory proposed by (Todorov 2004; Todorov and Jordan 2002), where they showed that the optimal strategy in the face of uncertainty/noise is to allow variability in task-irrelevant dimensions. iv. OPTIMALITY PRINCIPLE IN MOTOR CONTROL Humans have a complex body with more degrees of freedom than absolutely necessary to perform any particular task. The redundant design of the human sensorimotor system provides flexible and adaptable motor behavior, but at the same time requires well designed controller that can intelligently choose among almost infinite number of possibilities. Optimization theory has become an important aid to discovering organizing principles that guide the generation of goal-directed motor behavior, specifying the results of the underlying neural computations without requiring specific details of the way those computations are carried out. Different computational theories can be obtained by varying the specification of the physical plant controlled, the performance index optimized, and the control constraints imposed; and the application of optimal control framework forms two types of control laws: open-loop (feed-forward) and closed-loop (with feedback) formulations. In this section, we will first delineate the importance and the general 22 formulation of the optimization theory; then the evolvement of the cost function and performance indices in the literature will be reviewed; this is followed by a brief review of the open- or closed-loop optimal control laws; and finally further considerations on optimal control will be presented. a) What is and why need optimization? The optimality principle forms the basis of many physical sciences (Bryson and Ho 1975). The application of it is advantageous both in its natural origin and its power of solving ill-posed problems. In its origin, the concept obeys a process of “natural selection”, where organisms evolve to select favorable traits that satisfy the demands of life. In the context of motor control, the process of sensorimotor integration is to adjust continuously the coordination of the redundant states and effectors so as to improve the behavioral performance. Even if the skilled performance is not necessarily optimal, a good performance is sought in a process with asymptotic optimal limit. Thus optimality provides a natural starting point for computational investigations of sensorimotor function (Todorov 2004). Practically, it also provides a powerful tool in searching among the redundant solutions the best possible coordination that yields the desired performance. The method is autonomous and general, and forms the basis of various control schemes (open- or closed-loop), motor learning process and adaptive control techniques. As we have reviewed in Section I.ii.b, many alternative methods solving redundancy in engineering control and neural control of movement reduce the number of independently controlled variables by forming synergies favorable for a given task. These synergies have been mostly observed empirically in the particular task under study. For example, in a the multi-finger force production task (Latash et al. 2005), the individual finger forces are formed by a sharing pattern that is in favor of stabilizing the total amount of pressing force (Li et al. 1998; Zatsiorsky et al. 2000), but how this pattern would be reorganized with a new task, for example, maintaining the position accuracy, was not addressed. A novel task usually demands a new set of synergy, the learning of which presents the CNS a problem with the same dimensionality as the original control problem. The optimal control method provides a computational underpinning of this evolvement of the synergy: by just identifying a performance criterion that describes what the goal is, the details of the synergistic actions of redundant elements will evolve automatically. 23 Generally, to apply optimality principle in computational motor control, the modeler has to specify (i) a family of admissible control laws (open- or closed-loop), (ii) a musculo-skeletal dynamic model (f(x,u, θ)), describing the plant dynamics with specific parameters ( θ) that relate the control (u∈ℜ n , such as muscle activation) to its system state (x∈ℜ m , such as, muscle force, joint angle, velocity) and performance outputs (y(u,x) ∈ℜ s , such as hand positions or finger force); and (iii) a quantitative definition of task performance (or a value function J∈ℜ) which is time-integral of some cost function (cf∈ℜ, such as energy, jerk, etc.). The cost is an instantaneously defined scalar function that depends on the current set of control signals and the set of variables describing the current state of the musculo-skeletal system. Motor redundancy is reflected in the larger dimensionality of the control space than the state space, and of the state space compared to the task space (i.e.n>m>s). The goal of optimization is to find a unique control policy, u x (e.g. a sequence of u or synergy of u among many possibilities of u), that minimizes the time integral of the cost function or maximizes a performance index (e.g. position accuracy). In a symbolic form, the optimization problem can be formulated as: J min u x where () ∫ = T 0 dt cf J (1.1) Subject to () θ = , , f u x x & and 0 ≤u i ≤1, i = 1, …, n. Notice this control tuning process relies heavily on the cost function (cf), while it ignores the detailed plant loop parameters ( θ). The quality of the solution u x is thus determined by how good the cost function or the performance index is. b) What constitute performance criteria? In human motor control studies, “minimum X” is often a signature of the optimality models used in searching for the control policies. The selection of a motor control is governed by this “X”, where X can be energy, muscle or bone stress/strain levels, speed/dexterity requirements, or other measures of performance such as smoothness, and sometimes a composite cost specified by a weighted set of different performance criteria. Costs, in the essence of optimality theory, denote what the sensorimotor system is trying to achieve, 24 so an appropriate quantification is the most important step in applying the optimal control model. But it is not always an easy task. A rare case where the choice of cost is transparent are behaviors that require maximal effort – for example, jumping as high as possible (Pandy and Zajac 1991; Pandy et al. 1990). Here one simply optimizes whatever subjects were asked to optimize, and the model predicts good fit with data. In most cases, however, the intuitive characteristics of the task do not necessarily reveal the true cost that the sensorimotor system actually cares. They are often hidden among the various behavioral regularities (smoothness of hand path), numerous constraints of the physical plant (such as the stability criteria), or the inherent properties of the signals (such as sensory or motor noise). In such cases, the detailed form of the cost is often considered a relatively free parameter, the choice of which has attracted a lot of attention in the literature. In a lack of empirically derived costs, researchers experimented with a variety of definitions, and found that motor behavior is near-optimal with respect to costs that vary with the task. These costs can be subdivided into four general categories – energy, smoothness, variance and multi-attribute criteria. Energy cost Optimal control modeling has a long history in biomechanics, particularly locomotion, where most models minimize energy used by the muscles (Anderson and Pandy 2001; Collins 1995; Davy and Audu 1987; Hatze and Buys 1977; Rasmussen et al. 2001). Although precise models of metabolic energy consumption that reflect the details of muscle physiology are rare, a number of cost functions that increase supra-linearly with muscle activation yield realistic and generally similar predictions (Collins 1995). Some models use optimization with additional behavioral measures as their limitation: they measured the kinematics of the gait cycle and compute most efficient muscle activations or joint torques that could cause the observed kinematics. These costs can be total muscle force, squared muscle force, muscle stress, intra-articular contact force and instantaneous muscle power. As an example, several attempts (Rasmussen et al. 2001) have been made to compute muscle force with various polynomial forms of cost function: p 1 i pcsa F J ∑ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = n i i m (1.2) 25 where n is the number of muscles, i m F is the force produced by i th muscle, and p is the power of the polynomial form that has been experimented with different p (1, 2 or 3) (Crowninshield and Brand 1981; Crowninshield et al. 1978; Herzog 1987; Pedotti et al. 1978). Other models avoiding this limitation lead to more challenging dynamic optimization problems. For example a more recent work (Anderson and Pandy 2001) studying human ground walking, a body with 23 mechanical DOFs and 54 muscles was used for estimating the optimal walking kinematics, muscle force and EMG sequences, with empirical measures only specifying the initial and final posture of the gait cycle. The cost function is the total metabolic energy expenditure per unit distance traveled. This performance index follows a more complex form that includes whole body metabolic heat, muscle energy consumptions depending on the internal state and mechanical work rate [details see reference (Anderson and Pandy 2001)]. Smoothness cost Though rather successful in biomechanical and locomotion applications, energy minimization alone has failed to account for average behavior in arm movements (Nelson 1983), eye movements (Harwood et al. 1999) and some full-body movements such as standing from a chair (Pandy et al. 1995). The usual remedy in such cases is a smoothness cost, which penalizes the time-derivative of hand acceleration (or jerk) (Flash and Hogan 1985; Hogan 1984; Smeets and Brenner 1999), joint torque changes (Kawato et al. 1990; Nakano et al. 1999; Uno et al. 1989), or muscle force (Pandy et al. 1995). The smoothness optimization has been successful in predicting average trajectories, particularly in arm movement tasks. The idea was first introduced in the minimum-jerk model by (Hogan 1984). Such kind of cost function depends on the kinematics of the motion, and the variables of interest include the positions (e.g. joint angles or hand Cartesian coordinates), angular velocities, accelerations and higher derivatives (e.g. jerk). With an example of plannar two-joint arm movement, x, y denote the Cartesian coordinates of the hand position, the minimum jerk cost function is: () ∫ + = T o dt y x 2 1 J 2 2 & & & & & & (1.3) 26 where T is the duration of the movement, and jerk is the third order time derivative of position. By assuming zero initial and ending velocities and accelerations, the optimal trajectory predicted from the minimum jerk model accounts for the straight paths and bell-shaped speed profiles of reaching movements, as well as a number of trajectory features in ‘via-point’ tasks. This model has been extended to the task of grasping, by using a cost that includes the smoothness of each fingertip trajectory plus a term that encourages perpendicular approach of the finger tips to the object surface (Smeets and Brenner 1999). The optimization result explains many observations from grasping experiments, particularly the effects of object size on hand opening and timing. Smoothness optimization has also been formulated on the level of dynamics – by minimizing the time- derivative of joint torque (Kawato et al. 1990; Nakano et al. 1999; Uno et al. 1989). This formulation arose from observation of mild curvature in hand-reaching trajectories in a large range of movements. Such kind of cost function depends on the dynamics of the arm, and the variables of interest contain torque change, muscle tension and muscle command. Let dτ ι /dt denote the rate of change of torque at the i th joint, the cost function of minimum torque-change model is: ∫ ∑ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ τ = = T o n i i dt d dt 2 1 J 1 (1.4) Interestingly, despite the nonlinear dynamics of the arm, the hand-reaching trajectories predicted by this model are roughly straight in Cartesian space, with a mild curvature that agrees with experimental data. Variance cost Although the smoothness criteria have predicted average behavior that achieves the task goals accurately, but they can not account for the variability observed among many individual trials. Also, they only reproduced the smooth trajectories without revealing the advantage of producing smooth movement, and neither did it explain why the sensorimotor system evolved to optimize such a quantity. Furthermore, even if jerk or torque change is really what the nervous system cares to minimize, the neural mechanisms that could performance such complex estimation and integrate them over the movement duration remains unknown. Within the open-loop optimization, these issues were addressed by the minimum variance model – TOPS (Task Optimization in the Presence of Signal-dependent noise) (Harris and Wolpert 1998). This 27 model predicted the sequence of muscle activations (u) and movement trajectories by minimizing the resulting variance of final hand/eye positions. ()( ) [] 2 T 2 T y* y x* x 2 1 J − + − = (1.5) where (x T, y T ) denotes the final hand position, (x* , y*) is the target position. The TOPS model is based on the assumption that the neural control signals, u, are corrupted by noise whose variance increased with the size of the control signal. This SDN imposes a fundamental physical reason for the smoothness of the movement trajectories: non-smooth movements require abrupt changes of muscle force, which require larger activation levels, hence more noise. It also explains the speed-accuracy trade-off (Fitts’ law): fast movement requires large muscle forces that would carry more noise, and the elevated control noise will cause the trajectories to deviate from the desired path. These deviations will accumulate over the movement time and result in large errors in the final position. Accuracy could be improved by having low amplitude control signals, but the movement will be slower accordingly. In addition, the minimum variance model partially predicts the pattern of variability. Furthermore, the variance of the final position reflects the intuitive cost of the CNS in reaching tasks, and its quantification is directly available to the CNS (for example through visual feedback). Multi-attribute cost (mixed error/effort cost) Exploration of simple costs has illuminated the performance criteria relevant to different tasks. However, the true performance criteria in most cases are likely to involve a mix of cost terms. This has been manifested by the empirical findings in the single-joint elbow movements (Osu et al. 2004): Without the need for greater accuracy, subjects accepted worse performance. This contradicts the minimum variance principle of TOPS and suggests that minimizing the endpoint deviation or endpoint error alone may not be the sole consideration of task optimization. It is clear that there are extra dimensions to the accuracy criteria and energetic cost has been proven to factor into many tasks (Burdet et al. 2001; Osu et al. 2004). In an extended TOPS model [TOPS-α (Nagata et al. 2002)], a cost term of the motor command magnitude multiplied by weighting factor, α , was added to the task achievement term (variance) and gave better prediction of trajectories than the TOPS model for a large number of point-to-point movements. The 28 similar mixed criteria also predicted muscle directional tuning – how the activation of individual muscles varies with the direction of desired net muscle force (Todorov 2002). While multi-attribute cost clearly improves the predictions of the open-loop optimization model, this mixed formulation is necessary for the framework of stochastic optimal feedback control model (Todorov 2004, 2005; Todorov and Jordan 2002). This is because the boundary constraints (such as hand position, velocity and acceleration at the beginning and end of the movement) are inapplicable to stochastic problems in which the final state is affected by noise. The stochastic models have to use final accuracy costs in addition to whatever costs are defined during the movement (Todorov and Jordan 2002), mixed criteria similar to TOPS-α (Nagata et al. 2002). Though assigning the relative weights (α) to quantities that have different units (energy and position error) is a non-trivial task, it may reveal a higher level command signals that consciously specifies the relative importance of different aspects of motor behaviors. This is closely akin to the real biological motor construction and has been implemented in the later optimality models that have hierarchical structures c) Open-loop and closed-loop control In terms of the control structure, optimization models of movement control can be classified into open- loop and closed-loop. Open-loop optimization essentially is to predict average movement trajectories or muscle activity by optimizing a variety of cost functions (as we have reviewed before), thus can be also called “models of average behavior”; whereas closed-loop optimization combines various information from different sensory modalities to optimize the performance, thus can be named as “models of sensorimotor integration” (Ijspeert 2001; Todorov 2004; Todorov and Jordan 2002). Open-loop construction has been a powerful aid to discovering organizing principles of motor control and specifying the results of the underlying neural computations, but a closed-loop approach is more akin to the biological sensorimotor organization and can address many more phenomena that the open-loop models have ignored or avoided. One such feature is the sensory and motor noise, the treatment of which basically differentiates the two control structures. 29 An open-loop optimal controller produces the movement following a sequential process: firstly movement details (such as muscle activation sequence or movement trajectories) are planed by optimizing predefined performance criteria, and then this “desired” plan is executed by a servo mechanism which cancels the instantaneous deviations between the desired and actual state of the body. This servo controller, however, has a fixed structure and is not taken into consideration in the optimization phase. The key principle of open-loop control is that it depends on the ability of the CNS to predict the consequence of the perturbation and assume there is a separation between the movement planning and execution. In contrast, closed-loop optimization treats the feedback mechanism as being fully programmable, which is similar to the concept of adaptive control, where the controller parameters are constantly tuned toward the optimizing control objective (Astrom and Wittenmark 1995; Tin and Poon 2005). In another word, it constructs the best possible transformation from states (x) of the body and environment into control signals (u(x)), and the resulting controller does whatever is needed to accomplish the task, that is, control scheme independent (force, position or impedance control). Such flexibility, however hard to grasp, matches the flexibility and resourcefulness apparent in biological motor behaviors (Liu and Todorov 2007; Todorov 2004). Note, this optimal tuning process relies on mapping actual body states into control signals, but when the state of a stochastic plant is observable only through delayed and noisy sensors, which is the case for biological motor system, the controller has to rely on an internal estimate of the state. The optimal feedback control framework addresses this problem by integrating an optimal controller with an optimal estimator (based on Bayesian inference) (Kording and Wolpert 2006; Todorov and Jordan 2002). In modeling practice, one typically uses Kalman filter, an optimal estimator when the dynamics and sensory measurements are linear and the noise is Gaussian. The estimator takes into account both sensory data, motor commands, and applies knowledge of body dynamics as well as its earlier movement output, and weights all these sources of information regarding the current state in proportion to their reliability. In another word, an integral component of feedback control is the presence of a state estimator, i.e., a process that combines sensory inflow and motor outflow to provide an estimate of the current state of the system [for detailed review please refer to (Kording and Wolpert 2006; Wolpert et al. 1995)]. 30 The results produced by optimal feedback controllers for such system match strikingly well with the complex behavior observed in human motor tasks. With an error/effort mixed cost, an optimal feedback controller (Todorov and Jordan 2002) predicted not only the invariant properties of multi-joint arm reaching movements, accounted for speed-accuracy trade-off, skewed velocity profiles, but also the structure of motor variability. It constructed a unifying theory that brings together a number of key observations and ideas in motor coordination study of redundant, noisy motor system. In a redundant motor system, the abundance of solutions is beneficial because it affords the nervous system both the ability ito stabilize the important performance variables and flexibility to compensate the noises or unexpected changes in the environment (see the review on structured variability in Section I.iii.c ). But this abundance also presents a difficulty in understanding how the choice is made. Optimal feedback control provided a principled solution in resolving the redundancy by prescribing a ‘minimal intervention’ principle: make no effort to correct deviations away from the average behavior unless those deviations interfere with task performance (Todorov and Jordan 2002). d) Further considerations in optimal control Optimality models have provided a principled approach that gives a unified account of many complex motor behaviors. Sophisticated control structure (such as optimal feedback control) and refined performance criteria (such as mixed error/effort cost) have substantially enriched the optimality models. More considerations are still needed to further the understanding of motor control. The first consideration is automation of error/effort cost. The error/effort cost function contains parameters that weight the respective contribution of state errors and effort in the control. Since the two components in the cost function have different units, the modeler has a non-trivial task of assigning the relative weights according to the intuitive importance in different aspect of the task. And, because different sets of parameters lead to different behaviors, a model based on error/effort minimization cannot provide an unequivocal description of motor control. There have been efforts in developing algorithms to automate this trade-off by defining a global optimality (Zhu et al. 1997). But how this optimality evolves with tasks, within novel environments and under different constraints is still not clear. There may be leading criteria 31 that vary with the goal of the task, or are constrained by inherent properties and sensorimotor process such as noise, or dominated mechanical interaction with instable environment (stability criteria). The second consideration involves the separate control of movement and posture. In fact, the central difficulty for an optimal control approach to redundancy is to formulate the problem in such a way that it accounts for the simultaneous control of posture and movement (Guigon et al. 2007). Here the movement control requires a dynamic controller that produces desired movement by generating sequences of muscle activity counteracting limb inertia and interaction torque; whereas the posture control deals with the static forces such as gravitational or muscle elastic force. So far, most optimal control theories have not considered the case of static forces, and a few that tackled it have reported difficulties in solving optimal control problems in the presence of gravitational forces (Thoroughman and Feller 2003). Although the neural bases for separate movement and posture controller are still elusive, behavioral observations of (Brown et al. 2003) clearly suggest their existence. A recent computational effort (Guigon et al. 2007) combined the optimal feedback controller with a separation principle (i.e. posture/movement separate control) has successfully predicted various characteristic features of pointing and grasping movements in a 3D space. However, the posture controller so far in the literature (Guigon et al. 2007) is often a simple mapping between states and forces; it cannot directly contribute to the stability of the equilibrium. In a system where noise is inherent in both sensory encoding and motor signals, instability is an ever- lasting problem that the CNS has to combat. This leads to the notion of impedance control, i.e., a controller that exploits the viscoelastic properties of muscles to maintain stable equilibrium. We will review impedance control and stability in next section. v. IMPEDANCE CONTROL IN BIOLOGICAL MOTOR SYSTEM Humans live by interacting: we push/pull a door to enter a house, we manipulate a tool to open a can, we pick up a cup to drink coffee, and we hold each other when we dance… stable interaction with the environment is crucial to achieve any behavior goals, but unfortunately, in our daily living, many tasks are inherently unstable (e.g. when we use a screwdriver), and this unstable nature is further amplified by noise 32 in our motor system. The neuromuscular system must compensate for such instability to successfully perform the task. One strategy is impedance control. In this section we will first define impedance control and outline its importance; then several major empirical items of evidence for impedance modulation in human movement will be reviewed; finally we will point out impedance-noise interplay, a unique phenomenon in biological motor system. 33 a) What is and why need impedance control? To clearly explain the importance of impedance control, we will first need to define impedance and the general concept of impedance control. Mechanical impedance is a measure of opposition to motion of a structure subject to force. It can be attributed to a mechanical system characterized by mass/inertia (M), viscosity/damping (Β), and stiffness (K), and the combined feature determines the dynamic responses of 1 q 2 q Shoulder Elbow l 1 l 2 We tend to understand the kinematic effects of muscle-joint geometry by considering how the mechanical behavior of the muscles is transformed into an equivalent behavior in the intro-personal joint space and then in the extra-personal hand space. These transformations are strongly dependent on the muscle-joint connection and joint configurations. As muscle provided arm stiffness has been a common focus in motor control study and most accessible in the experimental measurements, we now describe how the muscle elastic behavior is transformed physically and mathematically into the joint and hand spaces. Figure on the left illustrates a 2-link 2-DOF arm actuated by 6 muscles across shoulder, elbow and both joints. q=[q 1, q 2 ] T is the a vector of generalized coordinates, and l 1 and l 2 are the lengths of the upper and forearm respectively. Given the muscle stiffness matrix: BOX.4| Transformation theory for stiffness – muscle to joint to hand space (Hogan, 1991) [ ] 6 5 4 3 2 1 m m m m m m K K K K K K diag = m K we can obtain joint stiffness K j by the following equation: () () () m ee es se ss K K K K F q q R K q R K m T j ⋅ + ⋅ ⋅ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = Γ where R(q) is the moment arm matrix of the six muscles that is determined by the partial derivative of the muscle length (L mt (q)) with respective to joint angle: () ( ) q q L q R mt ∂ ∂ − = and Γ(q)  is an array of second partial derivatives of muscle lengths with respect to joint angles; it quantifies the variation of the moment arm: ( ) ( ) 2 2 q q L q Γ mt ∂ ∂ − = The second term of joint stiffness transformation (Γ(q) F m ) is a “fictitious” joint stiffness due to the variation of the moment arms with configuration. This component is usually negligent at small muscle forces (F m ), but becomes progressively more significant at high level co-contractions. Transforming the joint stiffness (K j ) to the end-point or hand stiffness (K h ) follows this relationship: () () () 1 j T 1 h q J K q J K − − ⋅ ⋅ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = yy yx xy xx K K K K here J(q) is the arm Jacobian matrix: () () ( ) ( ) () () () ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + + + + − + − − = 2 1 2 2 1 2 1 1 2 1 2 2 1 2 1 1 cos cos cos sin sin sin q q l q q l q l q q l q q l q l q J 34 the system to any imposed force (generalized force), either control commands, external load, or environmental perturbations. In a multi-DOF musculoskeletal system, force and motion are vector quantities with direction as well as magnitude. The three impedance quantities will become matrix formation, and they also have directional properties. In particular, the force and displacement need not be co-linear; both the magnitude and the direction of the output force vector resulting from an input displacement vector of a given magnitude will depend on the orientation of the displacement vector (Hogan 1991). In motor control literature, there is a common misnomer in the term “impedance”, which has been abused when the work was actually referring to “stiffness” measure. Nonetheless, this misuse has its “biological” roots: Compared to the passive joint linkage such as ligament and cartilage, skeletal muscles have dominant contribution to the viscous and compliant property of the musculo-skeletal plant, and they are actively controlled by the nervous system, while the limb mass/inertia properties are mostly constant or not controllable independently from posture. Nevertheless, a strict definition of “impedance control” also involves mass/inertia, whereas in biological motor control literature, “impedance control” is presumably referring to viscosity and stiffness. Impedance control basically is to allow complete characterization of system dynamics through active impedance specified by the three matrices Ma Βa and Ka. It is one of the interaction control schemes in robotic literature (Sciavicco and Sicilliano) where manipulator has to interact with the environment. In such cases, the environment with passive impedance, Mp Bp and Κp will be integrated with the manipulator as a combined mechanical system constituted by the parallel of the two, and then its dynamic behavior is conditioned by the relative weight between them. In any mechanical structure, the impedance property determines the inherent stability state of the system, thus, impedance control is a control scheme closely akin to successful performance of the task with high stability requirements. This is the case for most human activities. Many manipulation tasks in our daily activities are inherently unstable (Rancourt and Hogan 2001). For example, when using a screwdriver, fluctuation of the force produced by muscles could cause the screwdriver to slip, creating an unstable interaction that may result in large errors and failed action. 35 Stability, thus, becomes one major constraint or a common goal of many theoretical control frameworks. For instance, optimal feedback control theory (Todorov and Jordan 2002) proposes an online correction mechanism to pull back the deviated limb from hurting the performance. But the long delay of neural feedback make the error correction mechanism itself possible to introduce more unstable factors (delayed negative feedback loop). Fortunately, muscles has preflex structure whose inherent force-length and force- velocity property creates an intrinsic response to perturbation with zero delay (Brown and Loeb 2000c). Muscles are specialized actuators with tunable stiffness and viscosity roughly proportional to the output force. So while the tasks may be unstable, this active impedance of the musculoskeletal system is stabilizing in nature. In fact, modulating mechanical impedance is a general strategy for controlling interaction between both a limb and its environment (Hogan 1985) and a flexible base of support and the limb (Andersson and Winters 1991). Due to its immediate response to unexpected perturbations, the preflex property and impedance control has been proposed as an integral component of hierarchical optimal controller (Todorov 2004). Considering from a different perspective, impedance regulation that is needed to satisfy stability in a motor task essentially is placing another dimension to the output space of motor task (force, position and stability), thus reduces the neuromuscular redundancy. Taken a simple case, a single joint kinematics or joint torque can be determined by activating one muscle (i.e. one motor) around it, a torque regulation task then is non-redundant with one-muscle, one-joint system. But since muscles can only pull not push, to produce torque or movement in both directions (flexion/extension) will need a pair of muscles opposing each other (antagonistic muscles). Because the impedance contributions of the muscle add while their moment contributions subtract, by carefully coordinating muscle activities, the net joint torque may be held constant while the net impedance increases. In this case, a joint torque can be produced by many possible activation levels of the two muscles, thus arises a problem neuromuscular redundancy. However, if the task requires a necessary impedance to meet its stability criteria, an appropriate co-contraction level will result unique solution of this apparent redundant system. A recent modeling study (Bunderson et al. 2008) on neuromuscular redundancy and postural stability provided computational evidence. The authors used a three-dimensional musculoskeletal model of the cat hind limb with 31 muscles to determine the muscle 36 contribution to limb stability during isometric force generation. The dynamic stability analysis demonstrated that within the large set of activation patterns that satisfy the force requirement for posture, only a reduced subset produced a mechanically stable limb configuration. This implies that muscular redundancy is reduced when muscle activation patterns are chosen with respect to intrinsic musculoskeletal stability as well as endpoint force production. Finding the solution may sound easy in the single-joint case, but certainly non-trivial for multi-joint system. The kinematics of the musculo-skeletal system is far more complex with many DOFs and redundant muscles. Interaction between different links in a multi-joint mechanics added directional property to endpoint impedance, and the geometry of the system (e.g. shoulder and elbow angles) can be a dominant factor in determining the spatial pattern of the endpoint impedance (Mussa-Ivaldi et al. 1985). From control point of view, modulating the directional character of the impedance is important for performing interactive tasks. For example, the limb might be made compliant in one direction to accommodate an external kinematic constraint and stiff in another direction to minimize the effects of disturbing forces (Hogan 1991). The desired changes in the orientation of the endpoint stiffness can theoretically be produced by coordinating muscle groups around the joints. Then the question is can neuromuscular coordination override the joint geometry to achieve such flexible modulation, and if it can, does the CNS exploit this capability to control impedance? Both conditions appear to be true, and the empirical evidences will be reviewed below. b) Has impedance control scheme been seen in human movements? Although the history of musculo-skeletal mechanical impedance provide by neuromuscular coordination is not so long [since (Hogan 1985)], the appeal of maintaining system stability through impedance modulation have attracted huge amount of research across many labs in the world of motor control (Burdet et al. 2001; Burdet et al. 2000; Burdet et al. 2006; Flash and Mussa-Ivaldi 1990; Franklin et al. 2007; Franklin et al. 2004; Hoffer and Andreassen 1981; Hogan 1985; Hunter and Kearney 1982; Ostry and Feldman 2003; Osu et al. 2003; Osu et al. 2002; Osu and Gomi 1999; Osu et al. 2004; Perreault et al. 2002; Selen et al. 2005; Selen et al. 2006). In fact, since the first observation (Mussa-Ivaldi et al. 1985) that showed invariant property of endpoint stiffness determined by arm configuration (shoulder and elbow joint 37 angles), it took a long while (about 15 years) for the numerous contradicting evidences to emerge – e.g. isometric force regulation (Gomi and Osu 1998; Perreault et al. 2001, 2002), movement (Burdet et al. 2001; Burdet et al. 2006; Franklin et al. 2007; Franklin et al. 2004), etc. – that proved the flexible modulation of impedance patterns by neuromuscular system. In a two-link limb endpoint force regulation task, (Perreault et al. 2001) found that arm stiffness is nearly posture-independent, and were due primarily to the actions of the single-joint muscles spanning the shoulder and elbow; a later study (Perreault et al. 2002) further demonstrated that stiffness orientation can be modulated if there is independent control of the two joints. To examine if the impedance modulation also occurs during movements, (Burdet et al. 2001) placed subject’s hand in an unstable environment characterized by a divergent spring force field along the axis perpendicular to the movement direction; the resulting stiffness measure demonstrated a selective increase in the direction of instability; and the increase in endpoint stiffness was interpreted as occurring primarily through the coactivation of biarticular muscles (Franklin et al. 2004). However, if a general impedance controller exists, the CNS should be able to match changes in endpoint stiffness to arbitrary changes in the orientation of the instability. A recent work (Franklin et al. 2007) provided such an evidence by showing subjects being able to increase the hand 200 N/m S1 S2 S3 S4 S5 S6 S7 S8 0 o -45 o 80 o DF DF BOX.4| Multi-joint stiffness adaptation to instable environment (Franklin et. al. 2007) End-point stiffness ellipses has been shown to be able to adapt to the environment with instable divergent force (DF) field perpendicular to the hand moving direction (Burdet et al, 2001). This study by Franklin et al. (2007) further confirmed the flexible control of limb impedance properties by examining the stiffness field subject to DF along arbitrarily selected directions (0 o , 80 o , and -45 o ). The figure demonstrates that the eight subjects all showed a primary increase in the stiffness along the direction of DF, whereas in the orthogonal direction, the magnitude was close to the null-force field condition (gray filled ellipses). Although not perfectly aligned with the DF direction, the ellipses are significantly rotated toward the direction of instability. These results confirmed that the CNS is able to control the arm impedance properties and selectively adapt them to the environment. The flexible modulation is likely achieved through coordinating the multiple muscles around the joints, and thus ensures stability with minimal energy cost. 38 stiffness selectively along the different directions of unstable force fields (see box 5). Furthermore, the EMG recording of single-joint and bi-articular muscles showed different co-activation patterns associated with the three adapted stiffnesses. So far, however, impedance modulation in multi-joint multi-muscle systems has been shown to stabilize the musculo-skeletal system in response to external perturbations. But biological motor system is full of internal fluctuations (e.g. motor noise). When desired impedance is set by a specific muscle coactivation pattern, the presence of motor noise will shift the pre-specified value. So, as we will argue in the following section, an impedance controller only concerning environment dynamics may be suboptimal. c) How do impedance control and signal-dependent noise interplay? Noise inherent in the sensorimotor system has provided constraints for a broad range of control schemes (reviewed in Section I.iii.c), impedance control is no exception. Impedance modulation has been shown to be an effective control method used by the CNS to gain stability in instable environments, but its control on internal destabilizing perturbations seems paradoxical. This is because, increased muscles 0 0.5 1 0 2 4 Co-activation level Kinematics variability (deg) BOX.5| Co-contraction in single-joint postural control with signal-dependent noise (Selen et al 2005) To combat environmental instability, co-contracting the antagonist muscles to increase the limb impedance is a commonly adopted strategy by the CNS; but in the presence of SDN, the increased impedance is also accompanied by the increased muscle force fluctuations. To stabilize a joint posture against SDN, we thus have two competing factors that create a paradoxical picture in human motor system. For movements, the limb dynamics and external perturbations might dominate the SDN in affecting limb stability and movement accuracy, but in postural tasks or many fine motor behaviors such as needling, SDN can have severe effect on the motor performance. Selen et al (2005) took the first step in examining this noise-impedance interplay through a simulation study on a simple single-joint system. Here the model is symmetric in that the two opposing muscles are identical with same architectural parameters, and same constant moment arms. The task was to keep the joint at a specified equilibrium angel while varying the co-activation levels of a pair of opposing muscles with SDN introduced in the muscle command signals (CV=0.2). As shown in the figure, the simulation results demonstrate a general trend of reduced joint variability with increased co-contraction level. Only at the lowest levels of co-activation, higher force resulted in larger joint fluctuation (solid line). By manipulating the noise amplitude and joint inertia parameters, the author was able to see the mono-tonic relationship along the whole range of co-activation (dashed line). These results support the idea that muscle co-activation and stiffened joint is effective for control SDN. 39 activity not only provides higher level impedance to damp the kinematic effects of noise, but it also is accompanied by larger fluctuations of generated force (i.e. SDN). We then have two competing mechanisms that determine the level of kinematic variability; which one, at what point dominates may determine what impedance control policy the CNS may choose. For instance, as muscle co-activation increases, if the increased impedance always surpass the effects of increased motor noise (i.e. a monotonic relationship between kinematic variability), then to achieve a maximum accuracy, a maximum co- contraction will be a desired policy. Whereas, if the relation is not monotonic but a convex curve, then there will be an optimal impedance with respect to minimum variance. Observations in single-joint (Osu et al. 2004; Selen et al. 2006) and multi-joint (Gribble et al. 2003) movements have demonstrated increased muscle co-activation when increased accuracy is demanded. A simulation study that took into account the stochastic nature of the motor unit recruitment process and the noise filtering properties of a biomechanical limb also suggested monotonic decreases of kinematic variability with co-activation. A recent modeling study (Selen et al. 2005) combined a motor-unit pool (Jones et al. 2002; Van Galen and De Jong 1995) with muscle dynamics (force-length and force-velocity properties) to model the opposing muscle actuators around a single joint, and the simulation results also demonstrated a monotonic trend, except for the lowest levels of co-activation. Although this monotonic relationship has been generally observed, the relation between variability and impedance could be far more complicated in complex systems with many DOFs. In a stochastic model of motor system, the biomechanical limb is essentially a noise filter, and the effect of the motor noise on the final endpoint kinematics is characterized by the dynamic property of this mechanical filter. Multi-joint mechanics results in complicated dynamics that depend on joint position (such as apparent inertia) and relative motion of multiple links (such as interaction torque). How the impedance could reduce the internally generated noise in such a system is still an open question. Furthermore, in a stochastic motor system, the limb impedance may be only vaguely defined because the reference trajectory, or the equilibrium position itself could be perturbed by internally generated noise. From these perspectives, we might postulate that the SDN inherent in neuromuscular systems indeed presents another constrant in addition to the stability criteria derived in a deterministic system 40 (Bunderson et al. 2008), thus further reducing the dimensionality of a neuromuscular redundant biomechanical system. Simulation results in our study provided support of this view. d) In view of optimal impedance control One drawback of impedance modulation strategy is that it is energetically inefficient. By opposing one another, the muscles perform no net mechanical work, yet they consume metabolic energy. “If we make the reasonable assumption that the biological system does not squander energy needlessly, one would expect synergistic activation of antagonists to be kept to a minimum except when the task demands it” (Hogan 1991). This trade-off has been experimentally observed in a single-joint goal directed arm movements (Osu et al. 2004): Without need for great accuracy, subjects accepted worse performance with lower co- contraction. This view of optimal impedance control is obeying the principle of maximum energy- efficiency, and the choice of optimal impedance is rather straightforward: find the minimum co-contraction level that corresponds to the largest allowed levels of kinematic variability (minimum accuracy). Clearly this is based on an assumption of a monotonic accuracy-impedance relationship. And the single-joint target reaching tasks does not involve any unstable interaction with the environment. In the human motor activities of daily living, however, neither of these two assumptions seems valid. First of all, interacting with dynamical environment places stability criterion into the control, hence a specified value of impedance may be needed to meet this requirement. Secondly, in a complicated mechanical system with many DOFs, the uncertain noise-impedance interaction would not allow a monotonic translation of minimum energy to least accuracy. In a realistic setting, a search for the “optimal” impedance can actually follows multiple surfaces depending on what is the “priori” in the task: If stability dominates, there may be a single value of impedance with respect to plant-environment system; if kinematics accuracy is to be maximized, locus of the impedance value will correspond to the minimum hand variance on the variance-impedance surface; and if the task does not need least possible hand deviation, rather just good enough accuracy, then the optimal impedance will be associated with minimization of energy expenditure. Therefore, the goal is to adapt the impedance controller to a maximized performance according to the task, a concept quite similar to the process of “reoptimization” (Izawa et al. 2008): Through practice in the novel environment, we learn internal models that predict 41 sensory consequences of motor commands; through reward-based optimization, we use the internal model to search for a better movement plan to minimize implicit motor costs and maximize rewards. In a complicated structure like human sensorimotor system, many factors intertwine (e.g redundancy, noise, motor variability, impedance, stability, energy). No matter the goal is searching control theories of a motor function, or explaining important features of motor behaviors, one must always keep in mind that these factors coexist and often interplay; and depending on the task, some of them may constrain the proposed control theories (such as noises on the optimal control), and others may facilitate us to find a solution (such as stability on muscle redundancy). Although many attempts have aimed at finding a cohesive control theory that could unify multiple features of human motor behaviors, none of them has successfully captured the whole picture. The reasons are multiple: it can be the muscle model being used [e.g. simple muscle model limited optimal feedback control (Todorov and Jordan 2002)]; or limitation of experiments [e.g. instable environments and other sources of motor variability could mask the effect of signal-dependent noise on impedance modulation during reaching movements (Franklin et al. 2007)]; or the system under study [e.g., a single-joint simulation limited the understanding of motor noise effect in redundant limbs (Selen et al. 2005)]. In the work presented in this thesis, the author aims at capturing a larger subset of the full scene of noisy, redundant human motor systems, exploring how the multiple features combined shed light on computational theories of motor control, which are essentially statements about the brain’s strategies for adjudicating among these multiple features. A key means to avoid experimental constraints that may cause us to overlook some of these features is to build a realistic sensorimotor system model and define a well-posed task, as described below. vi. REALISIM IN SYSTEM-LEVEL MODEL AND SIMULATED TASK Recently, scientists in motor control study often face an increasingly intricate picture, where it is hard to uncover the fundamental underlying principles that bear the real explanatory power. The panorama consists of a huge number of interconnecting system states existing in complex musculoskeletal system, neurophysiology of brain and spinal cord. Many of these states are likely to be essential in discovering the underlying control functions, but their inaccessibility in human and animal experiments placed limitations 42 in the derived control theories. Here lies the significance of modeling. Various motor control theories derived from either experimental observation or inspired by robotic control implementations have been tested and further investigated using models either of particular subsets in neuromusculoskeletal apparatus, or the whole systems with simplified linear or nonlinear dynamics in the individual components. But it has been realized that many discrepancies between the proposed theories and observed behaviors might be explained by the sophisticated nonlinear dynamics of each compartment at various levels of the motor hierarchy (Cheng and Loeb 2008). These properties may constrain some potential solutions or enable others. In order to reveal these possibilities, a mathematical model of sensorimotor system of human arm aiming at characterizing the realistic muscle, sensor and mechanical plant is constructed in this work. It plays critical roles in synthesizing various levels of analysis and in revealing the prescribed control laws. We also believe that the study of complex/realistic behaviors (e.g. SDN, multi-joint arm with redundant muscles) of a sensorimotor system should focus, at least initially, on relatively simple regimes (e.g. a simple posture task) in which the measurable or observable states and parameters of the task (kinematics variability) won’t obscure the consequences of the underlying principles. For example, a large part of the trajectory variability and endpoint errors of a moving limb might come from intersegmental dynamics and/or the deliberately changed motor commands in the planning phase of the movements. In such a case, the observed motor control strategy, either from the EMG activities of the mono- and multi- articular muscles (Gribble et al. 2003), or through measuring the end-point stiffness matrix (Franklin et al. 2007), tends to involve mixed messages: The phasic muscle co-contractions at the end of reaching movements are needed to brake fast moving limbs for stable and accurate ending, and the stiffness modulation are to counteract the environmental disturbances. Although the SDN plays important role in motor coordination as suggested by the TOPS model (Harris and Wolpert 1998), either of these two factors can disguise the effects of SDN on the postural variability, thus the strategy CNS uses to control it. We choose to investigate the noise effect on a relatively simple postural control task to avoid the confounding factors. But this ‘simplification’ must not sacrifice the ‘realism’ of human motor systems that might have placed, in the first place, substantial constraints on the CNS to decide its control strategies. Instead, simple 43 but realistic model is needed to augment our understanding on the neural strategies by enhancing their tractability. A simple posture regime has been used to investigate the neuromuscular control of SDN in a single- joint system with a pair of symmetric muscles and constant moment arms (Selen et al. 2005). The study was informative in revealing the competing effects of SDN and muscle impedance on controlling the final kinematics variability. But we are wary that the work has over-simplified the musculo-skeletal system thus omitted many inherent properties that might potentially constrain the control of SDN. For example, the non-monotonic moment arm variations with joint angles can either facilitate the stabilizing effect of the muscles [in the positive slope of moment arm curves (Bunderson et al. 2008)] or deteriorate the stability of the joint particularly at high co-contraction levels [in the negative slope of moment arm curves (Shadmehr and Arbib 1992)]. And, the co-contraction of bi-articular muscles in a multi-joint arm system might have different stabilizing effect on distal versus proximal joint (Lan 2002). In consideration of these, we decided to study multi-joint, multi-muscle postural control tasks. The realistic properties such as moment arm variation (matched to anatomical measurements) and muscle contractile properties (using VM model and physiological architectural parameters) were modeled to investigate their effects on selective recruitment of multiple muscles to control SDN. Also a two-joint system is more similar to daily motor tasks. Additionally, when deciding the tasks under study, we differentiated a “free” posture task when there was no external loading and a postural force production task that aimed at examining the effect of external force on control of postural accuracy in the face of SDN. Both tasks are closely akin to the daily activity but the control strategy might be different. Although seemingly simple, the postural control is rather sensitive to many parameters under control: we varied the equilibrium joint angles and external loading magnitudes and directions, all of which factor into the control strategies available to the CNS. Our choice of the motor task obeys a key principle in modeling complex behaviors suggested by (Holmes et al. 2006): “to develop and validate models in constrained (limiting) situations before attempting to explain everything” . Therefore, we can say that we used a well-posed problem as a “template” to reveal basic principles while avoiding confounding factors that are either inevitable or unaddressed in many previous experimental and computational works. 44 vii. SIGNIFICANCE Understanding normal sensorimotor coordination seems likely to find application in the diagnosis and treatment of clinical dysfunction. Loss of movement precision is a distinct feature of many pathological conditions. The modeling studies presented here provide insights into the mechanisms whereby noise in the system gives rise to observable kinematic errors. The models also reveal potential strategies whereby the nervous system could minimize those errors, strategies that can be inferred from experimental observations such as electromyography. By observing the kinematic performance, the activities of the muscles and their adaptation over time, it should be possible to understand and differentiate mechanisms of dysfunctions that appear to be clinically similar. Furthermore, it may be possible to develop mental and physical exercises that improve the ability of patients to cope with or correct for dysfunctions. For severe dysfunctions such as paralysis, it will be necessary for prosthetic interfaces and control systems to take over control of the musculoskeletal plant, so-called functional electrical stimulation (FES). In that case, the controllers are likely to benefit from an understanding of the principles underlying biological sensorimotor control, which can then be incorporated into the prosthetic system through biomimetic design where appropriate. 45 CHAPTER II: VIRTUAL MUSCLE MODEL i. ABSTRACT We have improved the stability and computational efficiency of a physiologically realistic, virtual muscle (VM 3.*) model (Cheng et al. 2000) by a simpler structure of lumped fiber types and a novel recruitment algorithm. In the new version (VM 4.0), the mathematical equations are reformulated into state-space representation and structured into a CMEX S-function in SIMULINK. A continuous recruitment scheme approximates the discrete recruitment of slow and fast motor units under physiological conditions. This makes it possible to predict force output during smooth recruitment and derecruitment without having to simulate explicitly a large number of independently recruited units. We removed the intermediate state variable, effective length (L eff ), which had been introduced to model the delayed length dependency of the activation-frequency relationship, but which had little effect and could introduce instability under physiological conditions of use. Both of these changes greatly reduce the number of state variables with little loss of accuracy compared to the original VM. The performance of VM 4.0 was validated by comparison with VM 3.1.5 for both single-muscle force production and a multi-joint task. The improved VM 4.0 model is more suitable for analysis of neural control of movements and for design of prosthetic systems to restore lost or impaired motor functions. VM 4.0 is available via the internet and includes options to use the original VM model, which remains useful for detailed simulations of single motor unit behavior. ii. INTRODUCTION The original Virtual Muscle TM software (VM 3.*) (Cheng et al. 2000) was based on a series of muscle characterization experiments designed specifically to develop an accurate set of mathematical equations representing the physiological processes of force production (Brown et al. 1999; Brown and Loeb 1999, 2000a, 2000b). The data from those and other experiments were used to fit coefficients to describe the details of recruitment and force generation in skeletal muscles with mixed fiber-types. VM 3.* includes a motor neuron pool with multiple motor units recruited discretely, in size-ranked order, and with physiological frequency modulation and 2nd order activation kinetics (Brown and Loeb 2000b). The 46 lumped passive forces and aggregate active forces of all motor units interact with the muscle mass and visco-elastic tendon and aponeurosis (Scott and Loeb 1995). The model differs from the usual Hill-type models in that it includes length dependency of the activation-frequency (Af) and force-velocity (FV) relationships as well as sag and yield behaviors that are fiber-type specific (Brown et al. 1999; Brown and Loeb 2000b). These processes have usually been ignored or treated independently in other muscle models. VM 3.* was built in MATLAB with inter-connected Simulink basic blocks, which is inefficient in terms of memory utilization and computational speed. In the new VM, the mathematical equations provided in (Cheng et al. 2000) were re-formulated in state-space representation by explicitly choosing a set of state variables. The state-space model was then implemented in a single CMEX S-function of SIMULINK. The S-function is a computational module of numerical integration, which accepts a set of inputs to the model, updates model states based on their derivative equations, and calculates a set of outputs efficiently. It is represented by a single, customized Simulink block whose inputs and outputs can be linked to other Simulink blocks representing the other parts of the mechanical and control system. Such model implementation is convenient for neurophysiologists to build their own models efficiently using VM blocks for each muscle. VM 3.* was intended as a generic modeling tool for public use in neuromuscular physiology and engineering applications. In this study, we evaluated its computational performance in realistic multi- muscle and multi-joint settings. These simulations revealed that smooth modulation of force output could only be obtained with a large number of discrete motor units, which resulted in excessively long computational time. There was also instability in the interaction of antagonistic muscles around a joint that was traced to the intermediate state variable, effective length (L eff ), which had been introduced to model the delayed length dependency of the activation-frequency relationship. These problems were fixed in VM 4.0 and validated against VM 3.1.5 as presented here. A new muscle recruitment algorithm called “Natural Continuous” approximates the discrete recruitment of individual units in the original VM 3.*. The new recruitment strategy lumps multiple motor units of the same muscle fiber type into a single motor unit that is frequency modulated (see below). This greatly reduces the number of motor units needed to generate a smooth recruitment response. The new 47 Natural Continuous recruitment algorithm weights the force output of the lumped motor units within a fiber type and between the fiber types to give rise to smooth force production over the entire range of activation levels. Preliminary results have been presented elsewhere (Song et al. 2007b). VM 4.0 includes both the old and the new recruitment schemes, making it suitable for analysis of neural control of movements (Alstermark et al. 2007; Chan and Moran 2006; Lan et al. 2005), and for design of prosthetics and FES systems to restore lost or impaired motor functions. VM 4.0 is available for downloading from http://ami.usc.edu/projects/ami/projects/bion/musculoskeletal/virtual_muscle.html iii. MATERIAL AND METHODS a) Musculotendon dynamics The musculotendon dynamics of the VM model is modeled as a second order mechanical system which includes the muscle mass driven by the difference of forces generated in contractile element (fascicle) and series element (tendon) (Fig II.1). Figure II.1: Mechanical structure of a single and lumped unit in the VM model. The contractile element (CE) operates in parallel with the passive element (PE), which consists of stretching (PE1) and compressing components (PE2), to represent the fascicles. The fascicle element is in series with muscle mass element and series-elastic element (SE), which represent the combined tendon and aponeurosis. F pe1 results from the well recognized non-linear spring K pe1 with a viscosity element ( η pe1 ) that resists quick stretch and compression in the passive muscle, while F pe2 results from a non-linear spring K pe2 resisting compression at the thick myofilaments during active contraction at short lengths, thus is proportional to the activity of the muscle; F se results from a non-linear spring K se with a low stiffness toe region. The F ce’ is the force produced by combined contractile and passive components in the fascicle or contractile element; the difference of F ce’ and F se operate on the muscle mass to drive the muscle contraction dynamics. 48 b) Recruitment types VM 4.0 allows the user to build four different models of a given muscle based on the recruitment strategy: Natural Discrete (Brown & Cheng) provides the original VM 3.* model with discrete motor units and effective length (L eff ) and produces an output consisting of nested Simulink blocks. Natural Discrete provides the VM 4.0 model (L eff omitted; single S-function output) with discrete motor units and size-ordered recruitment similar to VM 3.*. Natural Continuous provides the VM 4.0 model with a single motor unit of each fiber-type that is frequency modulated and weighted to approximate size-ordered physiological recruitment. Intramuscular FES provides the VM 4.0 model with a single motor unit of each fiber-type whose frequency of firing is specified by an input variable representing stimulus frequency. Force output is weighted equally among the motor units to approximate the random fiber-type recruitment reported for intramuscular electrical stimulation by (Singh et al. 2000). Activation, U (0-1) Firing Frequency , f env (f 0.5 ) U th 1 U th 2 U th 3 U th 4 U th 5 U th 6 =U r 1 f max f min Activation, U (0-1) Firing Frequency , f env (f 0.5 ) U th 1 U th 2 U th 3 U th 4 U th 5 U th 6 =U r 1 f max f min (a) Activation, U (0-1) Firing Frequency , f env (f 0.5 ) U th slow =0.001 U th fast =U r *F 1 pcsa U r 1 f max f min Activation, U (0-1) Firing Frequency , f env (f 0.5 ) U th slow =0.001 U th fast =U r *F 1 pcsa U r 1 f max f min (b) Figure II.2: Continuous and discrete recruitment algorithms in Virtual Muscle model. (a) Natural Discrete recruitment algorithm as applied to a muscle consisting of 3 simulated slow-twitch and 3 fast-twitch motor units, respectively. U th i is the recruitment threshold of i th motor unit; U r is the activation level at which all the motor units are recruited. Once a motor unit is recruited, the firing frequency of the unit will rise linearly with U between f min and f max . This recruitment scheme mimics biologic recruitment of motor neurons. (b) Natural Continuous recruitment algorithm as applied to two fiber types, slow and fast. U th i is the recruitment threshold of i th fiber type or unit; the first fiber type U th 1 =U th slow is fixed at 0.001. Similar to the discrete algorithm, once a fiber type is recruited, the firing frequency of the fiber will rise linearly with U between f min and f max . 49 One important feature of the Virtual Muscle model is the built-in recruitment scheme of multiple motor units for different fiber types. The Natural Discrete algorithm in VM 3.* was formulated to model the multiple motor neuron pool and natural orderly recruitment of motor units (Fig II.2a) in response to the neural activation input, U (Cheng et al. 2000) . Activation incrementally recruits the motor units based on the recruitment rank of the fiber type, and once the unit is recruited, the firing frequency input of this unit (f env ) start from a non-zero value (from f min = 0.5 f 0.5 to f max = 2 f 0.5 ) and is modulated linearly over a range of frequencies specific to the unit type. The default range is a function of f 0.5 , the frequency for which the steadily recruited motor unit produces 50% of its maximal tetanic force, F 0 , at its optimal fascicle length, L ce0. This neural excitation signal is modeled as two separate first-order dynamics (f env -f int -f eff ) to simulate the rise and fall time of calcium kinetics (Brown and Loeb 2000b). The level of effective activation of each fiber type results from a linear combination of multiple motor unit activations (Af) weighted by their respective fractional PCSA (physiological cross-sectional area). The differences between slow and fast fiber types are reflected in rise-fall times of excitation dynamics, sagging (S) or yielding (Y) properties, Af relation, and muscle force-length, force-velocity ( FL, FV) properties. The sudden recruitment of a motor unit at its initial firing rate causes a step increase of muscle force; the size of that step is determined by the fractional PCSA of that unit, which is a function of the number of motor units modeled in the pool. If the number of motor units is sufficiently large (e.g. ~100 as occurs in a typical muscle), the factional PCSA can be reduced to produce a relatively smooth recruitment of muscle force. In order to reduce computational time, however, VM models are often constructed with only a few motor units, so the step jump in force recruitment is unphysiologically large, particularly at low levels of activation. To solve this problem, we designed a ‘Natural Continuous’ algorithm (Fig II.2b) to match the average behavior of the Natural Discrete recruitment strategy. Instead of modeling each unit explicitly, the Natural Continuous algorithm lumps the multiple units according to the corresponding fiber types, thus requiring only one unit per fiber type. Each unit or fiber type becomes active at a threshold U th i that depends on the distribution of fractional PCSA (F pcsa ) among all the fiber types and their recruitment orders: ⎪ ⎩ ⎪ ⎨ ⎧ > ⋅ = = = ∑ − = 1 1 001 . 0 1 1 i F U U i U i k k pcsa r i th i th (2.1) 50 Once recruited, the lumped motor unit modulates its frequency according to the usual calcium dynamics as in the Natural Discrete algorithm. There are three important aspects of the new algorithm. Firstly, the total muscle activation and contraction dynamics depend on the proportions of slow and fast fibers. The continuous algorithm accounts for this effect by linearly combining Af and FL, FV properties of each fiber types (Eq 2.2) multiplied by a weighting factor, W i , representing proportions of the fiber type among the overall active muscle portion (Eq 2.3). Secondly, the discrete algorithm simulates force modulations in physiological muscles by sequentially adding or substracting fractional PCSA scaled forces of newly recruited or derecruited motor units. As the number of units increases, this additional process is equivalent to multiplication by the activation level (U eff ), which formulates the continuous version of neural modulation of muscle forces in the new Natural Continuous algorithm (Eq 2.2). Thirdly, the activation frequency (Af) relationship includes the dynamics of calcium activation and this effect must also be represented in the multiplication of activation. We thus introduced a first order dynamics to convert activation input U to effective activation U eff (Eq 2.4). The values of rising and falling time constants (T U ) were chosen to simulate the dynamic force responses of Natural Discrete system to sinusoidal activation inputs over the range of physiological neural modulations (0-10Hz). The normalized active force from contractile element (CE) with n fiber types (units) is calculated as: () [] ∑ = + ⋅ ⋅ ⋅ = n i pe i i i i eff ce F FV FL Af W U F 1 2 (2.2) Where n is the number of active muscle fiber types, and the W i is calculated based on the threshold of each fiber type or unit, U th i , and the effective activation, U eff ,: () i th eff i k k th eff i th eff i U U U U U U W ≥ ∀ − − = ∑ = , 1 (2.3) The amount of muscle actually recruited is specified by an intermediate muscle activation signal, effective activation, U eff , which incorporates first order dynamics simulating the rise-fall effect modeled in calcium dynamics: 51 U eff eff T U U U − = & , ⎩ ⎨ ⎧ < ≥ = eff eff U U U U U T (sec) 15 . 0 (sec) 03 . 0 (2.4) c) Contractile dynamics A realistic model of the human elbow and shoulder was constructed from six VM 3.* muscles as described in (Song et al. 2008a). As described under Results, 3.1. Stability, this model exhibited unphysiological instability that was traced to the intermediate state variable, effective length (L eff ), which was modeled in VM 3.* using a nonlinear 1 st order differential equation: () () () [] () i L i eff ce ce i eff ce Af T t L t L t L − − = 1 3 & (2.5) ( L eff described in the text and figures is an abbreviation of “effective length” and it is represented by eff ce L in equations, which denotes the normalized effective length of muscle contractile element). The construct of effective length (L eff ) was introduced in VM 3.* to model the delayed length dependency of the activation-frequency (Af) relationship (Brown et al. 1999). The activation process that determines which cross-bridges can contribute active force acts with a delay, so any mechanical effects of length on the activation process itself will be related to the length of the muscle at that earlier time, rather than the current length and velocity, which exert their effects on the cross-bridges themselves (Brown et al. 1999). This is a relatively subtle effect that can be revealed by appropriate experimental design under conditions of controlled activation and kinematics (Brown et al., 1999). To model this effect, (Cheng et al. 2000) replaced the length input to the Af relationship with a time-delayed version of length (effective length: L eff ) as seen in Eq 2.6: () () ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = = fast n a f S L f Af slow n a Yf V L f Af Af i f i f n f f i eff i i eff ce i eff i n f f i eff ce i eff ce i eff i i , exp 1 , , exp 1 , , ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + = 1 1 1 0 i eff ce f f i f L n n n (2.6) 52 The destabilizing effect is removed by omitting L eff , and the new Af relationship is then modeled as a function of current muscle fascicle length (Eq 2.7). The performances of these two versions of the muscle model are compared in Results, 3.1. Stability. () () ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = = fast n a f S L f Af slow n a Yf V L f Af Af i f i f n f f i eff i ce i eff i n f f i eff ce ce i eff i i , exp 1 , , exp 1 , , ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + = 1 1 1 0 ce f f i f L n n n (2.7) Table II.1 provides the complete equations and coefficients for human muscle fiber-types in the reformulated Virtual Muscle model. 53 Table II.1: Summary of model equations and best-fit constants Equations Slow-twitch fibers Fast-twitch fibers () ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = 1 exp ln T T r se T T se se k L L k c L F 8 . 27 T c 0047 . T k 964 . 0 T r L 8 . 27 T c 0047 . T k 964 . 0 T r L () ce r ce ce ce ce pe V k L L L k c V L F η + ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = 1 exp ln , 1 1 max 1 1 1 0 . 23 1 c 046 . 0 1 k 17 . 1 1 r L 01 0. η 0 . 23 1 c 046 . 0 1 k 17 . 1 1 r L () () [] {} , 1 exp 2 2 2 2 − − = r ce ce pe L L k c L F 0 2 ≤ pe F 02 . 0 2 − c 0 . 21 2 − k 70 . 0 2 r L 02 . 0 2 − c 0 . 21 2 − k 70 . 0 2 r L () ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = ρ β ω 1 ce ce L exp L FL 12 . 1 ω 30 . 2 β 62 . 1 ρ 75 . 0 ω 55 . 1 β 12 . 2 ρ 88 . 7 max − V 88 . 5 0 v c 0 1 v c 15 . 9 max − V 70 . 5 0 − v c 18 . 9 1 v c () () () [] () []() ⎩ ⎨ ⎧ > + + + − ≤ + + − = 0 , 0 , , 2 2 1 0 1 0 max max ce ce v ce ce v ce v v v ce ce ce v v ce ce ce V V b V L a L a a b V V L c c V V V L V FV 70 . 4 0 − v a 41 . 8 1 v a 34 . 5 2 − v a 35 . 0 v b 53 . 1 0 − v a 0 1 v a 0 2 v a 69 . 0 v b () () ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = fast n a f S L f Af slow n a Yf V L f Af i f i f n f f i eff i ce i eff i n f f i eff ce ce i eff i , exp 1 , , exp 1 , , ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + = 1 1 1 0 ce f f i f L n n n 56 . 0 f a 1 . 2 0 f n 5 1 f n 56 . 0 f a 1 . 2 0 f n 3 . 3 1 f n 54 Table II.1, Continued () () Y Y ce Y T t Y V V exp c t Y − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − − = 1 1 & 35 . 0 Y c 1 . 0 Y V 200 ) (ms T Y – – – () () S i S i eff i T t S a f , t S − = & , ( ) () ⎪ ⎩ ⎪ ⎨ ⎧ ≥ < = 1 0 1 0 2 1 . t f , a . t f , a a i eff S i eff S S – – – 76 . 1 1 S a 96 . 0 2 S a 43 ) (ms T S () () ( ) i f i int i env ce i env i int T t f t f L , f , t f − = & () () () i f i eff i int ce i int i eff T t f t f L , f , t f − = & ( ) () ⎪ ⎩ ⎪ ⎨ ⎧ < + ≥ + = 0 , 0 , 4 3 2 2 1 i eff ce i f f i eff i env f ce f i f f L Af T T f t f T L T T & & 3 . 34 ) ( 1 ms T f 7 . 22 ) ( 2 ms T f 0 . 47 ) ( 3 ms T f 2 . 25 ) ( 4 ms T f 6 . 20 ) ( 1 ms T f 6 . 13 ) ( 2 ms T f 2 . 28 ) ( 3 ms T f 1 . 15 ) ( 4 ms T f U eff eff T U U U − = & , ⎩ ⎨ ⎧ < ≥ = eff U eff U U U U T U U T T 2 1 30 ) ( 1 ms T U 150 ) ( 2 ms T U 30 ) ( 1 ms T U 150 ) ( 2 ms T U Notes: Top bar⎯x denotes the normalized variable x (forces by maximum isometric tetanic muscle force F 0 , lengths and velocities by optimal fascicle length or optimal tendon length [L ce0 or L se0 ]); superscript x i denotes the i th motor unit specific variable x. 55 Table II.2: Symbols and definition of VM parameters Symbols Definitions U Activation input (0-1) i pcsa F Fractional PCSA of i th motor unit (0-1) i th U Recruitment threshold for i th motor unit (0-1) r U Fractional activation level at which all motor units for a given muscle are recruited (0- 1) eff U Effective activation, an intermediate muscle activation signal simulating calcium dynamics (0-1) i W Proportion of active i th fiber type/motor unit of total active muscle (0-1) 0 F Muscle maximal tetanic force (N) 0 ce L Optimal fascicle length (cm) 0 se L Optimal tendon length (cm) se F Series elastic element (tendon) force (F 0 ) 1 pe F Stretching contractile passive element (fascicle) force (F 0 ) 2 pe F Compressive contractile passive element (fascicle) force (F 0 ) ce F Active contractile element force (F 0 ) ' ce F Total contractile element force (F 0 ) se L Tendon length (L se0 ) ce L Fascicle length (L ce0 ) max ce L Maximum fascicle length of the muscle at its maximum anatomic musculotendon length (L ce0 ) ce V Fascicle velocity (L ce0 /s) i FL Force-length function of i th muscle fiber type i FV Force-velocity function of i th muscle fiber type i Af Activation-frequency relationship of i th motor unit Y Yielding factor for slow motor units i S Sagging factor for i th fast motor unit i env f Firing frequency input to 2 nd order excitation dynamics of i th motor unit (f 0.5 ) i f int Intermediate firing frequency of 2 nd order excitation dynamics of i th motor unit (f 0.5 ) i eff f Effective firing frequency of i th motor unit (f 0.5 ) i eff ce L Effective fascicle length(L ce0 ) 56 d) Validation To test the effectiveness of the Natural Continuous strategy, we compared its performance with that of the Natural Discrete VM model for both static and dynamic conditions of activation and kinematics. A human muscle model (deltoid posterior) with motor neuron pool composed of 60% slow and 40% fast fiber types was used. The muscle mass was 97.5 g (i.e.F 0 =266 N), optimal fascicle length was 11 cm, optimal tendon length was 2 cm, and maximal musculotendon length was 18.4 cm. Static force-activation relationships (F-U curves) of Natural Discrete muscle models with a small number of motor units (3 slow and 3 fast) and with a more physiological motor pool (30 slow and 30 fast) were compared with those of the Natural Continuous muscle model at constant muscle length input (L mt =16cm) over the full range of activation levels (U=0-1). The Natural Discrete algorithm models multiple units which differ in the distributions of firing rates that they produce in the various muscle fibers at different levels of activation. Frequency of firing has complex effects on the dynamics of the contractile apparatus and its FL, FV properties that are also fiber- type specific. Therefore, we compared the predicted force for the Natural Discrete model with 30 units of each type to the Natural Continuous model during sinusoidal modulation of activation (U=0-1) at constant muscle length (L mt = 16cm) and during sinusoidal length input from anatomical minimum to maximum (L mt =14 – 18 cm corresponding to range of joint motion) at constant activation level (U=0.5). Because the activation dynamics include delay terms, we systematically varied the frequency of sinusoidal modulation of activation (U) over the range 1-10Hz. iv. RESULTS a) Stability VM 3.* (“Natural Discrete (Brown & Cheng)” in VM 4.0) introduced an intermediate state variable, effective length (L eff ), to model the effect of delayed length dependency on activation-frequency relationships (Af) (Eq 2.6). This forms an internal negative feedback loop between muscle activation and musculo-tendon length with an Af dependent delay, which may produce instability of the musculo-skeletal 57 dynamics even without reflex feedback (Fig II.3a). Removing this variable (Eq 2.7) resulted in a stable open-loop response as shown in Fig II.3b. Computation time was also significantly reduced after removing the state variable of effective length (see Fig II.7 below). 0 2 4 6 8 10 40 50 60 70 80 90 100 Time (s) Joint angle (deg) el-flex sh-flex (a) 0 2 4 6 8 10 40 50 60 70 80 90 100 Time (s) Joint angle (deg) el-flex sh-flex (a) 0 5 10 15 20 25 30 20 30 40 50 60 70 80 sh-flex el-flex Joint angle (deg) Time (s) (b) Figure II.3: Effects of L eff on musculoskeletal dynamics. Simulations in (a) and (b) are actuated by discretely recruited muscle models with and without L eff . (a) A demonstration of the instability caused by L eff in Virtual Muscle through an open-loop dynamic simulation using a 2-DOF-6-Muscle arm model. The muscle activation commands of the antagonist muscles across shoulder (sh), elbow (el) and both joints (bi): U sh =0.9, U el =0.4, U bi =0.0. (b) A demonstration of stable open-loop response at the shoulder and elbow joints after removal of the effective length state. b) Accuracy for static conditions The Natural Continuous and Natural Discrete recruitment schemes were compared firstly in isometric contraction of the deltoid posterior muscle models at constant muscle length, 16 cm. The F-U curves of discrete model with 3 and 30 units for each fiber type were superimposed with that of continuous muscle in Fig II.4. The continuous recruitment strategy demonstrated a similar profile of F-U relationship over the full range of activation (U=0-1 in Fig II.4a). However, for the discrete model with fewer motor units (3 slow and 3 fast) the finer scale of the curve (U= 0 – 0.3) revealed abrupt jumps in force (~15N) associated with stepwise recruitment of the slow units (Fig II.4b). This effect was mitigated as number of units increases. As expected, the Natural Continuous recruitment strategy provided essentially a smoothed version of the same output forces, closer to what would be expected for a real muscle consisting of a much larger number of discrete motor units. 58 0 0.1 0.2 0.3 0 5 10 15 20 25 30 35 40 45 u (0-1) 0 0.3 0.5 1 0 20 40 60 80 100 120 140 160 180 u (0-1) Isometric Force (N) Continuous Discrete 3-Unit Discrete 30-Unit (a) (b) 0 0.1 0.2 0.3 0 5 10 15 20 25 30 35 40 45 u (0-1) 0 0.3 0.5 1 0 20 40 60 80 100 120 140 160 180 u (0-1) Isometric Force (N) Continuous Discrete 3-Unit Discrete 30-Unit (a) (b) Figure II.4: Comparison of steady-state force-activation relationship of continuous and discrete recruitment schemes. Deltoid posterior was at L mt =16 cm between the Natural Continuous algorithm and Natural Discrete models with a small motor neuron pool (3 slow and 3 fast motor units) and a large motor neuron pool (30 slow and 30 fast motor units). All the models have 60% slow fibers and 40% fast fibers. (a) F-U curves spanning entire activation rage (U=0-1); (b) F-U curves zoomed in at the low activation levels (U=0-0.3). c) Accuracy for dynamic conditions To test the effectiveness of the Natural Continuous algorithm in modeling the dynamic response of muscle to neural modulation and muscle length changes, the force outputs of continuous muscle model were compared with those of Natural Discrete model under input conditions of sinusoidal U at constant L mt (Fig II.5a); and a sinusoidal L mt at constant U (Fig II.5b). The continuous model produced force transients approximating those of the discrete model in response to both dynamic muscle activation and muscle length changes. There was a slight phase lead for the continuous recruitment strategy during relatively rapid modulation of neural activation U (Fig II.5c). In neural modulation of movement dynamics, the important control variables of muscle output usually include the mean force and force modulation ranges (top panel of Fig II.5e). The ratios of continuous to discrete recruitments for these two variables were close to unity over the physiological range of modulation of activation U (1-10 Hz) (Fig II.5e). 59 0 1 2 0 0.5 1 U (0-1) 0 1 2 16 0 0.5 1 1.5 2 0 50 100 150 200 Time (s) Force (N) 0 1 2 0.5 0 1 2 14 16 18 Lmt (cm) 0 1 2 0 50 100 150 200 Time (s) Lmt Lmt U U Natural Discrete Natural Continuous (a) (b) (d) (c) Mean Forces Modulation Ranges 0 2 4 6 8 10 0.5 1 1.5 2 Frequency of Sinusoidal Activation (Hz) Mean Force and Range Ratios Mean Force: Con/Dis Modulation Range: Con/Dis (e) Figure II.5: Comparison of dynamic responses between the Natural Continuous muscle model and Natural Discrete muscle model with large motor neuron pool (30 slow and 30 fast motor units). Both models have 60% slow fibers and 40% fast fibers. (c) shows the dynamic force responses given the muscle inputs depicted in (a), sinusoidal activation (U: 0-1 at 3 Hz) and constant length (L mt = 16 cm); (d) shows the dynamic force responses given the muscle inputs depicted in (b), sinusoidal muscle length (L mt : 14-18 at 3 Hz) and constant activation (U=0.5). (e) shows the ratios of mean forces and force modulation ranges between the continuous and discrete models over the physiological range of modulation of activation U (1-10 Hz). The effect of removing L eff on model accuracy was assessed by comparing the same set of sinusoidal modulations of activation and length in the two versions of the model (Eq 2.6 vs. 2.7). Fig II.6 illustrated that the discrepancy of force production between the muscle models with and without effective length is less than 1%. 0 1 2 0 20 40 60 80 100 120 140 Time (s) Muscle Force (N) with L eff no L eff 0 1 2 0 50 100 150 Time (s) (a)(b) 0 1 2 0 20 40 60 80 100 120 140 Time (s) Muscle Force (N) with L eff no L eff 0 1 2 0 50 100 150 Time (s) (a)(b) Figure II.6: Comparison of dynamic force responses between the Natural Discrete muscle models with and without effective lengths (L eff ). The models include a large motor neuron pool (30 slow and 30 fast motor units). Both models have 60% slow fibers and 40% fast fibers. (a) shows the dynamic force responses to the sinusoidal activation (U: 0-1 at 3 Hz) and constant length (L mt = 16 cm) inputs; (b) shows the dynamic force responses to sinusoidal muscle stretching (L mt : 14-18 at 3 Hz) and constant activation (U=0.5) inputs. 60 d) Simulation speed The time elapsed to compute force output for models with different recruitment strategies and different model implementations were compared in Fig II.7. In the comparison, two second simulations with a single muscle were carried out in Matlab 7.1/Simulink 6.3 on a HP workstation (3.2GHz Intel ® Pentium ® 4CPU, 2.0 GB of RAM). In the two discrete models, the computation time increased linearly with number of motor units. Removing the effective length substantially reduced the CPU time. The natural continuous model achieved a 15-fold reduction of computation time compared to the discrete model with 30 motor units and without L eff . 1 3 6 15 30 0 20 40 60 70 Number of Motor Units in Each Fiber Type Computation Time for 2 sec Simulation (s) Natural Continuous (No L eff ) Natural Discrete (No L eff ) Natural Discrete (Brown&Cheng) (with L eff ) 1.08 s Figure II.7: The comparison of computation time of different recruitment models and implementation. The original Natural Discrete model by Brown and Cheng includes L eff and is implemented by interconnected Simulink basic-blocks; the Natural Discrete and Natural Continuous models do not include L eff and are implemented using Simulink CMEX S- functions. v. DISCUSSION Muscles are organized into motor units, each consisting of a motor neuron plus a homogeneous group of muscle fibers that are activated synchronously by that motor neuron. Motor units are generally recruited in an orderly sequence based on both the fiber type and size (PCSA) of the unit (Cheng et al. 2000). The original VM model (Natural Discrete) was conceived to simulate this behavior of individual motor unit recruitment in response to a common synaptic excitation to all of the motor neurons in the pool. This is an important feature not present in most other muscle models and is useful for modeling details of 61 intramuscular mechanics such as the responses of Golgi tendon organs (Mileusnic and Loeb 2006). Ideally, for the Natural Discrete recruitment to be realistic, the number of motor units included in the model should approach the actual number of motor units in the muscle being modeled, which is often more than 100. Each motor unit adds 3 states to the model (after removing L eff ), substantially increasing the computational time (Fig II.7). The Natural Continuous recruitment scheme introduced in VM 4.0 improves the computational efficiency of the model without significant loss of its physiological realism and accuracy. Lumping motor units into one single structure has been done successfully in another influential muscle model (Zajac 1989), but the physiologically and mechanically different fiber types were not represented at all. Results of static force generation show that the continuous strategy approximates the natural recruitment behavior of motor units very well in the full range of motor neuron pool recruitment (Fig II.4). Dynamic simulation also shows that the natural continuous strategy is accurate for the dynamics of force output while requiring only one unit for each fiber type in the muscle (Fig II.5). The slight discrepancy in force output dynamics (Fig II.5 c,e) arises because the continuous algorithm replaces the distribution of firing frequencies across individual motor units with a single frequency of the same fiber-type that is modulated smoothly according to activation. Recently, there has been considerable interest in the variability of force output actually produced by muscles under steady-state conditions of kinematics and activation, i.e. motor noise. The amplitude of the noise tends to be a constant fraction of the mean force being generated, i.e. signal-dependent noise or constant coefficient of variance. This has been attributed to the incremental recruitment and derecruitment of motor units whose size is related to their recruitment order (Jones et al. 2002). Minimizing motor variability has been proposed as an organizational principle for movement control (Hamilton and Wolpert 2002; Harris and Wolpert 1998). In a simulation study to investigate how muscle stiffness may be controlled to minimize kinematic variability in the presence of motor noise (Selen et al. 2005), the authors proved that the lumped Hill-type muscle model with a single motor unit cannot replicate the realistic force variability. Instead, a motor neuron pool model with 60 motor units combined with Hill-type muscle dynamics must be used. In the VM model, the Natural Discrete recruitment algorithm is appropriate to 62 replicate the signal dependent feature of motor noise if the range of sizes of motor units is constructed appropriately and the activation signal has constant noise added to simulate membrane noise. Alternatively, the Natural Continuous model can be used with signal-dependant activation noise on its input, but the actual amplitude of the output force noise will depend on the low-pass filtering properties of the contractile apparatus. VM 4.0 also includes a model of recruitment by electrical stimulation (Intramuscular FES) that is designed to simulate the recruitment produced during functional electrical stimulation using intramuscular electrodes (Singh et al., 2000). As the stimulation pulse strength (equivalent to the activation input) is increased, motor units tend to be recruited at random because threshold depends mostly on distance of the motor axon from the electrode, which is unrelated to unit type or size. Motor units that are activated fire one and only one action potential for each stimulation pulse. Electrodes placed on the muscle nerve or on the surface of skin may produce different patterns of motor unit recruitment. A general recruitment scheme for FES could not be formulated to account for all stimulation techniques. A discrete model similar to the Natural Discrete recruitment scheme could be constructed if sufficient characterization information were available regarding the expected recruitment order of the individual motor units. vi. CONCLUSION The second generation of virtual muscle software (VM 4.0) described in this paper substantially reduces the computational load, and improves the robustness and stability of force computation, particularly when modeling multiple muscles connected to a skeletal linkage. The discrete model of recruitment (Natural Discrete) is useful to study detailed effects of the recruitment of individual motor units, such as motor noise (Jones et al. 2002). The lumped model with the new Natural Continuous scheme is much more efficient for computation with little loss of accuracy. The ability to easily modify and accurately represent various types of fibers in mixed and/or pathological muscles should be useful for applications of the VM software to neurophysiology, rehabilitation (FES) and sports medicine. 63 CHAPTER III: SENSORIMOTOR SYSTEMS MODEL OF HUMAN ARM i. ABSTRACT An integrated, sensorimotor virtual arm (VA) model has been developed and validated for simulation studies of control of human arm movements. Realistic anatomical features of shoulder, elbow and forearm joints were captured with a graphic modeling environment, SIMM. The model included 15 musculotendon elements acting at the shoulder, elbow and forearm. Muscle actions on joints were evaluated by SIMM generated moment arms that were matched to experimentally measured profiles. The Virtual Muscle TM (VM) model contained appropriate admixture of slow and fast twitch fibers with realistic physiological properties for force production. A realistic spindle model was embedded in each VM with inputs of fascicle length, gamma static ( γ stat ) and dynamic ( γ dyn ) controls and outputs of primary (I a ) and secondary (II) afferents. A piecewise linear model of Golgi Tendon Organ (GTO) represented the ensemble sampling (I b ) of the total muscle force at the tendon. All model components were integrated into a Simulink block using a special software tool. The complete VA model was validated with open-loop simulation at discrete hand positions within the full range of α and γ drives to extrafusal and intrafusal muscle fibers. The model behaviors were consistent with a wide variety of physiological phenomena. Spindle afferents were effectively modulated by fusimotor drives and hand positions of the arm. The VA model reproduced experimentally observed characteristics of arm mechanical impedance. These simulations validated the VA model as a computational tool for studying arm movement control. The VA model is available to researchers at http://pt.usc.edu/cel. ii. INTRODUCTION Sensorimotor dynamics and musculoskeletal biomechanics present inevitable constraints on motor control strategies that are evolved in the brain. The effects of these peripheral constraints on central control could not be evaluated with models in many previous studies, which have generally been limited to the particular subsets of the neuro-musculoskeletal system. These previous models were developed such that they were most relevant to the experimental phenomena under study with simplifying assumptions in order to fit the experimental results. Recently, it was recognized that systems models that were based on the 64 neurophysiology and biomechanics of the sensorimotor system could expand the utility of modeling approach in parallel with experimental investigation (Alstermark et al. 2007; Lan et al. 2005). Hierarchical models of the neuron-musculoskeletal system were employed to investigate control strategies of arm postures and reaching movements (Lan 1997; Lan and Crago 1994; Loeb et al. 1999), as well as neural control of locomotion (Yakovenko et al. 2004). A simulation study with a systems α−γ model suggested a plausible role for proprioceptive reflexes in controlling joint equilibrium position by way of γ static fusimotor command (Lan et al. 2005). Such a systems model may also help explain the behaviors observed in deafferented patients (Ghez et al. 1990; Gordon et al. 1995), and provide insights into the roles of proprioceptive afferents in neural control of movement that was revealed with muscle percussion (Capaday and Cooke 1983; Cordo et al. 1993). Realistic models of neuro-sensorimotor system could also provide a computational tool for neuroscientists to understand plausible neural strategies that are usually inferred indirectly from behavioral, psychophysical and electrophysiological data. There have been many recent advances in the sophistication and completeness of individual model components (Cheng et al. 2000; Lin and Crago 2002; Mileusnic et al. 2006; Mileusnic and Loeb 2006; Zajac 1989), and in the tools necessary to embody specific model structures in the computational environment in which the models run (Blemker and Delp 2005; Delp and Loan 1995; Holzbaur et al. 2005; Loeb et al. 1999). All of these now make it feasible to develop a realistic multi-joint, multi-muscle virtual arm (VA) model for computational studies of human motor control and learning. The VA model described here has shoulder, elbow and forearm degrees of freedom (DOF) with 15 muscles. The spindle and GTO models are embedded in each muscle to provide afferent information about the state of muscle contraction. The VA model is able to calculate muscle stiffness analytically from the virtual muscle model. The virtual arm model integrates individual component models of skeleton, virtual muscle, spindle, and GTO into a systems model (Fig III.1). Each of the component models is validated in isolated physiological conditions, and the parameters of the component models are determined in their prior validation in the development. Only those that are related to operation in the integrated system need to be adjusted, e.g. muscle length- tension property. The validation of the integrated model is, thus, focused on whether the gross behavior of system inputs and outputs is consistent to sensorimotor physiology, when operating condition and muscle 65 activation, rather than model parameters, are changed. The VA model is validated in this study with open- loop dynamic simulations at different hand positions in space and with the full range of α and γ inputs to the muscles and spindles. This model could be further integrated with neural control elements of the central nervous system (CNS) for simulation studies in motor control and learning. It could also be modified to simulate abnormal behaviors of the human motor system under various pathological conditions, e.g. deafferentation, stroke and spinal cord injury. Preliminary results have been reported elsewhere (Lan et al. 2006; Song et al. 2006 (a)) iii. MATERIAL AND METHODS The VA model was constructed from a set of component models, i.e. a musculoskeletal model, the virtual muscle (VM), the spindle model, and a GTO model. A systems model of the complete virtual arm is shown in Fig III.1, in which these model components are integrated to give rise to the outputs of joint kinematics and proprioceptive afferents under descending muscle (α) and spindle (γ) inputs. In the following sections, the component models and their integration are described. α −Commands γ − Commands SIMM Model of Arm Virtual Muscle Spindle GTO GTO K m u I a γ stat I b L ce L mt F m Joint Kinematics Proprioceptive Afferents II γ dyn Figure III.1: Diagram for the sensorimotor systems model integration. The Virtual Muscle receives α-command (u) and produces force outputs to drive SIMM musculoskeletal model for joint kinematics. The joint dynamics in turn passes the changes in musculotendon length (L mt ) back to the Virtual Muscle models. GTO model receives muscle force output from Virtual Muscle blocks and produce I b afferent firings. Muscle Spindle model is driven by dynamics, static γ-commands ( γ stat , γ dyn ) and muscle fascicle lengths (L ce ) from Virtual Muscle blocks, and outputs primary (I a ) and secondary (II) firings. 66 a) Musculoskeletal model Skeletal Structure The musculoskeletal arm model was developed in a graphical modeling environment, SIMM (version 3.2), which is a widely used modeling software (Delp and Loan 1995). Part of the musculoskeletal structure of shoulder complex was based on the model developed by (Holzbaur et al. 2005). The rest of the bony structure of the VA was based on an elbow model in SIMM (Lan and Baker 2004; Lan and Murakata 2001). Fig III.2a shows the planar view of the right arm in SIMM. The skeletal system consists of the thorax complex, including the thorax, the sternum, the clavicle and the scapula. The arm contains the humerus, the radius, the ulna and the wrist/hand bones. Fig III.2b shows the 15 muscles across the shoulder, elbow and forearm joints: deltoid anterior (DA), deltoid posterior (DP), the clavicular portion of pectoralis major (PC), supraspinatus (SS), infraspinatus (IS), biceps long head (Blh), biceps short head (Bsh), triceps long head (Tlh), triceps lateral head (Tlt), triceps medial head (Tmd), brachialis (BS), brachioradialis (BR), pronator teres (PT), pronator quadratus (PQ), and supinator (SP). The thorax complex stands as a ground for the arm and as a reference of the clavicle and scapular bones, which provide the attaching sites for shoulder muscles. The model does not yet include muscles operating the wrist and finger joints, so these joints are fixed in the posture illustrated and the segments contribute only to the mass and inertial properties of the forearm, a common constraint of many experimental paradigms for studying reaching movements (Abend et al. 1982; Gordon et al. 1995; Hogan 1985). 67 BR BS Blh Bsh PC DA DP IS SS Tlh Tlt Tmd PT PQ SP Shoulder F/E Elbow F/E Forearm P/S BR BS Blh Bsh PC DA DP IS SS Tlh Tlt Tmd PT PQ SP Shoulder F/E Elbow F/E Forearm P/S Figure III.2: The virtual arm musculo-skeletal model in SIMM. (a). The planar view of right arm skeletal model demonstrating the segments (ribs, sternum, clavicle, scapular, humerus, radius, ulna and bones of the wrist and hand), the joint coordinate systems, the degrees of freedom (DOF), shoulder flexion/extension (F/E), elbow F/E and forearm pronation/supination (P/S). (b). The origin, insertion points and paths of the 15 major muscle elements across the three joints (shoulder, elbow and forearm) of skeletal model. The 15 muscle elements include: deltoid anterior (DA), deltoid posterior (DP), clavicle portion of pectoralis major (PC), supraspinatus (SS), infraspinatus (IS), biceps long (Blh), biceps short (Bsh), triceps long (Tlh), triceps lateral (Tlt), triceps medial (Tmd), brachialis (BS), brachioradialis (BR), pronator teres (PT), pronator quadratus (PQ), and supinator (SP). The thorax functions as a ground and was included as the reference of the clavicle and scapular, which provided the attachment sites for shoulder muscles. The 2 DOF and 6 muscles used in later simulation are bold labeled in (a) and (b) Joint Kinematics The coordinate systems for each segment as shown in Fig III.2a are defined according to the convention recommended by the International Society of Biomechanics (Wu et al. 2005). The local coordinate system of the thorax coincides with the world coordinate system. The positions of clavicle and scapula relative to thorax are fixed at 30 o based on the shoulder rhythms when humerus is elevated (abducted) at 90 o in the frontal plane (Holzbaur et al. 2005). The articulation between the scapula and 68 humerus (glenohumerus or GH joint) is modeled as a ball-and-socket joint, with the center of the rotation located at the origin of the humerus local coordinate system. In the present model, the motion of the GH joint is constrained to the horizontal flexion/extension, because of available measurements of shoulder muscle moment arm are obtained in this plane, and the majority of studies in movement control are constrained to planar motion at shoulder level. However, this constraint can be easily removed for more general purpose of motion simulation. The elbow (or humeroulnar) joint is modeled as a hinge joint with the rotation axis passing between the center of the trochlear sulcus of the ulna and the center of the capitulum of the humerus with a carrying angle at about 5 o . The forearm rotates about an axis that passes through the center of the proximal radius and the center of the styloid process of the distal ulna. The range of elbow flexion is from 0 o (fully extended) to 130 o (fully flexed), and forearm rotation from –10 o (fully supinated) to 180 o (fully pronated). The neutral position of the forearm (0 o of pronation) is defined to keep the hand in the elevation plane as the arm stays in the horizontal plane at shoulder level. Muscle Origin/Insertion Points and Musculotendon Paths Muscle insertion and origin points (I/O points) were defined with a single point contact to bone surface. Muscle I/O points were initially chosen according to their anatomical landmarks; via-points in the musculotendon path and wrapping surfaces attached to the underlying bone segments were defined to represent the anatomical constraints at joints and in musculotendon paths over the full range of joint motion. I/O, via points and wrapping surfaces were adjusted as necessary to match experimentally measured values of moment arms from the literature (Haugstvedt et al. 2001; Murray et al. 1995). Measurements of shoulder moment arms in cadaver specimens are available only for horizontal flexion (Kuechle et al. 1997). Because experimental moment arms in elbow and forearm muscles were obtained with the shoulder joint at its neutral position (i.e. 0 o of elevation), SIMM generated moment arms of these muscles were fitted to experimental data in the similar neutral position of shoulder joint, where the humerus is in parallel with the Y axis of the thorax. Muscle I/O points in local coordinates are listed in Table III.1. 69 Table III.1: Muscle origin (O) and insertion (I) points Bones x y z Muscle Abbreviation O – origin I – insertion (cm) (cm) (cm) Shoulder Deltoid Anterior DA Clavicle (O) -0.95 0.82 6.75 Humerus (I) 0.60 -11.38 0.68 Posterior DP Scapula (O) -5.57 0.12 -2.51 Humerus (I) 0.21 -7.60 1.05 Supraspinatus SS Scapula (O) -5.59 -0.03 -8.10 Humerus (I) -1.54 -0.11 1.56 Infraspinatus IS Scapula (O) -7.79 -3.59 -6.62 Humerus (I) -1.13 -1.33 1.31 Pectoralis major Clavicular PC Clavicle (O) 0.80 -0.22 3.29 Humerus (I) 1.22 -5.83 0.48 Elbow Biceps Long Blh Scapula (O) -3.12 -2.35 -1.31 Radius (I) 0.96 -3.68 0.38 Short Bsh Scapula (O) 1.10 -3.92 -2.79 Radius (I) 1.00 -3.68 0.37 Triceps Lateral Tlt Humerus (O) -0.44 -5.95 0.70 Ulna (I) -1.89 1.27 0.02 Long Tlh Scapula (O) -4.86 -4.68 -1.71 Ulna (I) -1.89 1.27 0.02 Medial Tmd Humerus (O) -0.84 -13.70 -0.91 Ulna (I) -1.89 1.27 0.02 Brachialis BS Humerus (O) 0.40 -16.96 0.04 Ulna (I) -0.85 -3.03 0.46 Brachioradialis BR Humerus (O) -0.41 -20.88 0.07 Radius (I) 0.49 -20.86 2.63 Major forearm Pronator teres PT Humerus (O) 0.99 -28.17 -2.81 Radius (I) -0.54 -11.44 2.92 Pronator quadratus PQ Ulna (O) -2.27 -21.17 1.61 Radius (I) 0.14 -20.23 3.41 Supinator SP Ulna (O) -2.96 -1.38 -0.37 Radius (I) -0.74 -4.38 0.99 70 Segmental and Inertial Parameters Inertial parameters used in the VA model are listed in Table III.2. The size of bones was measured from SIMM, and was close to the cadaver measurements in (Veeger et al. 1997). Thus, the segment masses and inertias measured by (Veeger et al. 1997) was used for this model. The center of mass of each segment was estimated as a percentage of the segmental length from the proximal end of the bones (Seireg and Arvikar 1989). Since forearm mass and inertia were measured as properties of lumped radius and ulna in (Veeger et al. 1997), their values were evenly distributed to radius and ulna segments in the model (Lemay and Crago 1996). Table III.2: Segment mass and inertia parameters Segment Inertias (kg-cm 2 ) Bone Segments Bone Length (cm) Segment Mass (kg) Segment Mass Center (cm) * It Il Humerus 30.0 1.79 13.08 132.08 16.69 Ulna 25.2 0.545 10.36 28.17 3.48 Radius 23.3 0.545 9.72 28.17 3.48 Hand 18.5 0.46 4.3 28.29 4.07 b) Virtual Muscle model parameters The virtual muscle (VM) model implemented in this system is a modified version of the original VM model of (Cheng et al. 2000). For details of modification, audience is referred to an online appendix provided at website http://pt.usc.edu/cel The VM requires a large set of morphometric and architectual parameters including muscle mass (M m ), optimal fascicle length (L ce0 ), maximum musculotendon length (L mt max ), optimal tendon length (L se0 ), and the fraction of fiber type distribution (Tables 3 and 4). The muscle mass determines the maximum tetanic isometric force (F 0 ) of the muscle: 0 0 ce m L M F ρ ε = (3.1) with zero pennation angle, constant muscle density ( ρ=1.06 g/ cm 3 ) and specific tension ( ε=31.8 N/ cm 2 ). The virtual muscle model uses L se0 (the tendon length at maximal tetanic isometric force) instead of 71 tendon slack length, or L ses , (Zajac 1989), since L ses is less well defined than L se0 and tends to be about 5% shorter than L se0 , which was estimated as (Cheng et al. 2000) : L se0 =1.05* L ses (3.2) Table III.3: Muscle architectural parameters Maximum musculo-tendon length Muscle Mass Peak force Optimal fascicle length Optimal tendon length Muscle Abbre -viation L mt max (cm) M m (g) F 0 (N) L ce0 (cm) L se0 (cm) Shoulder Deltoid Anterior DA 21.07 420.93 1147.99 11 10.5 Posterior DP 18.37 122.36 265.99 13.8 4.2 Supraspinatus SS 14.68 120.87 490.00 7.40 7.5 Infraspinatus IS 16.06 413.37 1203.99 10.3 5.6 Pectoralis major Clavicular PC 18.06 206.27 363.99 17 0.4 Elbow Biceps Long Blh 40.54 335.99 629.99 16 24.5 Short Bsh 38.13 311.03 434.00 21.5 15.5 Triceps Lateral Tlt 26.54 285.60 611.17 13.6 13 Long Tlh 38.85 452.20 763.49 17 22 Medial Tmd 18.69 193.20 629.99 9.2 9.56 Brachialis BS 16.66 341.27 994.27 10.3 6.6 Brachioradialis BR 30.57 177.33 265.93 20 11 Major forearm Pronator teres PT 15.92 91.47 560.00 4.9 11.51 Pronator quadratus PQ 5.22 10.2 76.50 4 1.2 Supinator SP 6.97 66.01 476.00 4.16 2.94 72 Table III.4: Muscle fiber type fractional distribution and number of motor units Fractional PCSA (%) Johnson et al. (1972) Number of Motor Units Muscle SS S F SS S F Shoulder Deltoid Anterior 0 0.6 0.4 0 3 3 Posterior 0 0.6 0.4 0 3 3 Supraspinatus 0 0.6 0.4 0 3 3 Infraspinatus 0 0.6 0.4 0 3 3 Pectoralis major Clavicular 0 0.4 0.6 0 2 4 Elbow Biceps Long 0 0.3 0.7 0 2 4 Short 0 0.5 0.5 0 2 3 Triceps Lateral 0 0.5 0.5 0 2 3 Long 0 0.3 0.7 0 2 3 Medial 0 0.5 0.5 0 2 3 Brachialis 0 0.5 0.5 0 2 3 Brachioradialis 0 0.4 0.6 0 2 3 Major forearm Pronator teres * 0 0.5 0.5 0 2 3 Pronator quadratus * 0 0.5 0.5 0 2 3 Supinator * 0 0.5 0.5 0 2 3 The steps used to determine muscle architectural parameters in VM model were given as follows: 1. L mt max for each muscle element was estimated from musculoskeletal model in SIMM. The measurements made within planar range of joint motion may under estimate the true values of L mt max . However, this parameter can be tuned to maximize the force production capacity of muscles in different applications. 2. Initially, L ce0 and L ses were obtained from literature (Holzbaur et al. 2005) for each muscle, and then L se0 was calculated by Eq 3.2. Once L mt max , L ce0 and L se0 were given, the operating range of the normalized muscle fascicle length ( ce L =L ce/ L ce0 ) operating range were then determined by the internal algorithm in VM (Cheng et al. 2000). Based on desired ce L operating range, the initially 73 chosen L ce0 and L se0 could be inversely determined through trial-and-error methodology and adjusted at the same time to lie in the tolerable range of empirical measured values. In our VA model development, since muscles operate mostly in the ascending part of the length-tension curve in the vicinity of L ce0 , the initially chosen L ce0 and L se0 were adjusted so that ce L is constrained by 0.45 ≤ ce L ≤1 in the full range of joint motion. 3. F 0 for each of the 15 muscles were based on literature (Holzbaur et al. 2005). Then M m was calculated using Eq 3.1. 4. The percentage of fiber distribution for slow- and fast-twitch types was obtained from literature (Johnson et al. 1973) for the shoulder and elbow muscles. For forearm muscles, no data on fiber type distribution are available in literature. Thus, an even distribution of fiber type was assumed between slow- and fast-twitch fiber types. The total number of motor units included in each muscle ranged from 5 to 7, in order to achieve relatively smooth force recruitment without a heavy computational burden for simulation. c) Implementation of proprioceptor models We implemented the muscle spindle model developed previously (Mileusnic et al. 2006) in a Simulink S-function block shown in Fig III.3a. The spindle model is composed of three intrafusal fiber types, i.e. bag1, bag2 and chain fibers. Each spindle model has three inputs: fascicle length (L ce ), static ( γ stat ) and dynamic ( γ dyn ) fusimotor drives, and two outputs representing activity in type I a and type II sensory neurons. 74 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 γ stat (0~1) fiber activation level (0~1) saturation freq: f bag2 =100 Hz f chain =200 Hz bag2 chain Figure III.3: S-function implementation of the spindle model. (a) Matlab/Simulink S-function implementation of each intrafusal fibers in the spindle model. The spindle model which consists of three intrafusal fiber models (bag1, bag2, chain); two afferent firing summation nodes (I a /II afferent firing models); and the partial occlusion effect in primary afferent firing. Each intrafusal fiber model responds to two inputs: fascicle length and the relevant fusimotor drive. The spindle model generates two outputs: I a and II afferent activity. (b) Frequency-to-activation conversion and saturation effect for bag2 and chain fibers. The S-function spindle model skips the frequency-to-activation conversion in the original model by directly scaling the gamma static commands ( γ stat ) for bag2 and chain fibers according to their different saturation frequencies. Receiving the same γ stat drive (0-1), bag2 saturate at 0.5 corresponding to its saturation frequency of 100 Hz, whereas chain fiber saturate at 1, corresponding to its saturation frequency of 200 Hz. In the original spindle model (Mileusnic et al. 2006), the two fusimotor inputs, γ stat and γ dyn, are represented by the firing frequency (Hz) of activation, which are converted to an activation level between 0-1 for each intrafusal fibers within the spindle model. In the current implementation, the conversion from frequency to activation is replaced by directly using activation levels from 0-1 for γ stat and γ dyn inputs to the spindle model (Fig III.3b). Because γ dyn modulates the sensitivity of the bag 1 fiber only, the activation level of bag 1 fiber is directly equal to γ dyn . However, bag 2 and chain fibers receive the same γ dyn input, but saturate at different firing frequencies (f bag2 =100Hz; f chain =200Hz), the scaling relation between γ stat drive and activation level saturates first at γ stat = 0.5 for the bag 2 fiber, then at γ stat = 1.0 for the chain fiber (Fig III.3b). (a) (b) 75 The spindle model adopted in this VA system was validated by a large set of experimental data in (Mileusnic et al. 2006). Our implementation using a Simulink S-function block makes the model numerically more efficient in computation and can be inserted into multi-muscle systems modeled in the Matlab/Simulink environment. The values of model parameters for the spindle were taken from those fitted by experimental data used in (Mileusnic et al. 2006). Simulation of the S-function spindle block with the set of inputs identical to those in (Mileusnic et al. 2006) verified the correct implementation of the spindle model. The outputs of the S-function spindle model gave the same patterns of primary and secondary firings as produced in (Mileusnic et al. 2006). 0 5 10 15 20 0 500 1000 1500 2000 2500 3000 Force (N) Ib discharge rate (Hz) GTO m odel 1 Ib f(u) 1 Force 0 5 10 15 20 0 500 1000 1500 2000 2500 3000 Force (N) Ib discharge rate (Hz) GTO m odel 1 Ib f(u) 1 Force Figure III.4: Simulink implementation of a GTO model representing piece-wise linear static relation between the total muscle force and I b afferent firing to account for the ensemble response of the GTO in the muscle. The higher slope at the lower levels of muscle forces is due to the recruitment of tendon organs, which makes the predominant contribution to the total response at low forces. We use a piece-wise linear static relation (Fig III.4) between the total muscle force and I b afferent firing to account for the ensemble response of the GTO in the muscle. The higher slope at the lower levels of muscle forces is due to the disproportionate influence of the early-recruited slow-twitch muscle fibers to GTO output (Crago et al. 1982; Mileusnic and Loeb 2006). d) Model integration in Matlab/Simulink Integration of individual model components, such as the VM, the spindle model, the GTO model and the skeletal SIMM model, was accomplished using a software tool – musculoskeletal modeling in simulink (MMS) (Davoodi et al. 2001). MMS is a C program that automatically converts a SIMM model of a 76 sensorimotor system into a Simulink block that embodies its mechanical dynamics. The dynamics pipeline module of SIMM calls for SD/FAST to generate the set of equations of motion for the skeletal model, which is in turn converted into a kinetics block in Simulink by the MMS. The virtual muscle blocks are then interfaced with the kinetics block to effect joint movements. The MMS connects the spindle and GTO models with muscle fascicle lengths and forces. The α and γ commands are scaled neural inputs (0-1) to muscles and spindles respectively. The outputs are joint kinematics ( q , q & , q & & ), muscle activations (u), muscle forces (F m ), musculo-tendon lengths (L mt ), and the ensemble activity levels of spindle primary (I a ) and secondary (II) GTO (I b ) proprioceptors. e) Dynamic simulation using the VA model The integrated sensorimotor virtual arm was validated in the full range of muscle ( α) and spindle ( γ) commands with open-loop simulations. In the validation study, the forearm DOF was constrained at zero degree, leaving only two DOFs for the shoulder horizontal F/E and elbow F/E. Muscles selected here included Deltoid Posterior (DP) and Pectorailis Major (Clavicle portion, PC) for mono-articular shoulder joint muscles; Brachialis (BS) and Triceps laterial (Tlt) for mono-articular elbow joint muscles; and Biceps short (Bsh) and Triceps long (Tlh) as bi-articular muscles. The choice of this “subset” muscles was based on that: (1) a large number of experiments were performed with the similar joint configuration (Burdet et al. 2001; Gordon et al. 1995; Perreault et al. 2002); and (2) two pairs of mono-articular muscles and one pair of bi-articular muscles were commonly adopted in previous simulation approaches investigating multi-joint arm reaching movements (Flash 1987; Hogan 1985; Lan 1997; Schweighofer et al. 1998). In the first set of simulations, three patterns of open-loop α commands were selected to achieve three equilibrium hand positions in the workspace (A, B and C in Fig III.8a). At each hand position, additional five patterns of increasing muscle activation were then tuned to maintain the same equilibrium hand position with higher stiffness. In these simulations, constant spindle fusimotor inputs ( γ stat = γ dyn =0.3) were used to drive the spindle. One of the purposes of developing such a comprehensive model for sensorimotor system is to investigate the roles of muscle proprioceptive feedback in multi-joint arm posture and movement control. 77 The second set of simulations was designed to demonstrate the modulating effects of fusimotor drive to spindle primary (I a ) and secondary (II) afferents. Therefore, simulations were carried out with six incremental levels of γ commands ( γ stat = γ dyn ) to muscle spindles at each of the three hand positions, respectively. At a fixed hand position, any change in spindle I a and II firings would be due to the modulating effects of increased fusimotor commands. These simple open-loop simulations were designed to distinguish the effects of different factors on the sensorimotor responses, thus validating the VA model without undue complexity in the interpretation of simulation results. iv. RESULTS a) Parameterizing the VA model The VA model is the integration of the four sub-components of sensorimotor system: the musculoskeletal model, the VM, and the muscle spindle and the GTO models. Accurate parameterization will ensure proper operation of the VA model in simulation. The choice of parameters is not tailored to subject specific information, but rather is based on measurement values available in literature. However, this will provide a template of parameters for further specification. The torque generating capacity of muscles within the range of joint motion is the most important feature of the musculotendon unit, which is given by their moment arm profiles in Fig III.5. For shoulder muscles, i.e. PC, DA, DP, SS and IS in Fig III.5a, the solid lines show the SIMM generated moment arm profiles that demonstrate a good fit to experimental profile (Kuechle et al. 1997) in a large range of motion. Discrepancy in the extreme regions of joint motion may be attributed both to measurement errors and limitation of SIMM to capture complex musculotendon geometry. However, the peak of moment arms of these shoulder muscles are matched well to the measured values, suggesting that the model is able to replicate muscle capacity of moment generation in the large range of shoulder movement. There is a paucity of information on moment arm data available for both heads of biceps and long head of triceps at the shoulder joint (Wood et al. 1989). Their moment arm profiles are adjusted to match one measured value in Fig III.5a. The single point experimental value of moment arm indicated by a star falls closely to the moment arm profile from SIMM (solid red line) of the biceps (short head) at the shoulder joint. For elbow 78 muscles in Fig III.5b, SIMM generated moment arm profiles (solid lines) fit very well with those of experimental measurements (dotted lines) (Murray et al. 1995) . Similarly good fit is obtained for all forearm muscles (Fig III.5c) based on the data in (Haugstvedt et al. 2001). It is noted that the biceps have a significant supination action at the forearm. Therefore, activation of the biceps will produce significant actions on the forearm, elbow and shoulder joints. Such coupled actions cross multiple joints will have an important implication on the functional use of the biceps in motor tasks. 0 50 100 0 5 Pectoralis Major Shoulder flexion MA (cm) 0 50 100 0 5 Deltoid Anterior 0 50 100 -4 -2 0 Deltoid Posterior 0 50 100 -4 -2 0 Infraspinitas 0 50 100 -4 -2 0 Shoulder flexion angles (deg) Supraspinitas 0 50 100 -4 -2 0 2 4 6 Biceps Triceps long 0 50 100 0 5 Pectoralis Major Shoulder flexion MA (cm) 0 50 100 0 5 Deltoid Anterior 0 50 100 -4 -2 0 Deltoid Posterior 0 50 100 -4 -2 0 Infraspinitas 0 50 100 -4 -2 0 Shoulder flexion angles (deg) Supraspinitas 0 50 100 -4 -2 0 2 4 6 Biceps Triceps long 20 40 60 80 100 120 0 5 10 Biceps 20 40 60 80 100 120 -3 -2 -1 Triceps 20 40 60 80 100 120 0 5 10 Brachialis Elbow flexion MA (cm) 20 40 60 80 100 120 0 5 10 Brachiradialis Elbow flexion angles (deg) 20 40 60 80 100 120 0 5 10 Biceps 20 40 60 80 100 120 -3 -2 -1 Triceps 20 40 60 80 100 120 0 5 10 Brachialis Elbow flexion MA (cm) 20 40 60 80 100 120 0 5 10 Brachiradialis Elbow flexion angles (deg) (a) (b) 0 100 -1.5 -1 -0.5 0 Biceps Pronation MA (cm) 0 100 -1.5 -1 -0.5 0 Supinator 0 100 0 0.5 1 1.5 Pronator Quadratus 0 100 0 0.5 1 1.5 Pronator Teres Pronation angles (deg) 0 100 -1.5 -1 -0.5 0 Biceps Pronation MA (cm) 0 100 -1.5 -1 -0.5 0 Supinator 0 100 0 0.5 1 1.5 Pronator Quadratus 0 100 0 0.5 1 1.5 Pronator Teres Pronation angles (deg) (c) Figure III.5: Moment arm matching of shoulder flexion and extension (F/E) (a), elbow F/E (b), and forearm pronation and supination (P/S) (c). The dotted lines are the experiment measurements, the solid lines are from SIMM model. The muscle paths are validated through matching the model moment arm profiles of shoulder, elbow and forearm joints with those obtained from experimental measurements. Since the shoulder F/E moment arms of two heads of biceps and the long heads of triceps muscles are only measured at one single joint configuration, the model moment arm profiles are compared to the single data points in (a). A large set of parameters for each component model was adopted from literature (Crago et al. 1982; Holzbaur et al. 2005; Johnson et al. 1973; Mileusnic et al. 2006; Seireg and Arvikar 1989; Veeger et al. 1997). Tables 2, 3 and 4 list the parameter values of limb segment, muscle architecture and fiber fractional distribution used in the VA model. Among the parameters, muscle architectural parameters must be tuned, 79 so that the muscles operate in the region, in which their normalized fascicle lengths ( ce L ) are constrained to 0.45 ≤ ce L ≤1. Fig III.6 confirms that the procedure we used to fine-tune muscle parameters (see Methods, Section b) does yield normalized fascicle length that falls within this pre-specified range. Thus, it ensures that the muscles are operating in the ascending part of their length-tension relationship with positive intrinsic stiffness. 0 0.40.5 1 1.5 SP PQ PT BR BS Tmd Tlh Tlt Bsh Blh PC IS SS DP DA Normalized Fascicle Length L ce (L ce0 ) Muscles Deltoid Anterior Deltoid Posterior Supraspinitas Infraspinitas Pectoralis Major (Clavicle) Biceps long Biceps short Triceps lateral Triceps long Triceps medial Brachialis Brachiradialis Pronator teres Pronator quadratus Supinator Figure III.6: Diagram showing the operating range of fascicle length (L ce ) predicted for each muscle in the virtual arm model. The left and right edges of each dark bar indicate the values of min ce L and max ce L corresponding to the minimum and maximum physiological lengths of each musculotendon elements respectively, and the dark bars illustrate the portions of the length – tension curve on which muscle develops active force. There has been a wide range of literature reports on muscle peak force ( 0 F ), optimal fascicle length ( 0 ce L ), and tendon slack length ( ses L ). We compared the values of our VM parameters with those of experimental measurements and previous modeling approaches here in Table III.5 on 0 F , Table III.6 on 0 F and Table III.7 on ses L . There is a large variability in the literature on each of the three sets of parameters due to the possible differences in experimental preparations and specimen size (Garner and Pandy 2003; Garner and Pandy 2001; Gonzalez et al. 1996; Holzbaur et al. 2005; Langenderfer et al. 2004; Murray et al. 2000; Scott and Loeb 1994; Stein et al. 2004). However, our model parameters fall generally within the range of physiological values, and the virtual arm model reproduces the realistic force generating capability of human arm muscles. 80 Table III.5: Comparison of muscle peak forces with literature values Peak force F 0 (N) Muscle Model Holzbaur (2005) Garner (2003) Garner (2001) Gonzalez (1996 ) Shoulder Deltoid Anterior 1147.99 1142.60 277.48 Posterior 265.99 259.90 567.15 Supraspinatus 490.00 487.80 687.99 687.84 Infraspinatus 1203.99 1210.80 1100.13 1099.61 Pectoralis major Clavicular 363.99 364.40 342.46 Elbow Biceps 1063.99 849.29 143-251 Long 629.99 624.30 461.76 Short 434.00 435.60 392.91 Triceps 2004.65 2332.92 279-1040 Lateral 630.00 624.30 1268.87 Long 798.00 798.50 629.21 Medial 629.99 624.30 619.67 Brachialis 994.00 987.30 853.76 853.90 183-588 Brachioradialis 266.00 261.30 101.56 101.58 93-200 Major forearm Pronator teres 560.00 566.20 592.31 592.80 Pronator quadratus 76.50 75.50 Supinator 476.00 476.00 186.36 186.38 81 Table III.6: Comparison of optimal fascicle lengths with literature values Optimal fascicle length L ce0 (cm) Muscle Model Holzbaur (2005) Garner (2003) Garner (2001) Langenderfer (2004) Murray (2000) Gonzalez (1996) Shoulder Deltoid 12.80 Anterior 11 9.76 14.68 10.12(0.3) Posterior 13.8 13.67 17.02 14.18(2.52) Supraspinatus 7.40 6.80 4.28 4.28 7.07(0.4) Infraspinatus 10.3 7.60 6.76 6.76 8.74(2.46) Pectoralis major 19.00 Clavicular 17 14.42 22.65 14.95(3.1) Elbow Biceps 14.22 14.3-15.3 Long 16 11.60 15.36 15.61(0.3) 12.8(3.2) Short 21.5 15.00 13.07 18.09(0.38) 14.5(3.2) Triceps 8.77 6.7-14.5 Lateral 13.6 11.40 6.17 10.28(2.44) 9.3(2.8) Long 17 13.40 15.24 17.62(1.05) 12.7(2.1) Medial 9.2 9.20 4.90 14.46(0.87) Brachialis 10.3 9.00 10.28 10.28 9.42(2.32) 9.0(1.6) 9.0-18.5 Brachioradialis 20 16.40 27.03 27.03 17.53(1.79) 17.7(3.0) 14.2-23.0 Major forearm Pronator teres 4.9 4.90 4.48 4.48 5.5(1.2) Pronator quadratus 4 2.83 Supinator 4.16 3.30 6.04 6.04 82 Table III.7: Comparison of tendon slack length with literature values Tendon slack length L ses (cm) Muscle Model Holzbaur (2005) Garner (2003) Garner (2001) Langenderfer (2004) Murray (2000) Gonzalez (1996) Shoulder Deltoid 5.38 Anterior 10 9.30 1.64 2.6(2.3) Posterior 4 3.80 5.93 4.0(0.8) Supraspinatus 7.14 4.00 13.03 13.03 2.95(0.65) Infraspinatus 5.33 3.10 5.58 5.58 5.08(0.13) Pectoralis major 6.35 Clavicular 0.38 0.30 0.45 2.28(0.68) Elbow Biceps 22.98 21.00 Long 23.33 27.20 22.93 18.28(1.38) 22.9(1.6) Short 14.76 19.20 22.98 15.75(0.85) 18.3(2.5) Triceps 19.05 18.38 Lateral 12.38 9.80 19.64 16.7(0.65) 18.7(1.8) Long 20.95 14.30 19.05 19.95(0.6) 21.7(2.9) Medial 9.1 9.10 12.19 17.8(1.13) Brachialis 6.29 5.40 1.75 1.75 3.35(0.45) 11.6(1.3) 4.52 Brachioradialis 10.48 13.30 6.04 6.04 10.35(0.5) 16.9(1.7) 12.60 Major forearm Pronator teres 10.96 9.8 11.58 11.58 12.0(1.6) Pronator quadratus 1.14 0.50 Supinator 2.8 2.8 2.48 2.48 83 b) Dynamic validation of the VA model The dynamic behavior of the VA model was evaluated with simple patterns of open-loop activation. Fig III.7 illustrates an example of dynamic responses with open-loop constant α and γ inputs ( stat γ = dyn γ =0.3). It is shown that both the shoulder and elbow joints are stabilized quickly after an initial oscillation, because of the inherent stability of the neuromuscular system. The hand converges to a position in workspace (i.e. position C in Fig III.8a). The responses of internal variables of the system, such as muscle force (F m ), musculotendon length (L mt ), fascicle length (L ce ), and sensory responses of I a , I b and II afferents, reveal that their dynamic and steady state values are all within physiological ranges. This indicates that the model dynamics are propagating in time throughout the system from muscle (α) and spindle (γ) inputs to the output of arm configuration within the biomechanical and physiological constraints. 20 30 40 50 x y Hand position (cm) (a) 0 20 40 60 Joint angle (deg) sh-flex el-flex (b) 0 0.5 1 α (0~1) (c) 0.6 0.8 1 L ce (L ce0 ) (d) 0 20 Fm (N) (f) 0 5 10 15 0 50 100 150 Ia (Hz) (g) time (s) PC DP Bsh Tlh BS Tlt 0 5 10 15 0 50 II (Hz) (h) time (s) 0 5 10 15 0 2 4 Ib (kHz) (i) time (s) 10 20 30 40 (e) Lmt (cm) Figure III.7: Dynamic responses of one simulation at hand position C: (a) hand kinematics in the Cartesian coordinate originated at shoulder, (b) shoulder and elbow joint kinematics, (c) α command of the six muscles, (d) L ce , and (e) L ce , (f) muscle force F m , (g) spindle I a and (h) II firings, and (i) GTO I b firings of 6 muscles. Sensory responses of spindle and GTO at steady state are verified under different sets of open-loop simulations, in which the hand of the arm is positioned to three positions (A, B, and C in Fig III.8a) in space. In these simulations, when fusimotor drives are kept constant ( γ stat = γ dyn =0.3), the I a and II firing 84 frequencies display an insignificant change (SD/mean<10%) for all increasing α commands in each hand position (data not shown here). This is due to trivial changes in muscle fascicle lengths with increasing muscle activations when the hand position is fixed. It also implies that the serial elastic component of tendon has an insignificant effect on spindle afferent encoding of joint angle. However, I a and II firings exhibit a strong dependency on hand positions in space and muscle fascicle lengths. This is depicted in Fig III.8b, in which I a and II firing rates (the left axis) and fascicle lengths (L ce ) (the right axis) are plotted against the three hand positions for the six muscles. Such monotonic correlation suggests a plausible mechanism for the CNS to encode hand positions in the workspace by spindle afferents (Scott and Loeb 1994; Stein et al. 2004). At the three hand positions, however, when fusimotor drives are increased from 0 ~ 1, spindle I a and II afferents display a strong modulation by the γ drives (Fig III.9). The monotonic increase in I a and II firings with increasing fusimotor drives at the three hand positions illustrates the effective modulation of spindle afferents by the γ commands in the model. Thus, γ fusiomtor drive has a significant effect on spindle afferent encoding of joint angle, which should be taken into consideration in the encoding and decoding of joint angle from spindle afferents. A B C 0 80 160 PC Ia, II (Hz) 0.5 0.75 1 Ia II 0 80 160 DP 0.5 0.75 1 L ce (L ce0 ) Lce 0 80 160 Bsh Ia, II (Hz) 0.5 0.75 1 0 80 160 Tlh 0.5 0.75 1 L ce (L ce0 ) A B C 0 80 160 BS Ia, II (Hz) Positions 0.5 0.75 1 A B C 0 80 160 Tlt Positions 0.5 0.75 1 L ce (L ce0 ) (a) (b) Figure III.8: (a) Three hand positions (A, B, C) in space achieved with three sets of open-loop activations. (b) With the constant fusimotor drives ( γ stat = γ dyn =0.3), the steady state values of I a and II firings (left axis, red) and muscle fascicle lengths (right axis, black) are shown at the three hand positions (A, B and C) in the workspace. Each plot is for one of the six muscles. There is a monotonic relation between afferent firings with the arm positions and muscle fascicle lengths. 85 0 50 100 150 200 Ia (Hz) Position A 0 0.2 0.4 0.6 0.8 1 0 50 100 II (Hz) γ (0~1) 0 50 100 150 200 Position B 0 0.2 0.4 0.6 0.8 1 0 50 100 γ (0~1) PC DP Bsh Tlh BS Tlt 0 50 100 150 200 Position C 0 0.2 0.4 0.6 0.8 1 0 50 100 γ (0~1) Figure III.9: The relation between equilibrium I a (1 st row) and II (2 nd row) firing with γ ( γ stat = γ dyn ) activations for the six muscles at the three workspace locations (A, B and C). The figure demonstrates that once the primary and secondary afferents were activated, the firing frequency shows a monotonic relationship with the fusimotor drives at all of the three hand positions. v. DISCUSSION We are taking an integrative approach to modeling a multi-muscle, multi-joint, sensorimotor system of the human arm. Realistic component models of muscles and proprioceptors that are developed and validated previously have been integrated into a larger systems model in the Matlab/Simulink environment. This environment facilitates addition of models of neural controllers that have been hypothesized to account for experimentally observed movements and behaviors. The modular approach eases further development, improvement, maintenance and exchange of model components, and updating of model parameters based on new experimental data. The ability of a systems level model to replicate physiological phenomena with anatomically realistic components and parameters effectively “closes the loop” on a modeling process that can appear to be arbitrary when viewed from the limited perspective of individual components and subsystems. Systems models are beginning to demonstrate its usefulness in the analysis of neural control of motor behaviors. The realistic features of the VA model developed and validated here further support its validity as a tool for computational studies of motor control in parallel with experimental approaches. 86 The realistic feature of the VA is greatly enhanced by the digitized bony structure of the upper arm (Holzbaur et al. 2005), which provides a realistic constraint to muscle attachments and actions on joints. With such anatomical constraint, we are able to match muscle moment arms at joints of span to available experimental data closely (Fig III.5). Wrapping surfaces and path constraints are useful to define musculotendon paths, so that they do not cut into the bone surface. The profiles of muscle MA can be matched closely to those of experimental measurements by adjusting muscle insertion points and musculotendon paths. Results in Fig III.5 may be the best match we could achieve with a line representation of the musculotendon path. More realistic representation of musculotendon path will await 3D modeling of muscle geometry (Blemker and Delp 2005). The choice of muscles included in this model is dictated by the availability of experimental data in moment arm and architectural parameters. The anatomical structure of the arm in this model represents an average size of human upper extremity, and a subject specific model may be developed by scaling the average model to match the size of a specific subject. Another major realistic feature is the incorporation of the virtual muscle TM (VM) model (Cheng et al. 2000) to the VA model. The VM model is capable of modeling three types of muscle fibers, i.e. slow twitch (or fatigue), fast twitch (or fatigue) and slow-fatigue-fast twitch. The relative portion of each fiber type in the whole muscle can be specified. A linear combination of recruitment and frequency modulation scheme that is built in the VM model mimics the physiologic order of recruitment according to the size principle of the spinal motoneurons. Each fiber type has a similar structure of activation and contraction dynamics, but with a different set of parameter values for different twitch speeds fitted from a series of experiment measurements (Brown et al. 1999; Brown and Loeb 1999, 2000a, 2000b). In addition, the phenomena of muscle yielding, sag, rise and fall times in activation that have not been considered in previous muscle models are included in the VM model. These physiologic properties of muscle fibers make the VM the most realistic muscle model currently available. With its built-in recruitment scheme in the VM, a lumped motoneuron pool at the spinal cord can be used to represent the “final common path” of motor outputs from the CNS to muscle actuators. 87 The third important feature of the VA model is the integration of proprioceptor models into the virtual muscle. Particularly, the physiologically realistic spindle model (Mileusnic et al. 2006) used here includes the delicate structures and properties of a mammalian muscle spindle. Three types of intrafusal fibers (bag 1 , bag 2 and chain) are modeled and parameters are fitted to a large set of experimental data. The intrafusal fibers are innervated by γ stat and γ dyn inputs, and their outputs are combined to give rise to the primary (I a ) and secondary (II) afferents. In the model, a range of continuous modulation of fusimotor activity (from 0 to 1) is used to control the static and dynamic sensitivities of the spindle to length and velocity changes. This would allow investigation of biologically plausible roles of fusimotor control in motor task performance. In human muscles, a large number of spindles may be recruited in order to increase the range of length sensitivity (Cordo et al. 2002). However, this lumped spindle model (Mileusnic et al. 2006) has been calibrated using ensemble firing patterns recorded in animals to represent the collective response of a set of spindles. The continuous static and dynamic fusimotor drives to the lumped spindle model may be equivalent to the outputs from γ motoneurons in the spinal cord. In this sense, fusimotor control of spindle model will allow further investigation in proprioceptive control of multi-joint posture and movement. The VA model also produced consistent results in spindle afferents with respect to γ fusimotor control (Fig III.9) and muscle fascicle length (or hand position) change (Fig III.8). Sensory responses (Fig III.9) suggest that γ stat and γ dyn are generally effective in modulating I a and II firing rates of the spindle model. Fig III.9 shows the monotonic relation between spindle I a , II firings and muscle fascicle lengths (or hand positions) in the six muscles at different hand positions. The fact that spindle afferents vary proportionally with changes in fascicle length lends a direct support to the hypothesis that end-point positions of a multi- joint limb can be decoded from proprioceptive information originating from different muscles (Scott and Loeb 1994; Stein et al. 2004). However, the secondary afferent is more susceptible to unloading at short muscle length than the primary afferent, as seen in position A for BS and Bsh muscles and in position C for DP and Tlh muscles (Fig III.8). Thus, primary afferent would be more reliable to indicate muscle fiber length than secondary afferent. The results also indicate that if primary afferents of multiple muscles are used to encode or decode joint (or hand) positions of the arm, the effect of γ efferent to spindle must be taken into account. Furthermore, we noticed that at a fixed hand position, muscle fibers are maintained at 88 roughly constant lengths even when α commands increased substantially for each muscle. This is also a result of stiff series elastic component of the tendon. It indicates that spindle afferents remain nearly unaffected by increasing α commands. In other words, tendon elasticity may be a negligible factor in encoding or decoding joint and hand positions from spindle afferents. In applications, the VA model may be further tuned to suit simulation of specific tasks. The operating range of muscle length-tension property chosen in this study may not reflect the best functional range for a specific task, since we are mainly concerned here with validating the input-output behaviors of the VA model. Restricting muscles to operate in the ascending limb of their length-tension curve (i.e. 0.45 ≤ ce L ≤1) may be stringent for motor task control. We have relaxed this condition in a postural control task (Song et al. 2007a). It was demonstrated that a larger range of postural angles can be reached with same descending commands and spinal reflex gains, if a wider range of length-tension curve in the muscles is specified. Thus, in a task simulation, tuning muscle length-tension curves to obtain an optimum functional setting would further enhance the utility of the VA model. vi. CONCLUSION The VA model developed in this study is able to replicate a wide range of sensorimotor responses with simple patterns of open-loop, α and γ motor commands. The VA model represents a complete sensorimotor system that translates neural commands to the actuation forces (muscle model), the sensory afferents (spindle and GTO models), the multi-joint dynamics (skeletal model), and the ultimate motor behaviors at the hand. This neuromotor process produces a full set of physiologic variables that may correspond to those that are measurable in experiments (e.g. joint kinematics, joint torque/stiffness and hand force/stiffness, etc.), and also those that are not accessible to invasive instruments (fascicle lengths, muscle force, fusimotor drives, etc). In future studies, we will use this VA model to address motor control issues that have not been completely elucidated with data from human behavioral experiments. The model may also be expanded, or modified, to represent various neuropathologic conditions, such as deafferentation, stroke and spinal cord injury, and to address the strategies of neuro-rehabilitation. 89 CHAPTER IV: CONTROL STRATEGIES FOR SIGNAL DEPENDENT MOTOR NOISE i. ABSTRACT Signal dependent noise (SDN) in biological motor systems has been thought of as an important constraint on motor planning to achieve accurate movements. In a single-joint system, muscle co- contraction has been shown to be effective in reducing kinematic variability in the presence of SDN, but this strategy is energetically inefficient. The extension of this strategy to a multi-joint system is not obvious because of the intersegmental dynamics and mix of mono- and multi-articular muscles. Selective coactivation of such “redundant” muscles has been shown to be an efficient means for stabilizing limbs against external perturbation because it can be used to increase the mechanical impedance at the end-point of the limb. SDN in cocontracted muscles will give rise to end-point fluctuations in the absence of external perturbations, an effect that seems likely to influence impedance control strategies in multiarticular systems. This paper explores how various patterns of muscle activation affect both the magnitude and direction of such perturbations and the energy consumption under various limb position and loading conditions. We used a 2-DOF, 6-muscle human arm model and Monte Carlo dynamic simulations to determine the effect of different activation patterns and SDN on end-point postural variability. Using numerical stability analysis, we demonstrated that within a large set of coordination patterns that produced the required posture and end-point force, only a small subset satisfied both postural stability and energetic efficiency. Increasing energy consumption could generally be used to reduce postural variability but the best patterns of muscle activation at each energy level tended to change quite drastically rather than simply scaling in magnitude. The best patterns were also highly dependent on posture and external load conditions. If the measure of variability was constrained to a single, task-specific axis, the best solution could be one in which the variability in the orthogonal axis was actually larger than for other solutions. However, no strong correlation was found between the patterns (shape, orientation) of end-point variability and equilibrium hand stiffness, a commonly used indicator of end-point impedance. We conclude that in the presence of 90 SDN, the CNS can utilize muscle redundancy to maximize task performance but that the set of optimal strategies does not follow simple rules and probably must be learned rather than computed by the CNS. ii. INTRODUCTION Stabilizing limb posture is important for all aspects of motor behaviors: the accurate positioning of the hand at the starting point is essential for subsequent ballistic movement (Gordon et al. 1995); and stable and accurate posture of proximal segments is also a necessary condition of any manipulation tasks. However typically, limbs have fewer degrees-of-freedom (DOFs) than the number of muscles actuating them; different muscle activation patterns can result in the same joint torques and equilibrium positions of the limbs, suggesting a wide range of coordination strategies available to the CNS. To resolve this “muscle redundancy”, limb stability has been a central criterion for postural control tasks (Bunderson et al. 2008). This is because coordination of multiple muscles can be used to adjust limb impedance (viscosity and stiffness) through controlling the viscoelastic property of individual muscles (Osu and Gomi 1999; Perreault et al. 2002); and the impedance properties stabilize a mechanical system against external perturbations. However, perturbations also exist within human motor systems: There is noise in the motor commands and muscle forces, and the level of noise scales with the control amplitude (signal dependent noise or SDN) (Jones et al. 2002). The noise propagates through musculo-skeletal systems and gives rise to end-point postural variability even without any environmental disturbances. We hypothesize in this study that selective recruitment of “redundant” muscles is not only used to stabilize the limbs against unpredictable disturbances in the environment, but it can also fine-tune postural accuracy in the face of SDN. Recently, signal-dependent noise has attracted considerable attention in motor control studies. In situations where muscle force operates against fixed loads such as limb inertia, the kinematic variability associated with high force levels requires for rapid movements accounts well for the well-known inverse correlation between speed and accuracy described phenomenologically by Fitts’ Law (Fitts 1954). In situations where external perturbations occur unpredictably, antagonist muscles may generate forces that operate against each other instead of a fixed load, generating effects on limb impedance that are independent of limb kinematics. The enhanced limb impedance can obviously reduce the limb excursions 91 due to external perturbations, but the effects of cocontraction on internal motor fluctuations are less straightforward. This is because, the increased muscle activity not only provides higher level impedance to damp the kinematic effects of motor noise, but it is also accompanied by larger fluctuations of generated force (i.e. SDN). These two factors compete and constrain the choices of control policy for motor stability. In single-joint system, co-contraction of opposing muscles has been shown to be able to surpass the increased force noise, and thus was effective to reduce the end-point kinematic variability due to SDN in both postural and movement tasks (Osu et al. 2004; Selen et al. 2005; Selen et al. 2006). But the extension of this strategy to multi-joint systems is not obvious because of the intersegmental dynamics and mix of mono- and multi-articular muscles. In a stochastic model of the motor system, the biomechanical limb is essentially a noise filter, and the effect of SDN on the end-point kinematic output is characterized by the dynamic property of this filter. Multi-joint mechanics results in complicated dynamics that depend on joint position (such as apparent inertia) and relative motion of multiple links (such as interaction torque). How the impedance could reduce the internally generated noise in such a system is still an open question. Furthermore, in a stochastic motor system, the limb impedance may be only vaguely defined because the reference trajectory, or the equilibrium position itself could be perturbed by internally generated noise. From these perspectives, we might postulate that the SDN inherent in neuromuscular systems indeed presents another degree of freedom in addition to the stability criteria derived in a deterministic system (Bunderson et al. 2008), thus further reducing the dimensionality of a neuromuscular redundant biomechanical system. For example, among the control patterns satisfying stability criteria, only a small subset can achieve stable motor behavior with SDN. In many postural tasks, saving energy is probably as important a consideration as motor stability and positional accuracy. Higher impedance through muscle co-contraction may be effective in reducing detrimental effects of SDN, but it is energetically inefficient because the opposing action of muscles performs no net mechanical work, yet consumes metabolic energy (Hogan 1991). Assuming the biological system does not squander energy needlessly, one would expect co-contraction to be kept to a minimum except when the task demands it. Behavioral study on single-joint movements (Osu et al. 2004) showed that without accuracy demands, the subjects accepted worse performance with lower co-contraction. This 92 choice of energetic-optimal solution is rather straightforward because lowest co-contraction corresponds to the largest levels of kinematic variability, a simple mono-tonic relationship observed in single-joint systems (Osu et al. 2004; Selen et al. 2005; Selen et al. 2006). However, in a system with multiple DOFs, searching for this energy optimum with consideration of postural accuracy is a nontrivial task, because the limbs are actuated by a mixed mono- and multi-articular muscles, each of which has complex moment arm variation and muscle length-tension and force-velocity properties (Murray et al. 1995). As a result, individual muscle contribution to limb stability varies drastically with limb posture and motion, and selective activation of those with a biomechanical advantage presumably can save energy consumption without sacrificing task performance. For example, in a modeling study of postural control, (Bunderson et al. 2008) suggested that selective recruitment of muscles with higher local stiffness due to moment arm properties might be advantageous in stabilizing the limb with low level co-contraction. When observing point-to-point reaching movements, (Franklin et al. 2007) found that the subjects selectively stiffened the limbs along the direction of environmental instability by adjusting relative co-activations of mono- and bi-articular muscles, a strategy that is much more economical than simply scaling up overall size of stiffness. We propose that such a strategy might also be effective for optimal control of SDN, but the appropriate coordination of multiple muscles for efficient control may need a learning process. In redundant motor systems, motor coordination that aims at both performance achievement and energy efficiency has been extensively explored in many optimization models that have been successful not only in predicting invariant motor behaviors [such as the TOPS-α model by (Nagata et al. 2002)], but also in providing principled explanation of many key observations of redundant, noisy motor systems. One of such key features is ‘structured variability’, that is in redundant motor systems, the distribution of motor variability is unsymmetrical: the elementary fluctuations that degrade task achievements are constrained, whereas those that do not interfere with the performance are uncorrected. The feature is universal in a wide range of motor actions (Balasubramaniam et al. 2000; Cole and Abbs 1986; Danna-dos-Santos et al. 2007; Domkin et al. 2002; Haggard et al. 1995; Li et al. 1998; Scholz et al. 2003; Scholz and Schoner 1999; Scholz et al. 2000; Tseng et al. 2002; Valero-Cuevas et al. submitted; Winter 1989; Wright 1990) and they have been extensively explored and mathematically defined using “uncontrolled manifold” method (Scholz 93 et al. 2003; Scholz and Schoner 1999; Scholz et al. 2000). Motor variability in redundant motor systems is accompanied by multiple solutions for a given task. This abundance is beneficial because it affords the nervous system both the ability to stabilize the important performance variables and flexibility to compensate the noises or unexpected changes in the environment (Gelfand and Latash 2002, 1998; Latash et al. 2007), but it also presents a difficulty in understanding how the choice is made. By penalizing the energy consumption cost, optimal feedback control theory can resolve the redundancy by prescribing a ‘minimal intervention’ principle: “make no effort to correct deviations away from the average behavior unless those deviations interfere with task performance” (Todorov and Jordan 2002). Following this principle, we postulate that the postural variability due to SDN may also be structured in a task-relevant manner through selective co-activation of redundant muscles for maximizing performance variable and at the same time achieving energetic minimization. One shortcoming of many optimal control models (Harris and Wolpert 1998; Todorov and Jordan 2002) was that they only accounted for control of the dynamic movement, and ignored the postural control that has to deal with the static forces such as gravitational or muscle elasticity. Several behavioral observations (Brown et al. 2003; Nishikawa et al. 1999) have suggested the existence of separate movement and posture controllers. Even in a recent optimal control model (Guigon et al. 2007) where the static control was considered separately, the author only incorporated a simple mapping between states (the body configuration) and forces (the gravitational force). Such a static controller could not directly contribute to the stability of the equilibrium. In human motor systems with inherent SDN and unpredictable external environments, even a static postural control has to deal with dynamic interactions. A controller that can exploit the muscle viscoelastic properties has to be put in to consideration. As reviewed above, there have been many attempts to develop a comprehensive control theory that could unify multiple features of human motor behaviors, but to the knowledge of this author, none of them has successfully captured the entire picture. The reasons are multiple: it can be the muscle model being used (e.g. simple muscle models may limit the realistic noise property simulated in the optimal feedback control model in (Todorov and Jordan 2002)); or the inevitable limitations of experiments (e.g. external force perturbations could mask the effect of SDN in modulating impedance in reaching movements 94 (Franklin et al. 2007)); or the simplified systems under study (e.g. single-joint model limited the understanding of motor noise effect on the systems with multiple DOFs (Selen et al. 2005)). In the present study, we aimed at capturing a larger subset of the full range of noisy, redundant human motor systems to explore how multiple features of a redundant, noisy motor system may shed light on computational theories of motor control. To avoid many experimental and computational constraints that may cause us to overlook some of these features, we focused on a relatively simple task – loaded postural control in a 2-dimensional plane in presence of SDN. The system was simulated using a previously developed realistic musculo- skeletal model with 2 DOFs and 6 muscles. We hypothesized that in the presence of SDN, the CNS can utilize the “redundant” muscles to maximize task performance, and the optimal strategy relies on appropriate coordination of these muscles. Through Monte Carlo simulation, we explored how various patterns of muscle activation affect both the magnitude and direction of the end-point positional variability and the energy consumption under various limb position and loading conditions. Using numerical stability analysis, we demonstrated that within a large set of coordination patterns that produced the required posture and end-point force, only a small subset with appropriate coordination pattern satisfied both postural stability and energetic efficiency. Selective activation of multiple muscles could achieve postural accuracy along a task-specified direction, but the complicated interactions among biomechanical constraints, muscle property, and force production make the decision of the coordination less straightforward. iii. MATERIAL AND METHODS a) Methods overview We used a realistic musculo-skeletal model of the human arm developed in the previous work (Song et al. 2008a) to study motor redundancy and motor noise. The model was two-dimensional (2-D), constrained kinematically on the horizontal plane at the shoulder level (Fig IV.1). It had two DOFs (shoulder and elbow F/E) and six virtual muscles with mechanical properties based on published anatomical and physiological data (Kuechle et al. 1997; Murray et al. 1995; Seireg and Arvikar 1989; Veeger et al. 1997).We performed noise simulation and linear spectrum analysis to verify that the muscle model is able to produce the realistic isometric force variability. A reduced model with only elbow-joint and a pair of antagonist muscles was 95 then used to examine the effect of co-contraction on joint kinematic variability in a single-joint system. This effect was then examined in the 2-joint, 6-muscle model. The task was to maintain the arm at a posture in a 2-D horizontal plane while pushing against a zero or sub-maximal external force. To examine the use of neuromuscular redundancy on control of SDN in postural tasks, we first obtained the redundant muscle control patterns (or control points) that produced the identical postural end-point force through Monte Carlo simulation. The dynamical simulation with SDN was then conducted to measure the kinematic variability in response to various control patterns. With the large set of control points for each task configuration, we used a method of multi-objective (error and energy) optimization analysis to identify efficient and accurate control patterns in the presence of SDN. The task-relevant control of postural accuracy was investigated by measuring the end-point error cost along multiple specified directions (desired axis of accuracy or DAA). Equilibrium joint and end-point stiffness were calculated for further analysis on the effect of impedance modulation. SIMM Model of Arm SIMM Model of Arm Virtual Muscle Virtual Muscle L mt F m Joint Kinematics Bsh (Biceps Short) Tlh (Triceps Long) Tlt (Triceps Lateral) BS (Brachialis) DP (Deltoid Posterior) PC (Pectoralis Major- clavicle portion) 2 DOFs: 1.Shouder horizontal F/E 2.Elbow F/E q 2 q 1 Noisy Control x y T2: [42 o 50 o 0N –] T1:[50 o 90 o 0N –] T3: [42 o 50 o 3N ‐90 o ] Noisy Control x y T2: [42 o 50 o 0N –] T1:[50 o 90 o 0N –] T3: [42 o 50 o 3N ‐90 o ] (a) (b) Figure IV.1: (a) Two-joint, 6-muscle virtual arm model. The model is kinematically constrained n the horizontal plane at shoulder level; the 2 DOFs includes shoulder horizontal flexion/extension, and elbow flexion extension, each of which is modeled as hinge joint. q 1 , q 2 are the generalized coordinates representing shoulder and elbow flexion angles respectively. The six muscles include two pairs of mono-articular antagonists and one pair of biarticular antagonists; the muscles are color-coded that is consistent with the muscle control patterns shown as Fm bar plots in Fig IV.7-IV.9. (b) A schematic figure showing the coordinate frame and the three task configurations. 96 b) Musculoskeletal model and posture tasks The arm model used in this study was a two-DOF, six-muscle Virtual Arm developed in the previous work (Song et al. 2008a) (Fig IV.1a). The two DOFs – shoulder horizontal flexion/extension and elbow flexion/extension were modeled as pin joints and linked to the shoulder girdle formed with by the clavicle and scapula. The equations of motion of the skeletal model were expressed in the generalized coordinate system, q=[q 1 q 2 ] T ∈ℜ 2 , where the subscript ‘1’ and ‘2’ represent shoulder and elbow joints respectively. Limb motion was described by the equations of motion: () ( ) () q G q q C τ q M q − − = − & & & , ) ( 1 msl (4.1) where M ∈ℜ 2x2 is a positive definite symmetric inertia matrix corresponding to the mechanical property of human arm (Seireg and Arvikar 1989; Veeger et al. 1997); C ∈ℜ 2 is a vector of centripetal and Coriolis forces; G represents the gravitational force, which is ignored in the horizontal plane; and τ msl =[τ 1 τ 2 ] T ∈ℜ 2 is the joint torque generated by the six muscles: () ( ) u q, F q R τ m msl = (4.2) where R(q) ∈ℜ 2x6 is the muscle moment arm matrix modeled to match data from anatomical experiments (Kuechle et al. 1997; Murray et al. 1995); and F m (q,u) or F m (L mt (q),u) ∈ℜ 6 is the 6-muscle force vector modeled as function of q that determines the musculo-tendon lengths (L mt ), and muscle control signals u=[u 1 u 2 u 3 u 4 u 5 u 6 ] T ∈ℜ 6 that represent the scaled muscle activation levels (0-1). Muscle dynamics, F m [L mt (q),u], is modeled using Virtual Muscle with a natural continuous recruitment method that simulates the orderly recruitment and frequency modulation of a wide range of motor unit sizes (Song et al. 2008b). The six muscles (see Fig IV.1a) include one pair of mono-articular muscles around shoulder joint – clavicle portion of pectoralis major (PC) and deltoid posterior (DP), one pair of bi-articular muscles crossing both joints – short head of biceps (Bsh), long head of triceps (Tlh), and one pair of mono-articular muscles around elbow joint – brachialis (BS) and lateral head of triceps (Tlt). The colors of the six muscle names shown in Fig IV.1a are consistent across all the plots in this paper. Muscle architectural parameters were chosen according to their physiological values reported in the literature (Garner and Pandy 2003; 97 Garner and Pandy 2001; Gonzalez et al. 1996; Holzbaur et al. 2005; Langenderfer et al. 2004; Murray et al. 2000). We only modeled the slow fiber type because it is known that human slow and fast fiber types have differential functionalities, and slow, fatigue resistant fiber type performs prolonged and delicate activities such as posture (Lieber 2002). Modeling one fiber type also saves computational time because it reduces the number of states at each step of integration. A given motor task (Ti, i=1,2, or 3) is to maintain arm posture (q 1, q 2) while producing zero or sub- maximal hand force (|Fh|) along a specific direction ( ∠Fh). Out of a large range of task configurations, we selected three for presentation and analysis (T1, T2, and T3 in Fig IV.1b). The arm position was controlled by various patterns of constant muscle activation signals (u) with Gaussian distributed signal-dependent noise (SDN). The standard deviation (SD) of this noise increased linearly with the mean value of u, with a slope of 0.2 (i.e. constant coefficient of variation, CV=0.2). c) Signal-dependent noise and force variability in Virtual Muscle To examine the effect of SDN on arm posture control, the first goal was to verify the ability of the natural continuous VM to produce realistic isometric force variability. Natural continuous VM is the VM model that uses a continuous algorithm that was designed to approximate the discrete recruitment of slow and fast motor units under physiological conditions (Song et al. 2008b). Unlike physiological discrete recruitment method where large numbers of motor units are modeled for each fiber type, the natural continuous VM models only one unit per fiber type. This greatly reduces the computation time and preserves the accurate prediction of muscle force production. But the model structure (single-unit) violates the muscle physiology underlying the origin of the signal-dependent features of muscle force variability. Motor noise (force variability) arises from random fluctuations in membrane potential that affect the recruitment of all-or-none action potentials in the motor units. This neural noise is amplified by the increasing force-generating capabilities of the larger motor units that are recruited progressively according to their size (Jones et al. 2002). This results in signal-dependent force variability whose magnitude tends to scale with the mean force being generated. This effect can be reproduced by adding constant noise to the activation signal, u, to the Virtual Muscle model with a pool of motor units whose range of sizes is 98 constructed properly and the forces are modulated through natural discrete recruitment schemes (Song et al. 2008b). The question is whether a VM with continuous recruitment scheme and only one slow motor unit is also suitable to simulate realistic force variability. We examined this behavior using a single muscle model (deltoid posterior or DP) over the full range of mean activation levels (u=0-1) and muscle path lengths (L mt =14-18 cm). Gaussian distributed white noise was introduced in muscle input (u ~ N (μ u , σ u 2 )) with SD proportional to its mean (σ u = 0.2μ u , i.e. CV=0.2). Τen seconds of simulations were used to quantify force output variability; the noise amplitude (SD) was plotted against the mean force level to examine the signal dependent feature of the motor output. The muscle force variability has lower frequency content than the input white noise due to muscle activation and contraction dynamics. The power of the noise was found to be mainly below 10 Hz frequency (Taylor et al. 2003). To verify the low-pass filtering property of the VM model, we used the linear spectrum analysis method to examine the frequency response of DP model during isometric force at 140 N and muscle length of 16 cm. The same SDN was added to a mean motor input of u=0.8, and linear frequency response of the muscle was plotted as magnitude and phase transfer functions at this operating point. The power spectrum of both activation and force noise would be presented. d) Co-contraction in a single-joint system for postural accuracy Muscle co-contraction has been shown to be effective in reducing single joint kinematic variability (Selen et al. 2005; Selen et al. 2006). We reproduced this control policy in a single joint (elbow) system with two opposing muscles (BS and Tlt). The task was to produce an isometric torque (0.5 Nm) while maintaining an 80 o joint posture. With SDN in motor inputs (CV=0.2), the dynamic simulation resulted in joint variability. We plotted the SD of the joint variability with respect to the muscle co-activation levels to examine the effect of co-contraction. e) Monte Carlo dynamic simulations We performed Monte Carlo simulations with deterministic muscle activation constrained by 0 ≤ u ≤1. The purpose was to solve the “ inverse ” problem (muscle redundancy) to identify the muscle coordination patterns (u) that produce the specified hand force. The simulation consisted of about one 99 million iterations, and the hand force accuracy was achieved by constraining the solution that produced the force vector within 1 o of directional accuracy and 0.1 N of magnitude accuracy (negligible deviation compatible with the experimental conditions). In the forward dynamic simulations, deterministic control u was first used to stabilize the arm at an equilibrium position. Muscle stiffness matrix (K m ) was calculated analytically from Virtual Muscle for each control pattern and the corresponding joint and end-point stiffness matrices (K j , K h ) were obtained through stiffness transformation (Hogan, 1991). The hand stiffness field was represented geometrically by a stiffness ellipse (Mussa-Ivaldi et al. 1985), whose major and minor axes were determined by the eigenvectors and eigenvalues of K h . The major axis represents the direction along which the arm is stiffest, and the minor axis is the direction of minimal stiffness. To facilitate interpretation, hand stiffness ellipses were characterized by three parameters: size, shape and orientation. The size corresponds to the area of the ellipse proportional to the determinant of K h matrix; the shape is the ratio of the major to minor axis length or larger to smaller eigenvalues (M/m) and represents the eccentricity of hand stiffness field; the orientation is defined as the counter-clockwise angle of the major axis direction (the principal eigenvector) with respect to x axis. To examine the control of SDN for postural accuracy, Gaussian random noise with 20% CV was introduced in each u vector once the system stabilized at the equilibrium position. Fifteen-second stochastic simulation resulted in random hand fluctuations. The spatial distribution of the end-point variability was characterized by an ellipses-like contour following the method defined in (Gordon et al. 1995). The orientation of the contour was determined by the principal axis of the hand distribution following the principle components procedure; the size and shape was based on the interquatile range in major and minor axes. The process follows: 1) the major axis was determined through principle component method, and the minor axis is perpendicular to it; 2) the median of the distribution along the two orthogonal axis was defined as the center of the contour; 3) the distance of each line from the center was calculated as the distance of the farthest point in each direction along that axis that did not exceed 1.5 times the interquatile range (a procedure used to draw ‘box plots’ for one-dimensional data); 4) finally the adjacent ends of the axes were connected with curved lines each of which was equivalent to a quarter of ellipses. If the hand 100 fluctuation is distributed symmetrically along two axes, the contour is equivalent to an ellipse. The major axis gives the direction with largest variability. Similar to stiffness ellipses, the hand random distribution is also described by the three parameters – size, shape and orientation. f) Postural variability versus energy consumption Like most problems in nature, human motor tasks usually have multiple (possibly conflicting) objectives to be satisfied. In posture control tasks, improved accuracy and stability is often at the cost of high energy consumption. These conflicting criteria – error and energy – form a Multiobjective Optimization Problem (MOP), where the term “optimize” essentially means finding a solution that generates values of all the objective functions acceptable to the decision maker (Osyczka 1985). The notion of optimum in multiobjective optimization is termed as “Pareto optimum”. In such a problem, there is almost always not a single optimal solution, but rather a set of solutions called the Pareto optimal set. The plot of the objective functions whose decision variables (the muscle control patterns u) are in the Pareto optimal set is called the Pareto front and we call the points in the Pareto front as Pareto frontiers (PF). In the current arm postural control study, we explored the high dimensional control space (decision variable, u) through Monte Carlo dynamic simulations and measured the two corresponding performance objectives – positional error and energy. The positional error is defined as the size of the hand random distribution in the Cartesian space (err=Size Hxy ), and the energy cost by a sum of squared muscle forces (E= ∑ i F mi 2 ). Here we modeled energy consumption that increases supra-linearly with muscle force, because muscle is less efficient at higher force levels where fast, fatigable motor units need more energy for higher rate of calcium turn over. Through a nondominated sorting procedure, the Pareto front was identified from the set of points in the multi-objective space (err-E space) that consisted of the desired control patterns (u * ) for efficient control of postural accuracy with the presence of SDN. Note here the Pareto optimal set was identified from the limited number of Monte Carlo simulations thus the “optimum” mentioned in this work only represents the best control patterns among the discrete Monte Carlo set, and may have missed the true optimum. However, the Monte Carlo method can avoid local minima by exploring the entire control space and the Pareto frontiers can be good initial conditions for subsequent continuous optimization process. 101 The Pareto front also reflects the error-energy trade-off, where the low energy point is often associated with high variability. This is consistent with the error-coactivation relationship identified in the single-joint system (Selen et al. 2005; Selen et al. 2006), but to achieve higher accuracy, a simple co-contraction may not apply in a multi-joint system, particularly with redundant muscles. We tested this idea by proportionally increasing the control pattern, i.e. the muscle force vector of the low energy point in the Pareto front, and the resulting kinematic variability was compared to the minimum variance point in the best set. Further more, the hand random drifting of a two-joint arm has spatial patterns that may be controlled through coordinating multiple muscles and modulating limb impedance. To examine this point, we presented the control patterns as bars of muscle force in the Pareto front with the corresponding Kh ellipses and hand drifting in the end-effector space. g) Control of postural variability along desired axis of accuracy We studied the task-relevant control of arm posture by sorting the best control patterns along the desired axis of accuracy (DAA). In this case, the error cost in the multi-objective space was defined as the standard deviation of hand random drifting along a specified DAA (err=SD DAA ). We defined two sets of DAA for this measurement. The first set was determined by the task requirement, and is independent of the joint configuration. Two orthogonal axes among the evenly distributed four axes in the Cartesian space (a1- a4 in Fig IV.2) were selected so that they were along and perpendicular to the end-point force vectors. In the case of zero hand force (Fh=0N), a1 and a3 that are along the x, y world coordinate frame were used. The second set of DAA was defined by the major and minor axes of the arm manipulability ellipsoid (m1, m2 in Fig IV.2). The purpose was to examine the flexible posture control constrained by the joint geometry. The manipulability ellipsoid was formulated by the limb manipulability matrix M(q) that is defined by arm Jacobian: M(q)=J(q)J T (q) (Sabes and Jordan 1997). It measures the sensitivity of the end- point to the uncertain control in the joint space, i.e. same amount of joint noise will result in largest end- point variability along the major axis of M(q) ellipsoid. Thus the major axis indicates the most difficult direction for controlling movement accuracy. 102 m2 m1 a1 a2 a3 a4 Figure IV.2: A schematic diagram showing the two sets of desired axis of accuracy (DAA): (1) four axes (a1-a4) are evenly distributed on the extra-personal space, separated by 45 o ; (2) two orthogonal axes are defined in the intro- personal space according to the arm manipulability ellipsoid: m1 is the major axis and m2 is the minor axis of the ellipsoid. Using MOP analysis, we identified the Pareto front for each of the two orthogonal axes in the two DAA sets. To reduce the noise effect along a particular axis, the CNS could choose the strategy that scales down the overall size of end-point variability (Size Hxy or |Hxy|). Alternatively, it could reshape or reorient the Hxy distribution to achieve accuracy along the task-specified direction, a strategy that might be more efficient to achieve postural accuracy. To test this point, we plotted the SD in the Pareto front along each DAA with respect to the overall size of hand variability (SD DAA vs |Hxy|). A monotonic increase of SD DAA with |Hxy| indicates a less flexible control along that axis or direction. Furthermore, the CNS might achieve such flexible controls through coordinating redundant muscles and modulating limb impedance. Muscle control patterns and Kh ellipses were presented for each control point in the Pareto front for this analysis. The shape and orientation of hand distributions (Hxy) were plotted against shape and orientation of Kh ellipses to examine the correlation between the two. iv. RESULTS a) Natural Continuous Virtual Muscle simulates realistic force variability Fig IV.3 shows how the SD of the muscle force of a natural continuous VM increases with increasing mean force levels. The results is plotted for the entire range of muscle operating lengths of deltoid posterior (L mt =14~18cm). The consistent slope (CV=2%) is in accordance with the measured values using physiological muscles (2-5%) (Jones et al. 2002; Taylor et al. 2003). 103 0 50 100 150 200 250 0 1 2 3 4 Mean Force (N) SD Force (N) CV=SD/mean~0.02 Figure IV.3: Force variability expressed as the SD as a function of mean force levels at full range of Deltoid Posterior (DP) muscle operating lengths (L mt =14~18cm). Given Gaussian distributed SDN with the noise amplitude proportional to the mean level of activation (CV u =0.2), the natural continuous recruitment scheme of VM can generate signal- dependent force variability whose SD is linearly scaled with force magnitude with a slope, CV=2%, within the range of experimental measurement. In addition to the magnitude property of force variability, we also examined the power spectrum of the motor noise. Fig IV.4 illustrates the transfer of the broad band white noise in muscle activation input (Fig IV.4a) to the noise in muscle force output (Fig IV.4c). Fig IV.4b shows the magnitude and phase response of the linear transfer function, T(f) = Fm(f)/u(f), at the operating point of u=0.8, L mt =16cm and Fm=140N. The plot demonstrates low-pass filtering property of the VM model (Fig IV.4b1). The corresponding amplitude spectrum of force noise (Fig IV.4c2) shows a concentrated power below 15 Hz, consistent with the experimental measurements (Taylor et al. 2003). 104 0 2 4 6 -1 0 1 Time Series of Noise in Muscle Activation (u) noise of u (0-1) 0 500 1000 1500 2000 2500 0 2 4 x 10 -3 Amplitude Spectrum of u Frequency (Hz) |u(f)| 0 2 4 6 -5 0 5 noise of Fm (N) Time Series of Noise in Muscle Force (Fm) 0 50 100 150 0 0.2 0.4 Amplitude Spectrum of Fm Frequency (Hz) |Fm(f)| 0 250 500 0 50 100 T(f)=Fm(f)/u(f) (magnitude) 0 250 500 -2000 0 T(f)=Fm(f)/u(f) (phase) Frequency (Hz) Degrees Virtual Muscle (Natural Continuous) Force (Fm N) Activation (u, 0-1) (a1) (a2) (b1) (b2) (c1) (c2) Figure IV.4: Linear spectrum analysis of the natural continuous Virtual Muscle model of deltoid posterior at operating point of u=0.8, L mt =16cm and Fm=140N. (a) Time series of broad-band white noise in muscle activation input u, and its amplitude spectrum. (b) Magnitude and phase frequency transfer function of the Natural Continuous Virtual Muscle. (c) Time series of noise in muscle force output, and its amplitude spectrum. b) Co-contraction can suppress SDN for single-joint postural accuracy Fig IV.5 shows the kinematics variability (SD) of elbow joint angles as a function of muscle co- activation. With increasing muscle co-contraction, the increased joint impedance surpasses the increasing levels of motor noise and led to a monotonic reduction of kinematics variability. This agrees with the general trend observed in previous simulation and experimental studies (Selen et al. 2005; Selen et al. 2006) 0 0.2 0.4 0.6 0.8 0 1 2 3 4 5 6 Co-activation level Kinematic variability (deg) Figure IV.5: Elbow joint kinematics variability (SD) as a function of muscle co-activation level at joint angle of 80 o , and isometric joint torque at 0.5Nm. 105 c) CNS coordinate redundant muscles for multi-joint postural accuracy In a set of muscle activation patterns that produced the same arm posture and end-point force, there were dynamical responses with varying levels of stability. We grouped them into three: non-stable, meta- stable, and stable (see Fig IV.6). The non-stable group was characterized by unstable responses even at the first 15 seconds of deterministic simulation. In this group, the arm oscillated in both joints with large amplitude and could not stabilize at the specified equilibrium posture. In this case, we did not proceed to stochastic simulation with SDN (Fig IV.6 (a) and (b)).The meta-stable group could stabilize in equilibrium position within 15 seconds, but once SDN was introduced, large oscillation resulted from dynamical interaction with the constant external forces applied at the end-point (Fig IV.6 (c) and (d)). The unstable behavior in this case masked the effect of SDN thus was excluded from the further analysis on postural control with SDN. The stable group demonstrated the stable behavior with both deterministic and noisy motor inputs (Fig IV.6 (e) and (f)); the joint kinematics variability during the 15-second stochastic simulation characterized the “true” filtering effect of the musculo-skeletal system, so this data group was used for quantifying end-point variability and the multi-objective optimization analysis. 0 20 40 60 50 100 40 42 44 40 50 60 0 10 20 30 41.5 42 42.5 0 10 20 30 49 50 51 SDN SDN Shoulder Elbow Non- stable: Meta- stable: Stable: Joint Angle (deg) Time (sec) Time (sec) (a) (b) (c) (d) (e) (f) Figure IV.6: Three types of dynamics responses of shoulder and elbow joint kinematics from the activation patterns that produces same equilibrium posture and static end-point force. The there types are classified by the level of their stability: (1) unstable – joints can not be stabilized during the first 15 second of deterministic simulation; (2) meta- stable – joints can be stabilized in the deterministic simulation, but once the SDN is introduced, non-stationary oscillation occurs; and (3) stable – joints can be stabilized in both deterministic and stochastic simulations. 106 Fig IV.7 presents the results of Monte Carlo dynamic simulations for unloaded task configuration (T1=[50 o 90 o 0N - ]). Subplot (a) shows all the data points in log-log scale of the two-objective space: end- point kinematics variability (Size Hxy ) versus energy consumption (E). The error was quantified by the size of the area circled by the contours that characterize the hand random distribution in a 2-D Cartesian space (Fig IV.7d), and the quantity does not account for the directional difference in its variability. As shown in (a), out of 330 control points ( ▪), four ( □) were identified as the PFs which represent the efficient control of postural accuracy; and many points with high energy consumption were accompanied by relatively large hand variability. Within the Pareto front, higher energy resulted in lower variability, an error-energy trade- off consistent with the single-joint system (Selen et al. 2005; Selen et al. 2006). But proportionally increasing the muscle control pattern of the low-energy PF (point 1) resulted in increased variability (see point 1* - ■). We also plotted the shape and orientation (i.e. pattern) of the hand random distribution (Hxy) in the Cartesian space with respect to the energy consumption (also in log scale) in (b) and (c). There was a large variation of Hxy patterns across the entire range of energy consumption, but relatively little changes among the PFs. These patterns can be viewed in (d) where the hand fluctuations in the 2-D plane with the characteristic contours were plotted for the 4 PFs and point 1*. To examine how the six muscles were coordinated to achieve the postural accuracy at various levels of energy consumption, we plotted the muscle force patterns normalized by total force of the six muscles as stacked bars corresponding to the 4 PFs and point 1* in (h). Comparing point 4 to point 1 and 1*, by increasing the activity of mono-articular elbow flexor (BS), and reducing the activity of bi-articular flexor (Bsh), higher postural accuracy was achieved with slightly greater energy cost compare to point 1, and much cheaper control and more accurate posture than point 1*. The absolute muscle forces are shown in (i), and the muscle activation patterns are shown in (j). Redundant muscles are believed to be used by the CNS for impedance modulation (Osu and Gomi 1999) to counteract external perturbations; we postulated that they may also contribute to the regulation of internal noise. To examine this, the end-point stiffness (Kh) ellipses corresponding to the 4 PFs and point 1* are illustrated in (g). The values at the center of the ellipses quantify the stiffness along the major axis of the ellipses. Their shape and orientation were plotted against energy abscissa with all the other control 107 points in subplots (e) and (f). There was a wide distribution of the Kh patterns across entire energy, but the PFs were constrained in a small range. Different muscle coordination patterns could control the Kh patterns as shown in the five Kh ellipses, but a comparison of Hxy contours and Kh ellipses did not reveal a strong correlation between the patterns of the two. 0 50 100 Fm (%) PC DP Bsh Tlh BS Tlt 0 50 Fm (N) 1 2 3 4 1* 0 1 Act (0-1) (h) 10 2 10 3 10 4 0 50 Shape Kh 80 100 120 θ o Kh 20Nm 25Nm 33Nm 43Nm 139Nm 10 2 10 3 10 4 10 -2 Size Hxy (cm 2 ) 5 10 15 20 Shape Hxy -50 0 50 θ o Hxy Fh=0N 1 2 3 4 1* (a) (e) (b) (c) (f) (d) Hxy distribution (g) Kh ellipses 10 2 10 3 10 4 0 50 Shape Kh 80 100 120 θ o Kh 20Nm 25Nm 33Nm 43Nm 139Nm 10 2 10 3 10 4 10 -2 Size Hxy (cm 2 ) 5 10 15 20 Shape Hxy -50 0 50 θ o Hxy Fh=0N 1 2 3 4 1* (a) (e) (b) (c) (f) (d) Hxy distribution (g) Kh ellipses (i) (j) Figure IV.7: Multi-objective optimization analysis of task configuration 1 (T1) when the error is measured by the size of hand random distribution (Size Hxy ). (a) Log-log plot of Size Hxy versus Energy consumption of all the control points ( ▪), the four identified Pareto frontiers ( □-1,2,3,4) , and one point ( ■-1*) that is the result of proportional scaled control pattern from Pareto frontier 1. The same sets of control points are plotted against same log scale of Energy for (b) shape, (c) orientation of the Hxy, and (e) shape, (f) orientation of the hand stiffness ellipses. (d) shows the hand random distribution (Hxy) traces and the contours of the five points. The four PFs are plotted on the same scale so relative size can be compared, whereas point 1* is in a much larger scale. The task configuration is also illustrated schematically on the first Hxy distribution. (g) shows the end-point stiffness ellipses of the five control points. The stiffness along the major axis of the ellipses is labeled at the center of each ellipsis. (h) shows the muscle control patterns of points 1,2,3,4 and 1* as stacked bar plots of normalized muscle force vectors. The length of each segment represents the relative contribution of each muscle in terms of muscle forces. From bottom up, the colored segments represent the relative contribution of shoulder flexor (PC-yellow), shoulder extensor (DP-green), biarticular flexor (Bsh-black), biarticular extensor (Tlh-blue), elbow flexor (BS-red), and elbow extensor (Tlt-magenta). (i) shows the absolute values of muscle forces. (j) shows the muscle activation levels. 108 d) Muscle redundancy allows task-relevant control of SDN We studied the task-relevant control by examining the Pareto control of end-point variability along two sets of orthogonal axes (DAAs): 1) a1 and a3 in the extra-personal space, and 2) m1 and m2 in the local (or intra-personal) space defined by the major and minor axis of joint manipulability ellipsoid. The measure of variability was constrained as the standard deviation along a single axis. In this section we uses task configuration 1 (T1) to first analyze the task-specified directional control along a1 and a3 axes (Fig IV.8); then the results along m1 and m2 axes are presented to examine the effect of joint geometry on the flexible control of postural variability (Fig IV.9); finally, the findings of the directional control are summarized for all the three task configurations (T1, T2 and T3 in Fig IV.10, Fig IV.11). Extra-personal space control – a1 and a3 axes Fig IV.8 shows the results of postural control along a1 (1 st column) and a3 (2 nd column) axes. Among the same control points from the dynamic simulation in T1 (shown in Fig IV.7), 4 PFs were identified along a1 axis (d1) and 6 PFs along a3 axis (d2). The corresponding size, shape and orientation of hand distribution were plotted with respect to energy consumption in subplot (e), (f) and (g). For a1 axis, as energy increases, the reduction of SD (see (d1)) is accompanied by the reduction of the total variability (Size Hxy in (e1)); this monotonic SD- Size Hxy relationship reveals a less flexible control along a1 axis. However, for axis a3, the reduction of this 1-D variability (SD) (see (d2)) sometimes was accompanied by larger overall variability (Size Hxy in (e2)), which indicates a flexible control along this direction. The flexibility was clearly achieved through modulating the Hxy shape and orientation as seen in (f2), (g2) and (h2). The subplot (a1) and (a2) show the control patterns as the normalized Fm bars for the PFs along a1 and a3 axes. There were large variations in the control patterns among PFs in both axes; and the patterns accompanying the minimal variance along the two axes were very different from each other (see point 4 in (a1) and point 6 in (a2)). To demonstrate the effect of impedance modulation on the end-point variability, we also plotted the shape (i) and orientation (j) of Kh patterns, and Kh ellipses of PFs (k) for both axes. Axis a3 demonstrated 109 larger variations of Kh orientation in comparison to axis a1, but again these variations were not closely related to those of Hxy distribution. Intro-personal space control – m1 and m2 axes One important biomechanical constraint is limb geometry. We took a further step in understanding the directional control under this constraint by performing the similar MOP analysis along major (m1) and minor (m2) axes of arm manipulability ellipsoid. Fig IV.9 (m1 vs. m2) has the same format of data presentation as Fig IV.8 (a1 vs. a3); and the control strategies along the two sets of axes demonstrated a good correspondence (m1 to a1, and m2 to a3). The four PFs identified for SD along m1 (Fig IV.9 (a1)) coincide with those of a1 (Fig IV.8 (a1)); alike a1 axis, there is a relatively rigid control along m1 revealed by a monotonic relationship between SD m1 and Size Hxy ; and similar to a3 axis, the larger fluctuations of SD m2 -Size Hxy curve indicate more flexibility in controlling the hand variability along the minor (m2) axis of manipulability ellipsoid. This correspondence between extra- and intro-space axes can be explained by the similar alignment between the two set of axes, which reveals the constraining effect of the joint geometry on postural control tasks. 110 0 50 100 Fm (%) PC DP Bsh Tlh BS Tlt 0 10 20 Fm (N) 1 2 3 4 5 6 0 1 Act (0-1) 0 50 100 Fm (%) PC DP Bsh Tlh BS Tlt 0 10 Fm (N) 1 2 3 4 0 1 Act (0-1) 10 -2 Size Hxy (cm 2 ) 5 10 15 20 Shape Hxy -50 0 50 θ o Hxy Fh=0N 10 -2 Size Hxy (cm 2 ) 5 10 15 20 Shape Hxy -50 0 50 θ o Hxy Fh=0N 10 2 10 3 10 4 10 -1 SD a1 (cm) 10 2 10 3 10 4 10 -1 SD a3 (cm) 10 2 10 3 10 4 0 50 Shape Kh 80 100 120 θ o Kh 20Nm 25Nm 34Nm 41Nm 10 2 10 3 10 4 0 50 Shape Kh 80 100 120 θ o Kh 20Nm 25Nm 27Nm 35Nm 41Nm 67Nm 1 23 4 1 2 3 4 56 (b1) (b2) (c1) (c2) (d1) (d2) (e1) (e2) (f1) (f2) (g1) (g2) (h1) Hxy distribution (h2) Hxy distribution (k1) Kh ellipses (k2) Kh ellipses (i1) (i2) (j1) (j2) (a1) (a2) a1 axis a3 axis 111 Figure IV.8: Multi-objective optimization analysis of task configuration 1 (T1) when the error is measured by the standard deviation (SD) along a1 (see column 1) and a3 (see column 2) axes. (a1) and (a2) are the muscle control patterns of the Pareto frontiers as stacked bar plots of normalized muscle force vectors. The length of each segment represents the relative contribution of each muscle in terms of muscle forces. From bottom up, the colored segments represent the relative contribution of shoulder flexor (PC-yellow), shoulder extensor (DP-green), biarticular flexor (Bsh-black), biarticular extensor (Tlh-blue), elbow flexor (BS-red), and elbow extensor (Tlt-magenta). (b1) and (b2) are the absolute values of muscle forces. (c1) and (c2) are the muscle activation levels. (d1) and (d2) are the Log-log plot of SD versus Energy consumption of all the control points ( ▪) and the four identified Pareto frontiers ( □). The same sets of control points are plotted against same log scale of Energy for (e) size, (f) shape, (g) orientation of the Hxy, and (i) shape, (j) orientation of the hand stiffness ellipses. (h1) and (h2) are the hand distribution (Hxy) traces and the contours of the Pareto frontiers. They are plotted on the same scale so relative size can be compared. The task configuration is also illustrated schematically on the first Hxy distribution. (k1) and (k2) are the end-point stiffness ellipses of the Pareto frontiers. The stiffness along the major axis of the ellipses is labeled at the center of each ellipsis. 112 0 50 100 Fm (%) PC DP Bsh Tlh BS Tlt 0 10 20 Fm (N) 1 2 3 4 5 6 7 8 0 1 Act (0-1) 0 50 100 Fm (%) PC DP Bsh Tlh BS Tlt 0 10 Fm (N) 1 2 3 4 0 1 Act (0-1) 10 2 10 3 10 4 0 50 Shape Kh 80 100 120 θ o Kh 20Nm 25Nm 27Nm 41Nm 41Nm 49Nm 56Nm 67Nm 10 2 10 3 10 4 0 50 Shape Kh 80 100 120 θ o Kh 20Nm 25Nm 34Nm 41Nm 10 -2 Size Hxy (cm 2 ) 5 10 15 20 Shape Hxy -50 0 50 θ o Hxy Fh=0N 10 -2 Size Hxy (cm 2 ) 5 10 15 20 Shape Hxy -50 0 50 θ o Hxy Fh=0N 10 2 10 3 10 4 10 -2 SD m2 (cm) 10 2 10 3 10 4 10 -1 SD m1 (cm) 1 23 4 1 2 3 4 5 6 (b2) (c2) (d1) (d2) (e1) (e2) (f1) (f2) (g1) (g2) (h1) Hxy distribution (h2) Hxy distribution (k1) Kh ellipses (k2) Kh ellipses (i1) (i2) (j1) (j2) 7 8 (b1) (c1) (a1) (a2) m1 axis m2 axis 113 Figure IV.9: Multi-objective optimization analysis of task configuration 1 (T1) when the error is measured by the standard deviation (SD) along m1 (see column 1) and m2 (see column 2) axes. (a1) and (a2) are the muscle control patterns of the Pareto frontiers as stacked bar plots of normalized muscle force vectors. The length of each segment represents the relative contribution of each muscle in terms of muscle forces. From bottom up, the colored segments represent the relative contribution of shoulder flexor (PC-yellow), shoulder extensor (DP-green), biarticular flexor (Bsh-black), biarticular extensor (Tlh-blue), elbow flexor (BS-red), and elbow extensor (Tlt-magenta). (b1) and (b2) are the absolute values of muscle forces. (c1) and (c2) are the muscle activation levels. (d1) and (d2) are the Log-log plot of SD versus Energy consumption of all the control points ( ▪) and the four identified Pareto frontiers ( □). The same sets of control points are plotted against same log scale of Energy for (e) size, (f) shape, (g) orientation of the Hxy, and (i) shape, (j) orientation of the hand stiffness ellipses. (h1) and (h2) are the hand distribution (Hxy) traces and the contours of the Pareto frontiers. They are plotted on the same scale so relative size can be compared. The task configuration is also illustrated schematically on the first Hxy distribution. (k1) and (k2) are the end-point stiffness ellipses of the Pareto frontiers. The stiffness along the major axis of the ellipses is labeled at the center of each ellipsis. Directional control in other task configurations We conducted the same set of MOP analysis for directional control on the other two task configurations. Fig IV.10 summarizes the major findings of the three tasks (T1,T2 and T3). The first column depicts the three task configurations with arm manipulability ellipsoids and the two sets of orthogonal axes (a1, a3 and m1, m2). From above analysis, SD-Size Hxy relationship is a good indication of flexible control along a specified direction of interest. Thus for each DAA set, we plotted the SD-Size Hxy curves for each of the two orthogonal axes, and the control patterns expressed as the Fm bars corresponding to the minimal SD along the two axes. We noticed that in another unloaded posture task (T2 =[45 o 50 o 0N - ]), there were more PFs along m1axis compared to those in T1, and the larger fluctuation in SD m1 -Size Hxy curve indicates a more flexible control along m1 axis in T2 than that in T1. When comparing m2 to m1 axis in T2 configuration, there are slightly larger fluctuations of m2 than m1 implying a more flexible control on minor axis of the manipulability ellipsoid, a trend consistent with T1 configuration. The minimal variance along m2 was achieved through a very different control pattern from that along m1, a general feature observed in most DAA sets of all the three tasks; but there is an exception for a1 versus a3 in T2. Here, same control patterns along a1 and a3 axes might be explained by their similar relative orientation with 114 respect to m1 and m2. In the case of postural task with non-zero external load (T3), there was a mono-tonic SD- Size Hxy curve along the direction of external force (a3), but a large fluctuation perpendicular to it (a1); and SD- Size Hxy curve along the m1 showed a much larger fluctuation than that along m2, a finding contrary to those in unloaded condition (T1 and T2). This indicates that external forces may override the joint geometry in constraining the postural control. To determine if the CNS could control the pattern of end-point variability though impedance modulation, we plotted the pattern (i.e. shape and orientation) of hand random distribution (Hxy) against those of equilibrium stiffness ellipses (Kh) of all the control points for each task (Fig IV.11). No strong correlation was found in either the shape (column 1) or the orientation (column 2) between Hxy and Kh, implying factors beyond impedance considerations may dominate the control of SDN. SD a1 Hxy Size SD a3 Hxy Size a1 a3 Fm Task DAA set 1 DAA set 2 SD a1 Hxy Size SD a3 Hxy Size a1 a3 Fm T1: a1 a3 Fm m1 m2 Fm SD a1 Hxy Size SD a3 Hxy Size a1 a3 Fm T3: T2: (42 o 50 o ) 0N a1 a3 m1 m2 (42 o 50 o ) 3N ∠-90 o a1 a3 m1 m2 SD m1 Hxy Size SD m2 Hxy Size m1 m2 Fm 0N (50 o 90 o ) a1 a3 m1 m2 SD m1 Hxy Size SD m2 Hxy Size m1 m2 Fm SD m1 Hxy Size SD m2 Hxy Size m1 m2 Fm Figure IV.10: Summary of directional control of the three task configuration: T1, T2, T3. First column shows the schematic diagram of the three task configurations. Their directional controls along the two sets of DAAs are presented in the same row. For each DAA set, the Pareto frontiers identified along each of the two orthogonal axes are presented as the Hxy size versus SD curves. The marker type represents the axis along which the SD is measured: ○-a1, □-a3, +-m1 and x-m2. The control patterns that correspond to the minimal SD (denoted as red marker in Hxy-size vs SD curves) along the two orthogonal axes are shown as the stacked bar plot. 115 T1: T2: T3: Shape Orientation Kh Hxy Kh Hxy Kh Hxy Kh Hxy Kh Hxy Kh Hxy Figure IV.11: Plots of the shape (column 1) and orientation (column2) of Hxy distribution versus end-point stiffness ellipses for the three task configurations. v. DISCUSSION Both experimental data and musculoskeletal models suggest that “redundant” muscles enable stable and flexible control of motor tasks in unstable environments (Burdet et al. 2001; Danna-dos-Santos et al. 2007; Osu and Gomi 1999). This work provided computational evidence that, even for externally unperturbed motor tasks, muscle “redundancy” is utilized by the CNS to maximize task performance in the presence of signal-dependent noise. We found that appropriate coordination of multiple muscles would be useful to achieve both energetic efficiency and postural accuracy when measured along a single, task- specified axis. No strong correlation was found between the patterns (shape, orientation) of end-point variability and equilibrium hand stiffness, a commonly used indicator of end-point impedance. Furthermore, the best control strategies were highly dependent on posture and external loading conditions. Thus, there seems to be no simple rule such as impedance control or co-contraction that can be used by the brain or by researchers to compute an optimal solution for multi-muscle, multi-joint systems. The optimal solution for any given performance criteria depends on complex details of musculoskeletal architecture such as the variation of moment arms with posture. Maintaining hand position is important for many motor tasks in daily living; to understand fully how multiple muscles are utilized for such a goal, a 3-dimensional arm 116 model with a complete set of muscles will be needed in future studies. Neurophysiological experiments on this system may need to collect data systematically from all or most of these muscles. This study extends the analysis of the effects of SDN in a uniarticular system operated by a simple pair of antagonist muscles (Selen et al. 2005). Unlike their simple model and performance criterion, we adopted a complex musculo-skeletal model of human arm with multiple articulation, realistic musculoskeletal architecture, and muscle models (VM) with physiological architectural parameters and physiological properties (Song et al. 2008a). The goal of the task was to stabilize the hand against inherent SDN in the motor system with minimal metabolic cost. The realism of the system could help reveal many details in the human motor system that might place important constraints on the CNS when selecting a desired strategy. Firstly, accurate modeling of muscle-joint geometry, i.e. muscle moment arm – angle relationship is important, because muscle MA variation affects its stabilizing nature at a given joint angle (Hogan 1991; Shadmehr and Arbib 1992). When deciding which muscles to be activated or co-activated, the CNS might have the knowledge of musculoskeletal geometry, and preferably chooses the muscles with positive MA(q) slopes, thus stabilizing the joint with less required co-contraction and energy cost (Bunderson et al. 2008). Secondly, muscle contractile properties such as force-length relationship, and muscle strength (by maximal isometric force) need to be represented accurately, because they determine its moment generating capabilities, the effects of the motor noise and energy consumption. Thirdly, we chose to explore a two- joint system, not only because it is more close to daily motor tasks, but also because it places many additional elements that are all crucial to the motor control. For example, the nonlinear intersegmental dynamics may affect the noise propagation in the musculo-skeletal systems; the joint geometry determines the sensitivity of the end-point kinematics to the variability in the joint space (e.g. the manipulability ellipsoid)(Sabes and Jordan 1997; Sabes et al. 1998) thus affect the distribution of hand variability (Fig IV.9, IV.10); and the mono- and bi-articular muscles may have differential functions in controlling postural stability (Lan 2002). Finally, we chose isometric postural tasks to study the effect of SDN as opposed to many multi-joint motor coordination studies on SDN that used reaching movement tasks (Franklin et al. 2007; Harris and Wolpert 1998; Nagata et al. 2002; Todorov and Jordan 2002). The purpose was to prevent large interaction torques in the movement from masking the relatively subtle effect of SDN on postural 117 accuracy. Although seemingly simple, the postural control was rather sensitive to many parameters under control: we varied the equilibrium joint angles and external loading conditions, all of which factor into the control strategies available to the CNS. One important message from this study is that, appropriate activation patterns of multiple muscles is the key determinant for constraining the detrimental effect of SDN in a multi-joint system. In a single joint system with a pair of antagonistic muscles, a unique solution of muscle activations can be determined by the required posture, net torque, and joint impedance requirement. To constrain the effect of SDN on single-joint kinematics, simply scaling up the co-activity has been shown to be effective both in previous modeling work by (Selen et al. 2005) and the present simulation results (Fig IV.5). But such a strategy did not apply to a multi-joint system with mixed mono- and bi-articular muscles. By proportionally increasing the muscle activation levels, our results showed a larger postural variability with higher energy consumption (point 1* in Fig IV.7a). And the best strategy with minimal variability (point 4 in Fig IV.7) demonstrated a very different control pattern (Fig IV.7h). This finding is quite revolutionary since many psychophysical studies have demonstrated that humans tend to stiffen their limbs when the motor tasks demand accurate control or when destabilizing factors are present (Franklin et al. 2007; Franklin et al. 2004; Gribble et al. 2003; Mussa-Ivaldi et al. 1985; Osu et al. 2004; Selen et al. 2006). For instance, when subjects produced point-to-point reaching tasks, (Gribble et al. 2003) observed a trend of increased muscle co-contraction with reduced target sizes in seven mono- and bi-articular muscles, reflected by their higher EMG activity. And after random force perturbation, (Mussa-Ivaldi et al. 1985) found that the subjects increased the overall stiffness in preparation for stabilizing the limb against potential perturbations. However, these energetically costly strategies were mostly used to combat the perturbing effects of joint interaction torques (Gribble et al. 2003; Gribble and Ostry 1999; Koshland et al. 2000), or the external forces (Burdet et al. 2006; Franklin et al. 2007; Franklin et al. 2004; Mussa-Ivaldi et al. 1985; Thoroughman and Shadmehr 1999). In the postural control tasks where these factors are almost absent, or negligible, our simulation results showed that minimizing the effect of SDN on the postural variability did not necessarily satisfy energy efficiency. Instead intelligently recruiting available mono- and biarticular 118 muscles to achieve a desired level of local stiffness as proposed in (Bunderson et al. 2008) tended to minimize kinematic noise at relatively low levels of cocontraction and energy consumption. Another important finding is that the redundant muscles can be coordinated in such a way that the variability along a single, task-specified axis is minimized, whereas that along other axes remains relatively large. Such a strategy can be associated with less energy cost (Fig IV.8, IV.9, and IV.10). This behavior conforms to a general observation in the redundant motor systems: the fluctuations of redundant effectors are not equally constrained; rather, they are coordinated in a way such that only the variability that hurts the performance is limited. Such a structure of motor variability is believed to arise from a “minimal intervention principle” adopted in the CNS, that is the control effort is only taken to correct deviations that affect task performances (Todorov and Jordan 2002). Such a strategy is of advantage because it saves control energy, at the same time “ensures both stability of important performance variables and flexibility of patterns to deal with other task components and possible perturbations” (Latash et al. 2007). The means to achieve these benefits is to form the synergistic actions of elementary variables and proper co-variation patterns (e.g. muscle coordination patterns) so that the variability or noise in the redundant elements is largely projected into an “uncontrolled manifold” (Scholz et al. 2003; Scholz and Schoner 1999) (e.g. the dimension orthogonal to the task-specified axis for the postural accuracy) of the high-dimension state space. This sub-manifold, in linear algebra analysis, is essentially the null space of the Jacobian matrix (J(u)) that relates the small deviation of the high-dimension elemental variables (u=[u i ], i=1,2,…,m) to the error in the task performance space with lower dimension (p(u)=[p j (u), j=1,2.,…,n], n<m). To get co-variation that is beneficial for stabilization of the task variables, the CNS may be aware (Goodman and Latash 2006) or quickly learn [through a neuronal network (Bullock et al. 1993)] this coordinate transformation, in effect developing a forward mechanism that selects one solution out of a set of “motor-equivalent” solutions. In the case of postural variability due to SDN, such a “coordinate transformation” is through dynamic interaction of biomechanical limbs and redundant muscles, thus our system Jacobian is “dynamic”. To get muscle co-activation patterns that selectively control postural variability in a specified direction, the CNS must have the knowledge of this “dynamic Jacobian”. 119 Determining the dynamical Jacobian is not trivial because it arises from all the complexities of musculoskeletal architecture and dynamics, some of which are non-stationary (e.g. muscle potentiation and fatigue). Even in our relatively simple and stationary system, the simulation results showed that the best control strategies were highly sensitive to the arm postures and external loads (Fig IV.10).The flexible control of postural accuracy along some task-specified axis was constrained by arm geometry, which determines the sensitivity of the end-point kinematics to joint noise (the manipulability ellipsoid). In T1 task configuration (T1 in Fig IV.10), when the control was along or close to the direction of maximal sensitivity, i.e. the principle axis of the manipulability ellipsoid, the neuromuscular control could hardly override the passive limb properties to achieve flexible control of limb posture. In such a case, the positional accuracy along one axis was always accompanied by reduction of the overall variability (indicated by a mono-tonic SD-Size Hxy curve along a1 and m1 in Fig IV.10 T1). Limb biomechanical properties of such kind have been found to constrain the control strategies in many motor tasks: Valero- Cuevas (Valero-Cuevas 2000) has found that the redundant muscles are needed for dexterous control constrained by finger mechanical properties; Sabes and Jordan (Sabes and Jordan 1997; Sabes et al. 1998) identified a sensitivity model of motor planning that takes into account limb inertia and geometrical properties to avoid collision on obstacles during reaching movements. However, when the task involves force production, the external loads seemed to be a dominant constraint on top of biomechanical properties that constrain posture control: there was a much more flexible control along the direction perpendicular to the force vector than along it (Fig IV.10 T3: a1 vs. a3). Here we demonstrated the control dependence on arm postures and external loads through very limited task configurations (two postures and one end-point force). Further studies on extensive joint configurations and external loading conditions are needed to establish possibly principled relationships between the two. Learning of this “dynamical Jacobian” is further confounded by the finding that there was no close correlation between patterns of end-point variability and the stiffness ellipses (Fig IV.11); thus a simple mapping of stability requirements with impedance modulation for external perturbations can not easily extend to internal fluctuations. Impedance control is often assumed to be closely related to the stiffness ellipses that can be determined experimentally in multi-joint motor tasks (Franklin et al. 2004; Mussa- 120 Ivaldi et al. 1985). Increasing stiffness has generally been useful to suppress internal noise similarly to its ability to reduce the kinematic effects of external perturbations (Osu et al. 2004; Selen et al. 2005; Selen et al. 2006). Because the patterns of the stiffnesses ellipses can be adapted to the directions of environmental disturbances, it is natural to assume that the patterns of stiffness ellipses will also shape the variability of the end-point due to internal noise – the variability along the major axis of the ellipses is smallest. Our simulation results, however, contradict this assumption. Even though we have a large range of variation of stiffness orientation and shape changes in our models, these factors could not predict the patterns of end- point variability (Fig IV.7, IV.8, IV.9). This implies that the CNS might separately control the internal fluctuation and external perturbation, and the stiffness that is derived from deterministic motor command (without noise) can not be used to infer effect of internal noise. Some limitations of this work also have to be considered. The noise model used assumed a linear scaling of noise level with the magnitude muscle force (i.e. constant CV), a signal-dependent feature that has been well explored and documented (Fuglevand et al. 1993; Galganski et al. 1993; Jones et al. 2002; Laidlaw et al. 2000; Pastor et al. 1991; Schmidt et al. 1979; Selen et al. 2005; Slifkin and Newell 1999). But some recent studies have showed that CV may not always be constant; it may decrease exponentially over the force range (Moritz et al. 2005; Taylor et al. 2003). Further, we did not have a model of motor unit synchronization, but it has been reported that when motor unit synchronization occurs at high co-activation levels, it possibly leads to abrupt increases in force variability, which becomes apparent as tremor (McAuley et al. 1997). Conversely, synchronization was suggested to be able to dampen the effect of discharge-rate variability on force fluctuations (Taylor et al. 2003). These findings imply a rather complicated picture of the motor noise, which will further confound the search for an optimal control strategy. However, if CV decreases at higher levels of force, the reduction of variability with muscle co- contraction will still be accompanied by high energy consumption, making it an unattractive strategy given performance criteria that balance postural variability and efficiency. Nonetheless, future research with more accurate noise models should be considered to understand specific control strategies at various levels of force production tasks. 121 Caution should also be taken when inferring the impedance control directly from the measured joint and end-point stiffness. In addition to exhibiting spring-like properties as a result of the FL relationship, muscles provide viscous damping as a result of their FV relationships. To control the dynamical response of a system to either the fluctuations of its motor inputs or external perturbations, both of these components need to be modulated (Sciavicco and Sicilliano). In musculoskeletal system, damping and stiffness are coupled in the muscle intrinsic contractile properties, i.e. increased stiffness is often accompanied with increased viscosity, and as a result stiffness ellipsoid has been commonly depicted in modeling and experimental studies to infer impedance control (Burdet et al. 2001; Burdet et al. 2000; Burdet et al. 2006; Franklin et al. 2007; Franklin et al. 2004; Osu and Gomi 1999; Selen et al. 2005; Selen et al. 2006). But although the two components co-vary, they are not likely to scale proportionally with each other due to the non-linearity of muscle contractile properties. For small, high frequency perturbations such as those associated with motor noise, the FV properties of muscle play a dominant role in stabilizing the limb compared to the relatively flat FL relationship that obtains for postures near the middle of ranges of motion (Brown and Loeb 2000c). Thus there is a tendency to misrepresent impedance modulation if only stiffness ellipses are considered. Future work with a complete measure of limb impedance properties will be needed to draw decisive conclusions regarding the relationship between impedance modulation and noise control. Our model of postural controller was an open-loop system that did not account for reflexes. In reality, the autogenic stretch reflex (or short-latency reflex) adds to the intrinsic muscle stiffness to help provide limb stability (Allum and Mauritz 1984). Simulations using a sensorimotor system model (Song et al. 2008a) that incorporated the spindle and Golgi-tendon feedback loops revealed such an effect: with closed loop servo mechanisms, lower level co-contraction resulted in smaller variability caused by the internal SDN (Song et al. 2006). Further research into the function of proprioceptive feedback is needed to identify optimal control policies when there is noise in both the motor output and sensory encoding [such as the noise described in muscle spindle afferents, which depends complexly on fusimotor control (Loeb 1984; Loeb and Marks 1985)]. However, the current study provided a starting estimation of how the common path way that integrated both central and peripheral information affects motor variability and energy consumption. 122 vi. CONCLUSION In summary, we have shown that signal-dependent noise in a multiarticular system affects postural stability in complex ways that depend on many details of musculoskeletal architecture and physiology. Selective recruitment of multiple muscles can maximize the task performance while minimizing metabolic load, but this optimal solution strongly depends on the posture and external loading conditions. This study is the first attempt to explore how motor noise propagates in such a complex and realistic limb. It sets up a starting point for further investigation on how multiple features in human motor system – i.e. noise (sensory and motor), energy, redundancy, muscle and joint dynamics – combine to determine the optimal control strategies for multi-joint movements. It remains to be determined experimentally how well the central nervous system can actually deal with these complexities to optimize performance under various conditions. 123 CHAPTER V: SUMMARY AND FUTURE WORK i. SUMMARY Models of the musculoskeletal system are used in a broad range of scientific studies; the levels of their complexity and focus of their details vary widely depending on the questions to be addressed or the specific applications. In this dissertation, we developed a sensorimotor systems model of the human upper extremity – virtual arm (VA) – in order to understand how its mechanical dynamics influence the choice of control strategies for posture and reaching movements (Song et al. 2008a) (Chap III). “Virtual” in the model name reflects the realism and completeness of the sensorimotor representation: each component of the sensorimotor system was modeled in close accordance to its physiological properties that were identified and verified in a series of experimental and modeling work (Brown et al. 1999; Brown and Loeb 1999, 2000a, 2000b; Cheng et al. 2000; Mileusnic et al. 2006; Song et al. 2008b). The system level behavior of the integrated model should be consistent with the gross sensorimotor physiology under varying operating conditions and neural control inputs. VA is an advance over previous models of upper extremity in that it includes both sensory and motor functions; the integrated models of spindle and GTO afferents can facilitate understanding of the functions of proprioception and sensorimotor integration. On the motor side, virtual muscle (VM) is to date the most physiologically realistic muscle model (Cheng et al. 2000). The recent design of a new recruitment method and S-function implementation (Song et al. 2008b) (Chap II) greatly improves its computational efficiency, and makes it suitable for comprehensive simulation study in motor control of multi-joint multi-muscle systems. Combining VM with accurate models of the postural dependence of the muscle moment arms allows precise estimation of moment generating capability of each muscletendon actuator in the arm model. Overall, the model is a useful tool for computational studies of sensorimotor coordination and motor learning. The modular structure allows easy modification to study various neuropathological conditions, such as deafferentation, stroke and spinal cord injury, and to address the strategies of neuro-rehabilitation. The VA model was used for the computational study of control strategies for dealing with signal dependent motor noise (Chap IV). This study integrated a number of features (e.g. redundancy, noise, 124 intersegmental dynamics) and control objectives (energy, accuracy) of human redundant and noisy motor systems that were often ignored in the previous studies due to the simple design of the task or obscured by many interacting factors. In addition, most previous computational analyses of motor coordination started by enforcing certain rules (such as optimization) that are assumed to reflect limits on the computational strategies or controllable parameters available to the brain. The current study, however, first identified the achievable states of the musculoskeletal plants across the entire control space and then identified strategies that might be selected to achieve desirable tradeoffs between kinematic and energetic performance criteria. This approach allows us to visualize various biomechanical constraints that are otherwise buried in the automated tuning process (e.g. optimization). Exploring this space for a wide range of simulated tasks helps to segment the effects of motor noise from other potentially confounding factors (e.g. intersegmental dynamics and external forces). The conclusions are rather revolutionary: to achieve postural stability, the CNS should consider separately internal noise and external perturbations, rather than addressing them together through modulation of the system impedance, a strategy often suggested or assumed by many previous motor control studies. Compared to the previous efforts in studying motor coordination in the face of SDN, the current work has addressed a more complex and realistic musculoskeletal system, revealing aspects of control that are qualitatively different from those apparent in highly simplified systems. It is only an initial step, however, in identifying and exploring the range of system and task parameters that influence the achievable performance criteria and their associated control strategies. ii. FUTURE WORK a) Generalization to 3-D musculoskeletal systems The current modeling work was constrained kinematically in a 2-D horizontal plane with only shoulder and elbow joints. This simplification of the complete limb has been employed in many other modeling studies (Lan 1997; Todorov and Jordan 2002) to simplify the specification and computation of the model system and because it comports with commonly adopted experimental configurations (Franklin et al. 2007; Osu et al. 2003; Osu and Gomi 1999). But we live in a 3-D world where any motor task requires 125 appropriate positioning and orientating of the multi-joint limbs. Moving the arm in a horizontal plane without benefit of a frictionless support surface entails activating many muscles to combat gravitational force, some of which are involved also in the planar movements. For example, shoulder abductor muscles such as the anterior deltoid might need to be recruited. However, the shoulder girdle is an elaborate 3-D articulation that is actuated by complex muscle groups (rotator cuff muscles). Each of the muscles has varying moment arms that cross multiple DOFs. For instance, the anterior deltoid functions as both shoulder abductor and flexor; as such, the torques produced for joint abduction at the frontal plane would also interfere with holding the posture in the horizontal plane. A high level of force required to elevate the arm would be associated with proportionately high levels of motor noise, which would affect the ability of this muscle to contribute to fine control of abduction for the movement of interest. To perform simultaneously the “multiple” functions (producing gravitational torque and controlling arm posture), the complete set of muscles around each joint are needed, and “actively” holding the arm at the horizontal (or any) plane is expected to present another constraint in addition to noise and external perturbations that the CNS must consider in selecting optimal control patterns. Extracting a separate static controller that maps the joint states (joint configuration) and gravitational force (Guigon et al. 2007) is possible but it may obscure interactions between postures, moment arms and force-length relationships of individual muscles (Garner and Pandy 2001) that affect their performance in the dynamic aspects of the task. b) Experimental testing of optimal control hypotheses The identification of optimal control strategies in modeled systems is useful primarily in comparison to the strategies that human subjects actually employ. Such comparisons are meaningful only if the experimental task meets certain criteria: 1. The mechanical configuration and constraints on the experimental task must be accurately represented in the modeled system. 2. The subject performing the experiment must be aware of and motivated to optimize the same mix of performance criteria that were used to constrain and identify the optimal solutions for the model. 126 3. The subjects performing the experiment must be provided with knowledge of how well they meet the performance criteria each trial and with sufficient time to practice the task to be able to discover the optimal solution, which is usually iteratively by trial-and-error. From the major findings of the current computational study, we hypothesize that in human multi-joint postural control tasks, end-point kinematic variability due to internal motor fluctuations (e.g, SDN) possesses geometric features that can be modulated flexibly through selective recruitment of redundant musculature; the optimal control strategies can be learned through practice in achieving specific task goals. To test this hypothesis on human subjects we would need to design experiment to examine the effect of different activation patterns of multiple arm muscles on the spatial characteristics of end-point variability. Though relatively simple to implement in a modeling study, such a testing on human subjects involves a number of complications that could potentially lead to misinterpretations of the data or even erroneous findings. Here, we will point out the potential constraining factors and propose corresponding approaches according to the aforementioned three criteria; the purpose is to circumvent the potential complications or minimize their detrimental effect on data collection and interpretations. In this modeling work, human arm was kinematically constrained in the horizontal plane at the shoulder level, where no active contraction of muscles was needed to abduct the shoulder joint to this configuration. This has to be implemented in the experiments through providing some form of mechanical support for the arm to sustain it in the horizontal plane without need for any muscle efforts. Furthermore, the arm was not affected by any passive joint forces or frictional disturbances. The mechanical support should also replicate this condition. We are also interested in how the impedance characteristics might affect the control of motor variability. To measure limb impedance, a robotic manipulandum is needed to provide small position disturbances. The challenge is to provide well-controlled displacements and accurate force measurements, at the same time to minimize the nonlinear forces from manipulandum dynamics that might interfere with the hand dynamics. Sophisticated instrumentation satisfying these requirements has been developed (Gomi and Kawato 1997). A PFM (parallel link drive air-magnet floating manipulandum) provides the support of the human arm on a horizontal plane to be free from the force of gravity and to reduce fatigue, the air-magnet floating 127 mechanism is friction-free, and the nonlinear forces applied at hand due to manipulandum dynamics is reduced through high frequency digital signal processor and an active control mechanism. The manipulandum slightly pushes and pulls the hand (6-8 mm) in eight randomized directions, and the restoring forces are recorded. To prevent the occurrence of voluntary response, the duration of displacement is kept short (0.12-0.3 sec). The limb impedance characteristics are then estimated through modeled arm dynamics and a linear regression method [details can be found in (Gomi and Kawato 1997)]. To examine whether humans behaves in a same way with the model, same motivation and task specification as those in the simulation have to be presented to the subjects during the experiments. Without motivation, subjects could choose a stereotypical control pattern that does not conform to the task requirements. The lack of task-related instruction has led to some misinterpretation of stiffness regulation mechanisms in free-posture tasks (Flash and Mussa-Ivaldi 1990; Mussa-Ivaldi et al. 1985). The optimal solution defined in the simulation study is based on multiple performance indices: the postural accuracy and energy efficiency. Accordingly, to identify whether humans adopt the same optimal strategies, identical mixed criteria has to be instructed to the subjects during the trial. Energy criteria could be instructed through telling the subjects to minimize the control effort and at the same time to display an energy bar on a visual monitor indicating the energy consumption based on muscle EMG activity. Here the computation of energy from EMG activity follows the equation used in the modeling study: E= ∑ i (EMG i *PCSA i ) where i index the individual muscle. EMG levels will have to be normalized to those recorded during maximal voluntary contractions. Subjects will be encouraged to minimize the length of the bar while satisfying the postural requirements. As to the postural accuracy, the criterion is measured according to the task requirement. If the task goal is to maintain the hand position without a consideration of directional accuracy, the error cost can be measured by the area subsumed by the fluctuations of hand position (Size Hxy ) in the simulation; this can be implemented in the experiments by providing a target with even geometric configuration, for example a circle with certain radius that denotes the required level of accuracy. Most previous behavioral studies on co-contraction and motor variability controlled the accuracy requirement using a circular target (Gribble et al. 2003). This was informative regarding the general trend of co-contraction on the motor accuracy, but could not reflect the effect of muscle coordination on the 128 spatial patterns of end-point variability. Many motor tasks in human daily activities, however, require accuracy along a specific direction or subspace. We can convey such postural criteria in the experimental setting by using a geometric target with uneven spatial dimensions. In a 2-D horizontal plane, ellipsoid or rectangle shape can be shown on the monitor to specify the direction that has the most stringent accuracy requirement for a given task. Subjects will be instructed to maintain the hand within this geometric target while minimizing the energy bar. The learning process will be monitored by the EMG levels and kinematic performances. To examine the controllability of the end-point variability through muscle coordination, separate trials can be designed in which the subjects are asked to selectively recruit groups of muscles at the various levels while monitoring the kinematic variability of the hand. One effective way to enforce specific coordination pattern is through real-time feedback of muscle EMG activity. Subjects will be asked to coordinate the muscle activities that meet the specified pattern shown as reference EMG bar graphs on a computer monitor, at the same time, the levels of energy expenditure will be informed to the subjects through the energy bar that is also used for the optimization trials. Both controllability and optimization trials require recording the subtle effect of SDN on end-point variability. An accurate and fast motion capture system is needed. An active marker-based motion track system (e.g. Optotrak Certus) with high position resolution (as small as 0.01 mm) and high speed tracking (sampling rate up to 4600 Hz) can be used for this purpose. The largest challenge in this experiment, probably in most psychophysical studies of motor control, is to accurately characterize the muscle activities and their mechanical contributions to the prescribed motor tasks. Most commonly used technique is EMG recording, but actual torque at a joint is determined by multiple factors (neural excitation, muscle length, velocity, moment arm and history of use). Furthermore, the estimation of the neural activation commands from the EMG signals is also questionable due to a number of technical limitations (Loeb and Gans, 1986). Firstly, typical skin surface recordings tend to over-sample the activities of only superficial subgroups of muscle fibers (typically fast-twitch fibers); they underestimate overall muscle activities at lower force levels where only the deeper-lying, slow-twitch fibers are recruited. Secondly, we are aiming at finding the recruitment patterns of a complete set of arm 129 muscles; the potential cross-talk, movement artifacts, or misplacement of the electrodes can considerably affect the data analysis and interpretations for our study. To date, no technique has been developed to resolve all of these complications. A recent study that combined a comprehensive musculoskeletal model and experimental recordings of the monkey single-joint movements has provided a systematic way to assess the completeness of EMG signals and their relevance to the motor functions (Cheng and Loeb 2008). In light of this study, we propose a modeling assisted experimental paradigm. In this paradigm, the recorded EMG envelopes are extracted as the activation signals to drive the simulation of the observed joint kinetics and kinematics under varying movement and loading conditions; the muscle actions at various levels of musculo-skeletal motor output (joint torque, joint position and angular velocity) would then be compared to observed motor behaviors to validate the recorded EMG signals. If substantial disparity is observed, fine tuning the electrode placement and recording parameters might be needed before the EMG data are used to infer control strategies. However, another challenge emerges from the validity of the model itself because of the inter-subject variability in the muscle and skeletal parameters. Muscle force outputs are very sensitive to architectural parameters (e.g. optimal fascicle lengths, muscle PCSA, specific tensions, fiber type distributions etc.), and the moment arm variation also determines its moment generating capabilities. In consideration of these, exact matching between simulation results and experimental measurements is not expected, but the qualitative agreement provides confidence for any conclusions from the experimental findings. Searching for the optimal solution given a particular set of task criteria requires a period of practice and learning process. In this process, sufficient knowledge of how well they meet the performance criteria should be provided for each trial and enough number of trials and learning sessions are also needed. In each learning session, warming-up trials are needed to potentiate the muscles. This is because the muscle contractile properties are significantly different before potentiation (Brown and Loeb 1999). Also the duration of the session and the intersession intervals should be specified to provide a reasonable learning time and at the same time to avoid muscle fatigue. For each trial, subjects will be asked to maintain the hand position according to the task specification (along certain axis or overall accuracy, produce or not to produce an external force vector) for a designated period (e.g. 15 seconds according to the simulation 130 study). At the end of the trial, the average energy consumption and the size of the error cost will be fed back to the subjects as an incentive for improvement in the next trial. Hand postural control is critical for many motor tasks, and our simulation study has demonstrated the opportunity for selective control of the various muscles of the arm to structure its variability. With the advances of instrumentation and testing methodology, it is increasingly possible to design sophisticated experiments to demonstrate to what extent the CNS takes advantage of this opportunity. Learning from the human motor system will inform the development of robotic or prosthetic control algorithms that emulate elegant biological motor behaviors. 131 BIBLIOGROPHY Abend W, Bizzi E, and Morasso P. Human Arm Trajectory Formation. Brain 105: 331-348, 1982. Allum JHJ and Mauritz KH. Compensation For Intrinsic Muscle-Stiffness By Short-Latency Reflexes In Human Triceps Surae Muscles. Journal Of Neurophysiology 52: 797-818, 1984. Alstermark B, Lan N, and Pettersson L-G. Building a realistic neuronal model that simulates multi-joint arm and hand movements in 3D-space. HFSP Journal, accepted, 2007. Anderson FC and Pandy MG. Dynamic optimization of human walking. Journal Of Biomechanical Engineering-Transactions Of The Asme 123: 381-390, 2001. Andersson GBJ and Winters JM. Role of muscle in postural tasks: spinal loading and postural stability. In: Multiple muscle systems: Biomechanics and movement organization, edited by Winters JM and Woo SL-Y: Springer, 1991, p. 377-395. Astrom KJ and Wittenmark B. Adaptive control. MA: Addison-Wesley, 1995. Atkeson CG and Hollerbach JM. Kinematic Features Of Unrestrained Vertical Arm Movements. Journal Of Neuroscience 5: 2318-2330, 1985. Balasubramaniam R, Riley MA, and Turvey MT. Specificity of postural sway to the demands of a precision task. Gait & Posture 11: 12-24, 2000. Bernstein A. The coordination and regulation of movements. Oxford, UK: Pergamon Press, 1967. Bernstein A. On dexterity and its development. Mahwah, NJ: Lawrence Erlbarm Associates, 1996. Berthouze L and Lungarella M. Motor skill acquisition under environmental perturbations: On the necessity of alternate freezing and freeing of degrees of freedom. Adaptive Behavior 12: 47-64, 2004. Blemker SS and Delp SL. Three-dimensional representation of complex muscle architectures and geometries. Annals Of Biomedical Engineering 33: 661-673, 2005. Brown IE, Cheng EJ, and Loeb GE. Measured and modeled properties of mammalian skeletal muscle. II. The effects of stimulus frequency on force-length and force-velocity relationships. Journal Of Muscle Research And Cell Motility 20: 627-643, 1999. Brown IE and Loeb GE. Measured and modeled properties of mammalian skeletal muscle. I. The effects of post-activation potentiation on the time course and velocity dependencies of force production. Journal of Muscle Research and Cell Motility 20: 443-456, 1999. 132 Brown IE and Loeb GE. Measured and modeled properties of mammalian skeletal muscle: III. The effects of stimulus frequency on stretch-induced force enhancement and shortening-induced force depression. Journal of Muscle Research and Cell Motility 21: 21-31, 2000a. Brown IE and Loeb GE. Measured and modeled properties of mammalian skeletal muscle: IV. Dynamics of activation and deactivation. Journal of Muscle Research and Cell Motility 21: 33-47, 2000b. Brown IE and Loeb GE. A reductionist approach to creating and using neuromusculoskeletal models. In: Biomechanics and neural control of posture and movement, edited by Winters JM and Crago PE: Springer, 2000c, p. 164-176. Brown LE, Rosenbaum DA, and Sainburg RL. Limb position drift: Implications for control of posture and movement. Journal Of Neurophysiology 90: 3105-3118, 2003. Bryson A and Ho Y. Applied Optimal Control. New York: Taylor and Francis, 1975. Bullock D, Grossberg S, and Guenther FH. A Self-Organizing Neural Model Of Motor Equivalent Reaching And Tool Use By A Multijoint Arm. Journal Of Cognitive Neuroscience 5: 408-435, 1993. Bunderson NE, Burkholder TJ, and Ting LH. Reduction of neuromuscular redundancy for postural force generation using an intrinsic stability criterion. J Biomech, 2008. Burdet E, Osu R, Franklin DW, Milner TE, and Kawato M. The central nervous system stabilizes unstable dynamics by learning optimal impedance. Nature 414: 446-449, 2001. Burdet E, Osu R, Franklin DW, Yoshioka T, Milner TE, and Kawato M. A method for measuring endpoint stiffness during multi-joint arm movements. Journal Of Biomechanics 33: 1705-1709, 2000. Burdet E, Tee KP, Mareels I, Milner TE, Chew CM, Franklin DW, Osu R, and Kawato M. Stability and motor adaptation in human arm movements. Biological Cybernetics 94: 20-32, 2006. Capaday C and Cooke JD. Vibration-induced changes in movement-related EMG activity in humans. Experimental Brain Research 52: 139-146, 1983. Chan SS and Moran DW. Computational model of a primate arm: from hand position to joint angles, joint torques and muscle forces. J Neural Eng 3: 327-337, 2006. Cheng EJ, Brown IE, and Loeb GE. Virtual muscle: a computational approach to understanding the effects of muscle properties on motor control. Journal Of Neuroscience Methods 101: 117-130, 2000. Cheng EJ and Loeb GE. On the use of musculoskeletal models to interprete motor control strategies from performance data. J Neural Eng 5: 1-22, 2008. Churchland MM, Afshar A, and Shenoy KV. A central source of movement variability. Neuron 52: 1085-1096, 2006a. 133 Churchland MM, Santhanam G, and Shenoy KV. Preparatory activity in premotor and motor cortex reflects the speed of the upcoming reach. Journal Of Neurophysiology 96: 3130-3146, 2006b. Clamann HP. Statistical Analysis Of Motor Unit Firing Patterns In A Human Skeletal Muscle. Biophysical Journal 9: 1233-&, 1969. Cole KJ and Abbs JH. Coordination of three-joint digit movements for rapid finger-thumb grasp. J Neurophysiol 55: 1407-1423, 1986. Cole KJ and Abbs JH. Kinematic And Electromyographic Responses To Perturbation Of A Rapid Grasp. Journal Of Neurophysiology 57: 1498-1510, 1987. Collins JJ. The Redundant Nature Of Locomotor Optimization Laws. Journal Of Biomechanics 28: 251- 267, 1995. Cordo P, Gandevia SC, Hales JP, Burke D, and Laird G. Force and displacement-controlled tendon vibration in humans. Electroencephalogr Clin Neurophysiol 89: 45-53, 1993. Cordo PJ, Flores-Vieira C, Verschueren SMP, Inglis JT, and Gurfinkel V. Position sensitivity of human muscle spindles: Single afferent and population representations. Journal Of Neurophysiology 87: 1186-1195, 2002. Corning PA. Synergy and self-organization in the evolution of complex systems. Systems research 12: 89- 121, 1995. Crago PE, Houk JC, and Rymer WZ. Sampling Of Total Muscle Force By Tendon Organs. Journal of Neurophysiology (Bethesda) 47: 1069-1083, 1982. Crowninshield RD and Brand RA. A Physiologically Based Criterion Of Muscle Force Prediction In Locomotion. Journal Of Biomechanics 14: 793-801, 1981. Crowninshield RD, Johnston RC, Andrews JG, and Brand RA. Biomechanical Investigation Of Human Hip. Journal Of Biomechanics 11: 75-85, 1978. Danna-dos-Santos A, Slomka K, Zatsiorsky VM, and Latash ML. Muscle modes and synergies during voluntary body sway. Experimental Brain Research 179: 533-550, 2007. Davoodi R, Brown IE, Lan N, Mileusnic M, and Loeb GE. An integrated package of neuromusculo- skeletal modeling tools in SimulinkTM. 23rd IEEE/EMBS Ann. Int'l Conf., 2001. Davy DT and Audu ML. A Dynamic Optimization Technique For Predicting Muscle Forces In The Swing Phase Of Gait. Journal Of Biomechanics 20: 187-201, 1987. Delp SL and Loan JP. A Graphics-Based Software System To Develop And Analyze Models Of Musculoskeletal Structures. Computers In Biology And Medicine 25: 21-34, 1995. 134 Domkin D, Laczko J, Jaric S, Johansson H, and Latash ML. Structure of joint variability in bimanual pointing tasks. Experimental Brain Research 143: 11-23, 2002. Donders FC. Beitrag zur Lehre von den Bewegungen von den menschilthen Auges. Hollandeshen Beitragen zu den Anatomischen und Physiologischen Wissenschaften: 104-145, 1847. Faisal AA, Selen LPJ, and Wolpert DM. Noise in the nervous system. Nat Rev Neurosci 9: 292-303, 2008. Fitts PM. The Information Capacity Of The Human Motor System In Controlling The Amplitude Of Movement. Journal Of Experimental Psychology 47: 381-391, 1954. Flanders M, Pellegrini JJ, and Soechting JF. Spatial/Temporal Characteristics Of A Motor Pattern For Reaching. Journal Of Neurophysiology 71: 811-813, 1994. Flash T. The control of hand equilibrium trajectories in multijoint arm movements. Biological Cybernetics 57: 257-274, 1987. Flash T and Hogan N. The Coordination Of Arm Movements - An Experimentally Confirmed Mathematical-Model. Journal Of Neuroscience 5: 1688-1703, 1985. Flash T and Mussa-Ivaldi F. Human Arm Stiffness Characteristics During The Maintenance Of Posture. Experimental Brain Research 82: 315-326, 1990. Franklin DW, Liaw G, Milner TE, Osu R, Burdet E, and Kawato M. Endpoint stiffness of the arm is directionally tuned to instability in the environment. Journal Of Neuroscience 27: 7705-7716, 2007. Franklin DW, So U, Kawato M, and Milner TE. Impedance control balances stability with metabolically costly muscle activation. Journal Of Neurophysiology 92: 3097-3105, 2004. Fuglevand AJ, Winter DA, and Patla AE. Models Of Recruitment And Rate Coding Organization In Motor-Unit Pools. Journal Of Neurophysiology 70: 2470-2488, 1993. Galganski ME, Fuglevand AJ, and Enoka RM. Reduced Control Of Motor Output In A Human Hand Muscle Of Elderly Subjects During Submaximal Contractions. Journal Of Neurophysiology 69: 2108-2115, 1993. Garner BA and Pandy MG. Estimation of musculotendon properties in the human upper limb. Annals Of Biomedical Engineering 31: 207-220, 2003. Garner BA and Pandy MG. Musculoskeletal model of the upper limb based on the visible human male dataset. Comput Methods Biomech Biomed Engin 4: 93-126, 2001. Gelfand IM and Latash ML. On the problem of adequate language in biology. In: Progress in motor control, edited by Latash ML. Champaign, IL: Human Kinetics, 2002, p. 209-228. 135 Gelfand IM and Latash ML. On the problem of adequate language in movement science. Motor Control: 306-313, 1998. Ghez C, Gordon J, Ghilardi MF, Christakos CN, and Cooper SE. Roles Of Proprioceptive Input In The Programming Of Arm Trajectories. Cold Spring Harbor Symposia On Quantitative Biology 55: 837-847, 1990. Ghez C and Martin JH. The Control Of Rapid Limb Movement In The Cat.3. Agonist - Antagonist Coupling. Experimental Brain Research 45: 115-125, 1982. Gielen C, Vrijenhoek EJ, Flash T, and Neggers SFW. Arm position constraints during pointing and reaching in 3-D space. Journal Of Neurophysiology 78: 660-673, 1997. Gomez C, Canals J, Torres B, and Delgadogarcia JM. Analysis Of The Fluctuations In The Interspike Intervals Of Abducens Nucleus Neurons During Ocular Fixation In The Alert Cat. Brain Research 381: 401-404, 1986. Gomi H and Kawato M. Human arm stiffness and equilibrium-point trajectory during multi-joint movement. Biological Cybernetics 76: 163-171, 1997. Gomi H and Osu R. Task-dependent viscoelasticity of human multijoint arm and its spatial characteristics for interaction with environments. Journal Of Neuroscience 18: 8965-8978, 1998. Gonzalez RV, Hutchins EL, Barr RE, and Abraham LD. Development and evaluation of a musculoskeletal model of the elbow joint complex. Journal Of Biomechanical Engineering-Transactions Of The Asme 118: 32-40, 1996. Goodman SR and Latash ML. Feed-forward control of a redundant motor system. Biological Cybernetics 95: 271-280, 2006. Gordon J, Ghilardi MF, Cooper SE, and Ghez C. ACCURACY OF PLANAR REACHING MOVEMENTS.2. SYSTEMATIC EXTENT ERRORS RESULTING FROM INERTIAL ANISOTROPY. Experimental Brain Research 99: 112-130, 1994a. Gordon J, Ghilardi MF, and Ghez C. ACCURACY OF PLANAR REACHING MOVEMENTS.1. INDEPENDENCE OF DIRECTION AND EXTENT VARIABILITY. Experimental Brain Research 99: 97-111, 1994b. Gordon J, Ghilardi MF, and Ghez C. Impairments Of Reaching Movements In Patients Without Proprioception.1. Spatial Errors. Journal Of Neurophysiology 73: 347-360, 1995. Gribble PL, Mullin LI, Cothros N, and Mattar A. Role of cocontraction in arm movement accuracy. Journal Of Neurophysiology 89: 2396-2405, 2003. Gribble PL and Ostry DJ. Compensation for interaction torques during single- and multijoint limb movement. Journal Of Neurophysiology 82: 2310-2326, 1999. 136 Guigon E, Baraduc P, and Desmurget M. Computational motor control: Redundancy and invariance. Journal Of Neurophysiology 97: 331-347, 2007. Haggard P, Hutchinson K, and Stein J. Patterns of coordinated multi-joint movement. Experimental Brain Research 107: 254-266, 1995. Hamilton AFD, Jones KE, and Wolpert DM. The scaling of motor noise with muscle strength and motor unit number in humans. Experimental Brain Research 157: 417-430, 2004. Hamilton AFD and Wolpert DM. Controlling the statistics of action: Obstacle avoidance. Journal Of Neurophysiology 87: 2434-2440, 2002. Hannaford B and Stark L. Roles Of The Elements Of The Triphasic Control Signal. Exp Neurol 90: 619- 634, 1985. Harris CM and Wolpert DM. Signal-dependent noise determines motor planning. Nature 394: 780-784, 1998. Harwood MR, Mezey LE, and Harris CM. The spectral main sequence of human saccades. Journal Of Neuroscience 19: 9098-9106, 1999. Hatze H and Buys JD. Energy-Optimal Controls In Mammalian Neuromuscular System. Biological Cybernetics 27: 9-20, 1977. Haugstvedt JR, Berger RA, and Berglund LJ. A mechanical study of the moment-forces of the supinators and pronators of the forearm. Acta Orthopaedica Scandinavica 72: 629-634, 2001. Herzog W. Individual Muscle Force Estimations Using A Nonlinear Optimal-Design. Journal Of Neuroscience Methods 21: 167-179, 1987. Higuchi T, Imanaka K, and Hatayama T. Freezing degrees of freedom under stress: Kinematic evidence of constrained movement strategies. Human Movement Science 21: 831-846, 2002. Hoffer JA and Andreassen S. Regulation Of Soleus Muscle-Stiffness In Pre-Mammillary Cats - Intrinsic And Reflex Components. Journal Of Neurophysiology 45: 267-285, 1981. Hogan N. Mechanical impedance of single- and multi-articular systems. In: Multiple muscle systems: Biomechanics and movement organization, edited by Winters JM and Woo SL-Y: Springer, 1991, p. 149- 164. Hogan N. The Mechanics Of Multi-Joint Posture And Movement Control. Biological Cybernetics 52: 315- 332, 1985. Hogan N. An Organizing Principle For A Class Of Voluntary Movements. Journal Of Neuroscience 4: 2745-2754, 1984. 137 Hogan N. Planning And Execution Of Multijoint Movements. Canadian Journal Of Physiology And Pharmacology 66: 508-517, 1988. Holmes P, Full RJ, Koditschek D, and Guckenheimer J. The dynamics of legged locomotion: Models, analyses, and challenges. Siam Review 48: 207-304, 2006. Holzbaur KRS, Murray WM, and Delp SL. A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control. Annals Of Biomedical Engineering 33: 829- 840, 2005. Hore J, Watts S, and Tweed D. Arm Position Constraints When Throwing In 3 Dimensions. Journal Of Neurophysiology 72: 1171-1180, 1994. Hu XT, Jiang HH, Gu CL, Li CY, and Sparks DL. Reliability of oculomotor command signals carried by individual neurons. Proceedings Of The National Academy Of Sciences Of The United States Of America 104: 8137-8142, 2007. Hughlings Jackson J. On the comparative study of disease of the nervous system. British medical journal: 355-362, 1889. Hunter IW and Kearney RE. Dynamics Of Human Ankle Stiffness - Variation With Mean Ankle Torque. Journal Of Biomechanics 15: 747-752, 1982. Ijspeert AJ. A connectionist central pattern generator for the aquatic and terrestrial gaits of a simulated salamander. Biological Cybernetics 84: 331-348, 2001. Izawa J, Rane T, Donchin O, and Shadmehr R. Motor adaptation as a process of reoptimization. J Neurosci 28: 2883-2891, 2008. Johnson MA, Polgar J, Weightman D, and Appleton D. Data on the distribution of fibre types in thirty- six human muscles, An autopsy study. Journal of the Neurological Sciences 18: 111-129, 1973. Jones KE, Hamilton AFD, and Wolpert DM. Sources of signal-dependent noise during isometric force production. Journal Of Neurophysiology 88: 1533-1544, 2002. Karst GM and Hasan Z. Initiation Rules For Planar, 2-Joint Arm Movements - Agonist Selection For Movements Throughout The Work Space. Journal Of Neurophysiology 66: 1579-1593, 1991a. Karst GM and Hasan Z. Timing And Magnitude Of Electromyographic Activity For 2-Joint Arm Movements In Different Directions. Journal Of Neurophysiology 66: 1594-1604, 1991b. Kawato M, Maeda Y, Uno Y, and Suzuki R. Trajectory Formation Of Arm Movement By Cascade Neural Network Model Based On Minimum Torque-Change Criterion. Biological Cybernetics 62: 275-288, 1990. Kording KP and Wolpert DM. Bayesian decision theory in sensorimotor control. Trends Cogn Sci 10: 319-326, 2006. 138 Koshland GF, Galloway JC, and Nevoret-Bell CJ. Control of the wrist in three-joint arm movements to multiple directions in the horizontal plane. Journal Of Neurophysiology 83: 3188-3195, 2000. Kudina LP. Analysis of firing behaviour of human motoneurones within 'subprimary range'. Journal Of Physiology-Paris 93: 115-123, 1999. Kuechle DK, Newman SR, Itoi E, Morrey BF, and An KN. Shoulder muscle moment arms during horizontal flexion and elevation. Journal Of Shoulder And Elbow Surgery 6: 429-439, 1997. Lacquaniti F, Grasso R, and Zago M. Motor Patterns in Walking. News Physiol Sci 14: 168-174, 1999. Lacquaniti F, Terzuolo C, and Viviani P. The Law Relating The Kinematic And Figural Aspects Of Drawing Movements. Acta Psychologica 54: 115-130, 1983. Laidlaw DH, Bilodeau M, and Enoka RM. Steadiness is reduced and motor unit discharge is more variable in old adults. Muscle & Nerve 23: 600-612, 2000. Lan N. Analysis of an optimal control model of multi-joint arm movements. Biological Cybernetics 76: 107-117, 1997. Lan N. Stability analysis for postural control in a two-joint limb system. Ieee Transactions On Neural Systems And Rehabilitation Engineering 10: 249-259, 2002. Lan N and Baker L. Biomechanical Couplings Between Elbow And Forearm Movements. 26th IEEE/EMBS Ann. Intl. Conf., San Francisco, CA, 2004. Lan N and Crago PE. Optimal-Control Of Antagonistic Muscle-Stiffness During Voluntary Movements. Biological Cybernetics 71: 123-135, 1994. Lan N, Li Y, Sun Y, and Yang FS. Reflex regulation of antagonist muscles for control of joint equilibrium position. Ieee Transactions On Neural Systems And Rehabilitation Engineering 13: 60-71, 2005. Lan N and Murakata T. A realistic human elbow model for dynamic simulation. 25th Ann. Meeting ASB, San Diego, CA, 2001, p. 267-268. Lan N, Song D, and Gordon J. Systems engineering approach to computational sensorimotor control. 28th International Symposium of Computational Neuroscience, Montréal, Canada, 2006. Langenderfer J, Jerabek SA, Thangamani VI, Kuhn JE, and Hughes RE. Musculoskeletal parameters of muscles crossing the shoulder and elbow and the effect of sarcomere length sample size on estimation of optimal muscle length. Clinical Biomechanics 19: 664-670, 2004. Lashley KS. Integrative functions of the cerebral cortex. Physiological Reviews 13: 0001-0042, 1933. 139 Latash ML. On the evolution of the notion of synergy. In: Motor Control, Today and Tomorrow, edited by Gantchev GN, Mori S and Massion J. Sofia, 1999. Latash ML, Aruin AS, and Zatsiorsky VM. The basis of a simple synergy: reconstruction of joint equilibrium trajectories during unrestrained arm movements. Human Movement Science 18: 3-30, 1999. Latash ML, Scholz JP, and Schoner G. Toward a new theory of motor synergies. Motor Control 11: 276- 308, 2007. Latash ML, Shim JK, Smilga AV, and Zatsiorsky VM. A central back-coupling hypothesis on the organization of motor synergies: a physical metaphor and a neural model. Biological Cybernetics 92: 186- 191, 2005. Lemay MA and Crago PE. A dynamic model for simulating movements of the elbow, forearm, and wrist. Journal Of Biomechanics 29: 1319-1330, 1996. Lestienne F. Effects Of Inertial Load And Velocity On The Braking Process Of Voluntary Limb Movements. Experimental Brain Research 35: 407-418, 1979. Levin O, Wenderoth N, Steyvers M, and Swinnen SP. Directional invariance during loading-related modulations of muscle activity: evidence for motor equivalence. Experimental Brain Research 148: 62-76, 2003. Li ZM, Latash ML, and Zatsiorsky VM. Force sharing among fingers as a model of the redundancy problem. Experimental Brain Research 119: 276-286, 1998. Lieber RL. Skeletal muscle structure, function & plasticity: Lippincott Williams & Wilkins, 2002. Lin CCK and Crago PE. Structural model of the muscle spindle. Annals Of Biomedical Engineering 30: 68-83, 2002. Liu D and Todorov E. Evidence for the flexible sensorimotor strategies predicted by optimal feedback control. Journal Of Neuroscience 27: 9354-9368, 2007. Loeb GE. The control and response of mammalian muscle spindles during normally executed motor tasks. Exerc & Sports Sci Rev 12: 157-204, 1984. Loeb GE, Brown IE, and Cheng EJ. A hierarchical foundation for models of sensorimotor control. Experimental Brain Research 126: 1-18, 1999. Loeb GE and Marks WB. Optimal control principles for sensory transducers. Int. Sympos: The Muscle Spindle, edited by Boyd. MacMillan Ltd. London, 1985, p. 409-415. Marsden CD, Obeso JA, and Rothwell JC. The Function Of The Antagonist Muscle During Fast Limb Movements In Man. Journal Of Physiology-London 335: 1-13, 1983. 140 McAuley JH, Rothwell JC, and Marsden CD. Frequency peaks of tremor, muscle vibration and electromyographic activity at 10 Hz, 20 Hz and 40 Hz during human finger muscle contraction may reflect rhythmicities of central neural firing. Experimental Brain Research 114: 525-541, 1997. McIntyre J, Stratta F, and Lacquaniti F. Viewer-centered frame of reference for pointing to memorized targets in three-dimensional space. Journal Of Neurophysiology 78: 1601-1618, 1997. Meyer DE, Kornblum S, Abrams RA, Wright CE, and Smith JEK. Optimality In Human Motor- Performance - Ideal Control Of Rapid Aimed Movements. Psychological Review 95: 340-370, 1988. Mileusnic MP, Brown IE, Lan N, and Loeb GE. Mathematical models of proprioceptors: I. Control and transduction in the muscle spindle. journal of the Neurophysiol, 2006. Mileusnic MP and Loeb GE. Mathematical models of proprioceptors: II. Structure and function of the Golgi tendon organ. Journal of the Neurophysiol 96: 1789-1802, 2006. Morasso P. Spatial Control Of Arm Movements. Experimental Brain Research 42: 223-227, 1981. Moritz CT, Barry BK, Pascoe MA, and Enoka RM. Discharge rate variability influences the variation in force fluctuations across the working range of a hand muscle. Journal Of Neurophysiology 93: 2449-2459, 2005. Murray WM, Buchanan TS, and Delp SL. The isometric functional capacity of muscles that cross the elbow. Journal Of Biomechanics 33: 943-952, 2000. Murray WM, Delp SL, and Buchanan TS. Variation of Muscle Moment Arms with Elbow and Forearm Position. Journal of Biomechanics 28: 513-525, 1995. Mussa-Ivaldi FA, Hogan N, and Bizzi E. Neural Mechanical And Geometric Factors Subserving Arm Posture In Humans. Journal of Neuroscience 5: 2732-2743, 1985. Nagata Y, Miyamoto H, Nakano E, and Kawato M. Quantitative examinations of trajectory planning criteria in multijoint movement - TOPS model. Tech Rep Inst Electr Inform Commun Engineer NC2002-70: 25-35, 2002. Nakano E, Imamizu H, Osu R, Uno Y, Gomi H, Yoshioka T, and Kawato M. Quantitative examinations of internal representations for arm trajectory planning: Minimum commanded torque change model. Journal Of Neurophysiology 81: 2140-2155, 1999. Nelson WL. Physical Principles For Economies Of Skilled Movements. Biological Cybernetics 46: 135- 147, 1983. Nishikawa KC, Murray ST, and Flanders M. Do arm postures vary with the speed of reaching? Journal Of Neurophysiology 81: 2582-2586, 1999. Ostry DJ and Feldman AG. A critical evaluation of the force control hypothesis in motor control. Experimental Brain Research 153: 275-288, 2003. 141 Osu R, Burdet E, Franklin DW, Milner TE, and Kawato M. Different mechanisms involved in adaptation to stable and unstable dynamics. Journal Of Neurophysiology 90: 3255-3269, 2003. Osu R, Franklin DW, Kato H, Gomi H, Domen K, Yoshioka T, and Kawato M. Short- and long-term changes in joint co-contraction associated with motor learning as revealed from surface EMG. Journal Of Neurophysiology 88: 991-1004, 2002. Osu R and Gomi H. Multijoint muscle regulation mechanisms examined by measured human arm stiffness and EMG signals. Journal Of Neurophysiology 81: 1458-1468, 1999. Osu R, Kamimura N, Iwasaki H, Nakano E, Harris CM, Wada Y, and Kawato M. Optimal impedance control for task achievement in the presence of signal-dependent noise. Journal Of Neurophysiology 92: 1199-1215, 2004. Osyczka A. Multicriteria optimization for engineering design. In: Design Optimization, edited by Gero JS. New York, NY: Academic Press, Inc, 1985, p. 193-227. Pandy MG, Garner BA, and Anderson FC. Optimal-Control Of Non-Ballistic Muscular Movements - A Constraint-Based Performance Criterion For Rising From A Chair. Journal Of Biomechanical Engineering- Transactions Of The Asme 117: 15-26, 1995. Pandy MG and Zajac FE. Optimal Muscular Coordination Strategies For Jumping. Journal Of Biomechanics 24: 1-10, 1991. Pandy MG, Zajac FE, Sim E, and Levine WS. An Optimal-Control Model For Maximum-Height Human Jumping. Journal Of Biomechanics 23: 1185-1198, 1990. Pastor AM, Torres B, Delgadogarcia JM, and Baker R. Discharge Characteristics Of Medial Rectus And Abducens Motoneurons In The Goldfish. Journal Of Neurophysiology 66: 2125-2140, 1991. Pedotti A, Krishnan VV, and Stark L. Optimization Of Muscle-Force Sequencing In Human Locomotion. Mathematical Biosciences 38: 57-76, 1978. Perreault EJ, Kirsch RF, and Crago PE. Effects of voluntary force generation on the elastic components of endpoint stiffness. Experimental Brain Research 141: 312-323, 2001. Perreault EJ, Kirsch RF, and Crago PE. Voluntary control of static endpoint stiffness during force regulation tasks. Journal Of Neurophysiology 87: 2808-2816, 2002. Person RS and Kudina LP. Discharge Frequency And Discharge Pattern Of Human Motor Units During Voluntary Contraction Of Muscle. Electroencephalogr Clin Neurophysiol 32: 471-&, 1972. Piotrkiewicz M. An influence of afterhyperpolarization on the pattern of motoneuronal rhythmic activity. Journal Of Physiology-Paris 93: 125-133, 1999. Rancourt D and Hogan N. Stability in force-production tasks. Journal Of Motor Behavior 33: 193-204, 2001. 142 Rasmussen J, Damsgaard M, and Voigt M. Muscle recruitment by the min/max criterion - a comparative numerical study. Journal Of Biomechanics 34: 409-415, 2001. Robertson EM and Miall RC. Multi-joint limbs permit a flexible response to unpredictable events. Experimental Brain Research 117: 148-152, 1997. Sabes PN and Jordan MI. Obstacle avoidance and a perturbation sensitivity model for motor planning. Journal Of Neuroscience 17: 7119-7128, 1997. Sabes PN, Jordan MI, and Wolpert DM. The role of inertial sensitivity in motor planning. Journal Of Neuroscience 18: 5948-5957, 1998. Schmidt RA, Zelaznik H, Hawkins B, Frank JS, and Quinn JT. Motor-Output Variability - Theory For The Accuracy Of Rapid Motor Acts. Psychological Review 86: 415-451, 1979. Scholz JP, Kang N, Patterson D, and Latash ML. Uncontrolled manifold analysis of single trials during multi-finger force production by persons with and without Down syndrome. Experimental Brain Research 153: 45-58, 2003. Scholz JP and Schoner G. The uncontrolled manifold concept: identifying control variables for a functional task. Experimental Brain Research 126: 289-306, 1999. Scholz JP, Schoner G, and Latash ML. Identifying the control structure of multijoint coordination during pistol shooting. Experimental Brain Research 135: 382-404, 2000. Schweighofer N, Arbib MA, and Kawato M. Role of the cerebellum in reaching movements in humans. I. Distributed inverse dynamics control. Eur J Neurosci 10: 86-94, 1998. Sciavicco L and Sicilliano B. Modeling and control of robot manipulators: Springer. Scott SH and Loeb GE. The computation of position sense from spindles in mono- and multiarticular muscles. Journal Of Neuroscience 14: 7529-7540, 1994. Scott SH and Loeb GE. Mechanical properties of aponeurosis and tendon of the cat soleus muscle during whole-muscle isometric contractions. Journal of Morphology 224: 73-86, 1995. Seif-Naraghi AH and Winters JM. Optimized strategies for scaling goal-directed dynamic limb movements. In: Multiple Muscle Systems: Biomechanics and Movement Organization, edited by Winters JM and Woo SL-Y: Springer-Verlag, 1990, p. 312-334. Seireg A and Arvikar R. Biomechanical Analysis of the Musculoskeletal Structure for Medicine and Sports: Hemisphere Publishing Corp., 1989. Selen LPJ, Beek PJ, and van Dieen JH. Can co-activation reduce kinematic variability? A simulation study. Biological Cybernetics 93: 373-381, 2005. 143 Selen LPJ, van Dieen JH, and Beek PJ. Impedance modulation and feedback corrections in tracking targets of variable size and frequency. Journal Of Neurophysiology 96: 2750-2759, 2006. Semmler JG, Kornatz KW, Kern DS, and Enoka RM. Motor unit synchronization reduces the steadiness of anisometric contractions by a hand muscle. In: Soc Neurosci Abstr 27, 2001, p. 3. Semmler JG and Nordstrom MA. Motor unit discharge and force tremor in skill- and strength-trained individuals. Experimental Brain Research 119: 27-38, 1998. Shadmehr R and Arbib MA. A Mathematical-Analysis Of The Force-Stiffness Characteristics Of Muscles In Control Of A Single Joint System. Biological Cybernetics 66: 463-477, 1992. Shadmehr R and Mussaivaldi FA. Adaptive Representation Of Dynamics During Learning Of A Motor Task. Journal Of Neuroscience 14: 3208-3224, 1994. Sherrington CS. Flexion reflex of the limb, crossed extension reflex, and reflex stepping and standing. Journal of physiology: 28-121, 1910. Singh K, Richmond FJR, and Loeb GE. Recruitment properties of intramuscular and nerve-trunk stimulating electrodes. IEEE Trans Rehabil Eng 8: 276-285, 2000. Slifkin A and Newell KM. Noise, information transmission, and force variability. Journal Of Experimental Psychology-Human Perception And Performance 25: 837-851, 1999. Smeets JBJ and Brenner E. A new view on grasping. Motor Control 3: 237-271, 1999. Song D, Lan N, and Gordon J. Biomechanical constraints on equilibrium point control of the multi-joint arm, a simulation study. ASB 2007, Stanford Univ., CA, 2007a. Song D, Lan N, and Gordon J. Simulated hand variability during multi-joint arm posture control. Neuroscience 2006, Georgia World Congress Center, Atlanta, GA, USA, 2006. Song D, Lan N, Loeb GE, and Gordon J. Model-Based Sensorimotor Integration for Multi-Joint Control - Development of a Virtual Arm Model. Annals of Biomedical Engineering 36: 1033-1048, 2008a. Song D, Mileusnic M, Lan N, and Gordon J. A Sensorimotor Systems Model for Dynamic Simulation of Arm Movement Control. 16th Annual Neural Control of Movement Meeting, Key biscayn Florida, USA, 2006 (a). Song D, Raphael G, Lan N, and Loeb GE. Improvement in Computational Efficiency of Virtual Muscle Model. BMES 2007, Los Angeles, USA, 2007b. Song D, Raphael G, Montazem PT, Lan N, and Loeb GE. Computationally efficient models of neuromuscular recruitment and mechanics. Journal of Neural Engineering 5: 175-184, 2008b. 144 Sosnoff JJ, Valantine AD, and Newell KM. Independence between the amount and structure of variability at low force levels. Neuroscience Letters 392: 165-169, 2006. Stein RB, Weber DJ, Aoyagi Y, Prochazka A, Wagenaar JBM, Shoham S, and Normann RA. Coding of position by simultaneously recorded sensory neurones in the cat dorsal root ganglion. Journal Of Physiology-London 560: 883-896, 2004. Taylor AM, Christou EA, and Enoka RM. Multiple features of motor-unit activity influence force fluctuations during isometric contractions. Journal Of Neurophysiology 90: 1350-1361, 2003. Thoroughman KA and Feller KJ. Gravitational effects on torque change and variance optimization in reaching movements. SFN, Washington, DC, 2003. Thoroughman KA and Shadmehr R. Electromyographic Correlates of Learning an Internal Model of Reaching Movements. J Neurosci 19: 8573-8588, 1999. Tin C and Poon C-S. Internal models in sensorimotor integration: perspectives from adaptive control theory. Journal of Neural Engineering 2: 147-163, 2005. Todorov E. Cosine tuning minimizes motor errors. Neural Comput 14: 1233-1260, 2002. Todorov E. Optimality principles in sensorimotor control. Nature Neuroscience 7: 907-915, 2004. Todorov E. Stochastic optimal control and estimation methods adapted to the noise characteristics of the sensorimotor system. Neural Comput 17: 1084-1108, 2005. Todorov E and Ghahramani Z. Analysis of the synergies underlying complex hand manipulation. 26th Annual International Conference of the IEEE EMBS, San Francisco, CA, USA, 2004, p. 391-398. Todorov E and Jordan MI. Optimal feedback control as a theory of motor coordination. Nature Neuroscience 5: 1226-1235, 2002. Tolhurst DJ, Movshon JA, and Dean AF. The Statistical Reliability Of Signals In Single Neurons In Cat And Monkey Visual-Cortex. Vision Research 23: 775-785, 1983. Tseng YW, Scholz JP, and Schoner G. Goal-equivalent joint coordination in pointing: Affect of vision and arm dominance. Motor Control 6: 183-207, 2002. Uno Y, Kawato M, and Suzuki R. Formation And Control Of Optimal Trajectory In Human Multijoint Arm Movement - Minimum Torque-Change Model. Biological Cybernetics 61: 89-101, 1989. Valero-Cuevas FJ. Predictive modulation of muscle coordination pattern magnitude scales fingertip force magnitude over the voluntary range. Journal Of Neurophysiology 83: 1469-1479, 2000. Valero-Cuevas FJ, Venkadesan M, and Todorov E. Structured variability of muscle activations supports the minimal intervention principle of motor control. Journal of Neurophysiology, submitted. 145 Valero-Cuevas FJ, Zajac FE, and Burgar CG. Large index-fingertip forces are produced by subject- independent patterns of muscle excitation. Journal Of Biomechanics 31: 693-703, 1998. van Beers RJ. The sources of variability in saccadic eye movements. Journal Of Neuroscience 27: 8757- 8770, 2007. van Beers RJ, Haggard P, and Wolpert DM. The role of execution noise in movement variability. Journal Of Neurophysiology 91: 1050-1063, 2004. van den Dobbelsteen JJ, Brenner E, and Smeets JBJ. Endpoints of arm movements to visual targets. Experimental Brain Research 138: 279-287, 2001. Van Galen GP and De Jong WP. Fitts Law As The Outcome Of A Dynamic Noise Filtering Model Of Motor Control. Human Movement Science 14: 539-571, 1995. Vangalen GP and Dejong WP. Fitts Law As The Outcome Of A Dynamic Noise Filtering Model Of Motor Control. Human Movement Science 14: 539-571, 1995. Veeger HEJ, Yu B, An KN, and Rozendal RH. Parameters for modeling the upper extremity. Journal Of Biomechanics 30: 647-652, 1997. Vereijken B, Vanemmerik REA, Whiting HTA, and Newell KM. Free(Z)Ing Degrees Of Freedom In Skill Acquisition. Journal Of Motor Behavior 24: 133-142, 1992. Vindras P and Viviani P. Frames of reference and control parameters in visuomanual pointing. Journal Of Experimental Psychology-Human Perception And Performance 24: 569-591, 1998. Viviani P and Schneider R. A Developmental-Study Of The Relationship Between Geometry And Kinematics In Drawing Movements. Journal Of Experimental Psychology-Human Perception And Performance 17: 198-218, 1991. Winter DA. In: Perspectives on the coordination of movement, edited by Wallace SA. Amsterdam: Elsevier, 1989, p. 329-363. Wolpert DM, Ghahramani Z, and Jordan MI. An Internal Model For Sensorimotor Integration. Science 269: 1880-1882, 1995. Wood JE, Meek SG, and Jacobsen SC. Quantitation Of Human Shoulder Anatomy For Prosthetic Arm Control.1. Surface Modeling. Journal Of Biomechanics 22: 273-&, 1989. Wright CE. In: Attention and Performance XIII: Motor representation and control, edited by Jeannerod M. Hillsdale, New Jersey: Lawrence Erlbaum, 1990, p. 294-320. Wu G, van der Helm FCT, Veeger HEJ, Makhsous M, Van Roy P, Anglin C, Nagels J, Karduna AR, McQuade K, Wang XG, Werner FW, and Buchholz B. ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion - Part II: shoulder, elbow, wrist and hand. Journal of Biomechanics 38: 981-992, 2005. 146 Yakovenko S, Gritsenko V, and Prochazka A. Contribution of stretch reflexes to locomotor control: a modeling study. Biological Cybernetics 90: 146-155, 2004. Zajac FE. Muscle And Tendon - Properties, Models, Scaling, And Application To Biomechanics And Motor Control. Critical Reviews In Biomedical Engineering 17: 359-411, 1989. Zatsiorsky VM, Li ZM, and Latash ML. Enslaving effects in multi-finger force production. Experimental Brain Research 131: 187-195, 2000. Zhu G, Rotea MA, and Skelton R. A convergent algorithm for the output covariance constraint control problem. SIAM J Control Optim 35: 341-361, 1997.

Abstract (if available)

Abstract Precise control of posture is necessary for the success of any motor actions. The strategies hereby precise movement is achieved reflect many complexly interacting factors, including the tendency of active muscles to generate substantial motor noise as well as to influence the mechanical impedance of a limb with multiple degrees of freedom. A realistic model of the human sensorimotor system would be instrumental in understanding the underlying principles whereby the CNS controls the kinematic effects of signal dependent noise (SDN) in the actuators of a redundant musculoskeletal system. In this work we first constructed a comprehensive sensorimotor system model of the human arm that accurately represents the physiological properties of human muscles and proprioceptors. A simplified two-joint six-muscle arm model was then used for the simulation study of hand postural control in the face of SDN. In single-joint systems, co-contraction has been shown to be effective to reduce the kinematic effect of SDN. But the strategy is energetically inefficient, and extension to a multi-joint system is not obvious because of the intersegmental dynamics and redundant musculature.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Spinal-like regulator for control of multiple degree-of-freedom limbs

PDF

Neuromuscular dynamics in the context of motor redundancy

PDF

Investigating the role of muscle physiology and spinal circuitry in sensorimotor control

PDF

Experimental and model-based analyses of physiological determinants of force variability

PDF

Understanding reactive balance control strategies in non-disabled and post-stroke gait

PDF

Inference of computational models of tendon networks via sparse experimentation

PDF

insideOut: Estimating joint angles in tendon-driven robots using Artificial Neural Networks and non-collocated sensors

PDF

The representation, learning, and control of dexterous motor skills in humans and humanoid robots

PDF

A non-invasive and intuitive command source for upper limb prostheses

PDF

Using nonlinear feedback control to model human landing mechanics

PDF

A medical imaging informatics based human performance analytics system

PDF

Regulation of linear and angular impulse generation: implications for athletic performance

PDF

On the electrophysiology of multielectrode recordings of the basal ganglia and thalamus to improve DBS therapy for children with secondary dystonia

PDF

Engineering scalable two- and three-dimensional striated muscle microtissues for human disease modeling

PDF

Identifying injury risk, improving performance, and facilitating learning using an integrated biomechanics informatics system (IBIS)

PDF

Upper extremity control and dynamics during manual wheelchair propulsion at different speeds and wheelchair configurations

PDF

Sensory acquisition for emergent body representations in neuro-robotic systems

PDF

A percutaneously implantable wireless neurostimulator for treatment of stress urinary incontinence

PDF

Investigation of preclinical testing methods for total ankle replacements

PDF

Economic model predictive control for building energy systems

Asset Metadata

Creator Song, Dan (author)

Core Title Model-based studies of control strategies for noisy, redundant musculoskeletal systems

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Publication Date 07/24/2008

Defense Date 06/05/2008

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag motor control,neural noise,OAI-PMH Harvest

Language English

Advisor Gerald E. Loeb (committee chair), Gordon, James (committee member), Lan, Ning (committee member), McNitt-Gray, Jill L. (committee member), Valero-Cuevas, Francisco J. (committee member)

Creator Email dansong@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m1410

Unique identifier UC1224170

Identifier etd-SongD-20080724 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-200088 (legacy record id),usctheses-m1410 (legacy record id)

Legacy Identifier etd-SongD-20080724.pdf

Dmrecord 200088

Document Type Dissertation

Rights Song, Dan

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu