Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Computational modeling of neuromusculoskeletal systems : from filaments to behavior Yeo, Sang Hoon 2012

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2013_spring_yeo_sanghoon.pdf [ 4.65MB ]
JSON: 24-1.0052217.json
JSON-LD: 24-1.0052217-ld.json
RDF/XML (Pretty): 24-1.0052217-rdf.xml
RDF/JSON: 24-1.0052217-rdf.json
Turtle: 24-1.0052217-turtle.txt
N-Triples: 24-1.0052217-rdf-ntriples.txt
Original Record: 24-1.0052217-source.json
Full Text

Full Text

Computational Modeling of Neuromusculoskeletal Systems: from Filaments to Behavior  by Sang Hoon Yeo B.Sc., Seoul National University, 1998 M.Sc., Seoul National University, 2001  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Computer Science)  The University of British Columbia (Vancouver) December 2012 c Sang Hoon Yeo, 2012  Abstract This thesis describes computational approaches to modeling and simulating aspects of the neuromusculoskeletal system. We make contributions to models at three different levels of detail. We first investigate the mechanics of shortening muscle and evaluate two forms of the traditional Hill-type muscle model, force scaling and f-max scaling, and show that the f-max scaling model is significantly better at predicting experimental results. We hypothesize a new model called the winding filament model that incorporates the role of titin during active force development. Based on the proposed hypothesis, we develop a computational model that is able to simulate residual force enhancement. The suggested model can qualitatively simulate the pattern of the force enhancement observed in previous studies. In order to model the higher levels of the system consisting of muscles and bones, we propose an optimal design framework for estimating parameters of the musculoskeletal model. The method finds a set of morphological and physiological parameters that can optimally simulate the measured force and moment at the point of action. We apply the suggested framework to modeling two rat hindlimb muscles, gracilis posticus and posterior part of biceps femoris, to see if the traditional  ii  line segment based muscle geometry model is valid for musculoskeletal system modeling. The result shows that even a complex muscle like biceps femoris can be well modeled as a line segment, but its estimated insertion point is far from that of the traditional model based on anatomy. Finally, this thesis addresses a behavioral aspect of biological movement; in particular, how a high level movement is planned and controlled, in coordination with perception. We present a fully generative model of object interception that can simulate realistic, human-like behavior of ball catching for given arbitrary ball trajectory. The model includes a simplified probabilistic model of vision, a model of eye movements combining saccades and pursuit, and corresponding head, hand and body movements. The movements are constructed from submovements. By combining these components, realistic interception behavior is simulated with minimal user intervention.  iii  Preface The works comprising this thesis are taken from papers that have been published or are planned to be submitted in the near future in collaboration with other people. The work described in Section 2.2 is from a paper currently being prepared for submission. The present author worked as the first author and the main writer of the paper. The text was later edited by Dr. Dinesh K. Pai and Dr. Kiisa C. Nishikawa who provided overall guidance and discussion on data collection and the analysis. Other coauthors, Dr. Jenna A. Monroy and Dr. A. Kristopher Lappin, worked on data collection on mouse and toad respectively. The winding filament hypothesis presented in Section 2.3 is from (K. C. Nishikawa, J. A. Monroy, T. E. Uyeno, S. -H. Yeo, D. K. Pai, and S. L. Lindstedt. Is titin a ‘winding filament’? A new twist on muscle contraction. Proc. R. Soc. B, 279(1730):981–990, 2012). The hypothesis was initiated and developed by Dr. Nishikawa, while coauthors including the present author participated in the discussion. In addition, the present author developed the mechanical model of titin winding and numerical simulation of the residual force enhancement in collaboration with Dr. Nishikawa and Dr. Pai. Except for some figures and legends cited from the paper, the text in this chapter is written by the present author based on the  iv  published paper and does not appear in any other manuscript. Musculoskeletal modeling work in Chapter 3 is based on (S. -H. Yeo, C. H. Mullens, T. G. Sandercock, D. K. Pai, and M. C. Tresch. Estimation of musculoskeletal models from in situ measurements of muscle action in the rat hindlimb. J. Exp. Biol., 213(5):735–746, 2011). The present author was the primary author of the work and performed majority of the data analysis. Christopher Mullens and Dr. Matthew C. Tresch performed data collection. Dr. Thomas G. Sandercock participated in data analysis and specifically worked on the analysis comparing the result with that of [88]. Dr. Pai participated in the early development of the idea and provided discussions on the data collection and analysis. The entire project was led by Dr. Tresch. The text and figures are based largely on the paper, which was written by the present author and Dr. Tresch. Chapter 4 is from (S. -H. Yeo, M. Lesmana, D. R. Neog, and D. K. Pai. Eyecatch: Simulating Visuomotor Coordination for Object Interception. ACM. Trans. Graph. (Proc. SIGGRAPH), 31(4), 2012). The text is from the paper, which was written by the present author and edited by Dr. Pai. The present author was the primary author who worked on the majority of the data processing and analysis and development of the simulation model. Martin Lesmana worked on the data collection, inverse kinematics module of the simulation model, character building and rendering. Debanga R. Neog worked on the data collection, hand movement analysis, character building and rendering. Dr. Pai provided overall guidance of the project and participated in writing.  v  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xiv  1  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  2  Modeling Muscle Mechanics . . . . . . . . . . . . . . . . . . . . . .  10  2.1  Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  10  2.1.1  11  2.2  Anatomy and Physiology of Skeletal Muscle . . . . . . .  Simulating the Dynamics of Shortening Muscle Using Hill-type Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  14  2.2.1  Models of FL and FV Relationships . . . . . . . . . . . .  15  2.2.2  Two forms of FLV relationship . . . . . . . . . . . . . . .  17  2.2.3  Data Collection . . . . . . . . . . . . . . . . . . . . . . .  20  vi  2.2.4  tionships . . . . . . . . . . . . . . . . . . . . . . . . . .  22  2.2.5  Results . . . . . . . . . . . . . . . . . . . . . . . . . . .  27  2.2.6  Discussion . . . . . . . . . . . . . . . . . . . . . . . . .  32  Winding Filament Hypothesis . . . . . . . . . . . . . . . . . . .  34  2.3.1  Currently Unsolved Problems in Muscle Mechanics . . . .  35  2.3.2  Winding Filament: a New Hypothesis . . . . . . . . . . .  40  2.3.3  Computational Model of the Winding Filament Hypothesis  48  2.3.4  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . .  58  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  62  Modeling Musculoskeletal Systems . . . . . . . . . . . . . . . . . . .  64  3.1  Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  64  3.2  Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  67  3.2.1  Experimental Preparation . . . . . . . . . . . . . . . . . .  67  3.2.2  Experimental Protocol and Data Collection . . . . . . . .  70  3.2.3  Data Analysis and Modelling . . . . . . . . . . . . . . . .  72  Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  79  2.3  2.4 3  Data Processing and Estimation of FV, FLV, and FL Rela-  3.3  3.3.1  Forces and Moments Evoked by In Situ Stimulation of Rat Hindlimb Muscles . . . . . . . . . . . . . . . . . . . . .  3.3.2  Force and Moment Fields Produced by Stimulation of Gracilis Posticus and Biceps Femoris . . . . . . . . . . . . .  3.3.3  79  81  Using Measured Force and Moment Fields to Create Muscle Models . . . . . . . . . . . . . . . . . . . . . . . . .  vii  85  3.3.4  3.4  Hindlimb . . . . . . . . . . . . . . . . . . . . . . . . . .  90  Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  91  3.4.1  Characterization of Muscle Actions by In Situ Stimulation  91  3.4.2  Developing Musculoskeletal Muscle Models Based on In Situ Measurements of Muscle Actions . . . . . . . . . . .  94  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . .  97  Modeling Sensorimotor Behavior . . . . . . . . . . . . . . . . . . . .  99  4.1  Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  99  4.2  Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100  4.3  Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103  4.4  Visual Estimation of Target Movement . . . . . . . . . . . . . . . 106  3.4.3 4  4.5  5  Comparison to Existing Musculoskeletal Models of the Rat  4.4.1  Vision . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107  4.4.2  Internal Model and Bayesian State Estimation . . . . . . . 108  Movement Generation . . . . . . . . . . . . . . . . . . . . . . . 111 4.5.1  Gaze Movements . . . . . . . . . . . . . . . . . . . . . . 113  4.5.2  Head Movements . . . . . . . . . . . . . . . . . . . . . . 116  4.5.3  Hand Movements for Interception . . . . . . . . . . . . . 118  4.5.4  Body Movements . . . . . . . . . . . . . . . . . . . . . . 122  4.6  Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124  4.7  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.1  Overview of Contributions . . . . . . . . . . . . . . . . . . . . . 129  5.2  Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 viii  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134  ix  List of Tables Table 2.1  Cross validation result of the FLV models for all data . . . . .  32  Table 2.2  Geometric parameters of sarcomere used in the modeling . . .  48  Table 4.1  Parameters for error standard deviations of foveal vision . . . . 107  x  List of Figures Figure 1.1  Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . .  2  Figure 2.1  Inner structure of myofibril and sarcomere . . . . . . . . . . .  13  Figure 2.2  The force-length relationship of sarcomeres in frog muscle . .  16  Figure 2.3  The force-length and force-velocity relationship of the skeletal muscle . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  16  Figure 2.4  Experimental setup for mouse soleus . . . . . . . . . . . . . .  21  Figure 2.5  Load-clamp data from mouse soleus muscle . . . . . . . . . .  23  Figure 2.6  Illustration of the load-clamp experiment from the viewpoint of the Hill-type muscle model . . . . . . . . . . . . . . . . .  Figure 2.7  Force-velocity relationship of four mouse soleus muscles estimated from the measured initial velocity . . . . . . . . . . . .  Figure 2.8  27  Estimated slow phase responses of four mouse soleus muscles using force scaling model . . . . . . . . . . . . . . . . . . .  Figure 2.9  24  28  Estimated slow phase responses of four mouse soleus muscles using f-max scaling model . . . . . . . . . . . . . . . . . . .  29  Figure 2.10 Estimated force-length curves for four mouse soleus muscles using force scaling and f-max scaling model . . . . . . . . . .  xi  30  Figure 2.11 Estimated slow phase responses and corresponding FL curves of the toad DM muscle . . . . . . . . . . . . . . . . . . . . .  31  Figure 2.12 Schematic of a skeletal muscle half-sarcomere illustrating the layout of titin . . . . . . . . . . . . . . . . . . . . . . . . . .  42  Figure 2.13 Schematic illustrating the winding filament hypothesis . . . .  46  Figure 2.14 Kinematics of titin winding . . . . . . . . . . . . . . . . . . .  50  Figure 2.15 Elasticity of PEVK in single titin . . . . . . . . . . . . . . . .  53  Figure 2.16 Simulation of winding in isometric contraction for different lengths . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  57  Figure 2.17 Simulation of winding in isometric contraction followed by stretch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  59  Figure 2.18 Residual force enhancement of single fibers from the frog tibialis anterior muscle . . . . . . . . . . . . . . . . . . . . . .  60  Figure 2.19 Simulation of activation-lengthening of the rabbit psoas muscle by winding filament hypothesis . . . . . . . . . . . . . . .  60  Figure 2.20 Simulated residual force enhancement . . . . . . . . . . . . .  61  Figure 3.1  Experimental setup and kinematic structure for the hindlimb .  68  Figure 3.2  Force and moments evoked by in situ muscle stimulation in the rat hindlimb . . . . . . . . . . . . . . . . . . . . . . . . . . .  Figure 3.3  80  Force and moment fields evoked by in situ stimulation of GRp for each animal . . . . . . . . . . . . . . . . . . . . . . . . .  82  Figure 3.4  Force and moment fields evoked by in situ stimulation of BFp  83  Figure 3.5  Comparison of force and moment fields measured in “restrained” and “released” trials . . . . . . . . . . . . . . . . . . . . . .  xii  84  Figure 3.6  Example prediction of the force and moment fields evoked by GRp stimulation . . . . . . . . . . . . . . . . . . . . . . . .  Figure 3.7  86  Example prediction of force and moment fields evoked by BFp stimulation . . . . . . . . . . . . . . . . . . . . . . . . . . .  86  Figure 3.8  Estimated parameters for GRp and BFp models for each animal 87  Figure 3.9  Prediction of force and moment fields of GRp and BFp using anatomical parameters . . . . . . . . . . . . . . . . . . . . .  90  Figure 3.10 Tendon structure in BFp . . . . . . . . . . . . . . . . . . . .  93  Figure 4.1  Captured catching behavior . . . . . . . . . . . . . . . . . . . 103  Figure 4.2  Eye movements for three different ball trajectories . . . . . . 104  Figure 4.3  Tangential velocity of captured hand and eye movements . . . 105  Figure 4.4  Eye coordinate frame and the corresponding uncertainty of an object . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106  Figure 4.5  Decomposed submovement of hand trajectory . . . . . . . . . 114  Figure 4.6  Overall control structure of gaze . . . . . . . . . . . . . . . . 115  Figure 4.7  Distance from head to the interception point . . . . . . . . . . 119  Figure 4.8  Decision variables and ball apex . . . . . . . . . . . . . . . . 120  Figure 4.9  Transition of the hand trajectory around the apex . . . . . . . 121  Figure 4.10 Submovement composition of the hand trajectory for a given ball trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Figure 4.11 Simulation of three dimensional ball catching . . . . . . . . . 125 Figure 4.12 Simulated gaze behavior for a given ball trajectory . . . . . . 126 Figure 4.13 Comparison between real and simulated hand trajectories for the same ball catching tasks . . . . . . . . . . . . . . . . . . 127  xiii  Acknowledgments I would like to thank my supervisor, Dinesh Pai, for his guidance and mentoring during my PhD. His broad and deep mind has been a huge inspiration. I would also like to thank the members of the supervisory and examination committee, Professors Elizabeth Croft, John Gosline, Michiel van de Panne, and Matthew Tresch for their helpful comments. I also would like to thank to my mentors and collaborators, Professors Kiisa Nishikawa, Thomas Sandercock, Miriam Spering, Ted Uyeno and Jenna Monroy for their help. I am very grateful to all the past and current members of the sensorimotor systems lab. Tim, Misha, Ye, Francois, Benjamin, Ben, Atsutoshi, Garrett, Andrew, Danny, Martin, Duo, Josh, Mahk, Dave, Guillaume, Debanga, Cole, Shin, Mitsunori, Bin, Qi, thank you for all your help and friendship. This work would have been impossible without love and support from my family. I thank my family members: Bongseo, Jiyoung, Jinwon, Yoonhui, Soowan, Jaewoo, Taewhan Hwang, Kumsil Won. I would like to address special thanks to my parents, Neungku Yeo, Sookyoung Kong for their endless support. I also would like to remember my grandmom, Guiok Kwon, who would have been proud of her grandson. Lastly, I thank my wife, Eun-Seung for always being there for me.  xiv  Chapter 1  Introduction Biological movement is a product of a complex mechanical and neural interaction among muscles, the skeleton and the brain. The term “neuromusculoskeletal” has been coined to refer to the entire system that is responsible for the generation of biological movement. This thesis addresses the questions that arise in building integrated, generative models of the neuromusculoskeletal system, and it studies modeling of three specific aspects of the system.  Thesis Outline From an engineering viewpoint, the neuromusculoskeletal system is a large scale dynamical system with tremendous complexity as shown in Figure 1.1. At the lowest level, movement originates from molecular level interactions among myofilaments. These nano-scale interactions are integrated into higher level components; sarcomeres, myofibrils, muscle fibers, and the whole muscle, each of which contains a large amount of architectural complexity and degrees of freedom. At the highest level, multiple muscles non-linearly interact with the skeleton, tendons and 1  target of simulation  model to be investigated  sensorimotor behavior neuromusculoskeletal system  Three dimensional visuomotor coordination pattern  Shared motor program in eye-hand coordination  Chapter 4  musculosksletal system  Six dimensional force and moment field  Point-to-point geometric representation of muscle  Chapter 3  Isotonic shortening profile  Hill-type muscle model  Chapter 2.2  Residual force enhancement  Winding filament model  Chapter 2.3  modeling layer  muscle  muscle fiber, fiber architecture  sarcomere, myofibril  myofilaments  Figure 1.1: Thesis overview ligaments to produce body movement. Moreover, this whole structure is monitored and controlled by neural signals that contain a similar amount of complexity. Due to this highly complex and hierarchical nature, building a model of the neuromusculoskeletal system is a challenging task. To model the system in an understandable form, we need to substantially reduce the complexity of each layer. However, we should also avoid oversimplification which could lead to serious misunderstanding of the real working principles of the system. For this reason, we need an appropriate strategy to model each layer of the neuromusculoskeletal system while maintaining the overall integrity. In this thesis, we implement two main frameworks to model the neuromusculoskeletal system: We build generative, computational models that can simulate the mechanical behavior of the system. At  2  the same time, we attempt to maintain an integrative view that covers the overall hierarchy of the system.  Computational Modeling and Simulation Since building a realistic model is key to understanding a system, we build a model of each component and layer them as reliably as possible. Therefore, throughout all the modeling studies presented in this thesis, we specifically aim to find generative, computational models that can simulate the experimental observations in a sufficiently wide range of time and space, and we define our methodology as a computational modeling approach. In contrast to descriptive models that aim to explain the data by logical reasoning and statistical inference, we use generative tools of numerical analysis, robotics, and control theory to build models to simulate the behavior of the neuromusculoskeletal system. The computational modeling approach gives us a systematic way to understand the working principles of the system. It is important to note that the major benefit of generative modeling comes from the interactive relationship between modeling and understanding. We build a model based on our understanding, but the generative modeling procedure can also provide insight into understanding the real design principles of the system. As a convincing example, it is well known that one of the earliest mechanical understandings of the bird wing was achieved by the pioneers of aviation such as Otto Liliental [113] and the Wright brothers [174] as they were building actual machines that fly. Building a working model is a direct way to identify challenges that the system faces and to discover the feasible solutions that can be implemented in the system. Another benefit of the approach comes from testability and generalization. 3  Since the model is able to produce simulation data that can be quantitatively compared to real data, testability is immediately achieved; anyone can conveniently verify or falsify the proposed hypothesis that the suggested model is based on. Since the simulation does not explicitly rely on specific data, its ability can be generalized beyond the scope of the data that the model is built on. This gives us a quantitative prediction of a wide range of system behavior and provides a useful guideline for future experiments. Together with the methodological advantages, computational models are also useful in practical science and engineering. Since biological systems have been improved through a long history of adaptation, some of their mechanical and control abilities easily outperform those of man-made machines. For this reason, incorporating the discovered design and control principles of the neuromusculoskeletal system into an engineered system can be a promising way to improve performance, and computational modeling provides a direct pathway to implementation. There have been successful studies on developing bio-mimetic robot systems [27, 97, 139] and bio-inspired controllers for character animation [178, 187, 199]. In clinical studies, computational simulation can replace expensive data collection and ad hoc reasoning in rehabilitation or surgical planning. Pioneering studies have shown that neuromusculoskeletal simulation can be efficiently used in the analysis and surgical planning of motor diseases of various body parts, including the eye [155, 190], finger [166], and lower extremity [4, 32].  Maintaining an Integrative View From an engineering perspective, the layers and components of the neuromusculoskeletal system can be easily mapped to those of a robot system and therefore 4  can be individually modeled and analyzed in a similar way. This is often referred to as a reductionist approach [16]. Reductionism provides a feasible and efficient way to construct complex models of the neuromusculoskeletal system. Like all other studies of this kind, this thesis also develops models heavily based on the reductionist approach. However, taking this highly engineering-oriented view can seriously underestimate the magnitude of the problem we are facing in neuromusculoskeletal modeling. The main difference between robot modeling and neuromusculoskeletal system modeling is the degree of consensus, i.e., the maturity of the field. As will be discussed in the following chapters, there is indeed no clear standard model that has reached consensus, and no model in any layer of the neuromusculoskeletal system, from single sarcomere contraction to whole body movement, that has been mechanically verified and reasonably guaranteed for use in modeling studies at higher layers. As will be discussed in the following chapters, model in any layer of the neuromusculoskeletal system, from single sarcomere contraction to whole body movement, has not reached consensus. Various models that are used in current studies have not been mechanically verified and are not reasonably guaranteed for use in modeling studies at higher layers. For this reason, it would be somewhat dangerous to regard the neuromusculoskeletal system as a layered collection of well-defined mechanical components. Therefore, modeling in any layer of the system cannot completely rely on known models in lower levels without re-investigation. Unfortunately, however, it is not difficult to find in various modeling studies excessive usage of so-called “standard” and “traditional” models without serious consideration of their mechanical validity. In musculoskeletal modeling studies, 5  the Hill-type muscle model and cross-bridge theory are widely accepted as a standard models of muscle mechanics. However, even the original pioneers of these models reported numerous results that contradict the predictions of their models. With regard to neural control, as pointed out by Buchanan [18], neuroscientists often find it surprising that biomechanics researchers sometimes replace the entire nervous system with a simple, unverified equation. For this reason, this thesis implements an integrative view of modeling the mechanics of the neuromusculoskeletal system, spanning the range from muscle mechanics to sensorimotor behavior. In various layers of the system hierarchy, we focus on whether the common principles or methodologies used in modeling are indeed able to predict the mechanical behavior of the system. If any theoretical or phenomenological discrepancies that can seriously affect the simulation of the entire system are found, we seek a possible modification of the existing model and, by extension, a new hypothesis.  Thesis Organization From the viewpoint of the integrative and computational modeling approach, this thesis investigates aspects of three layers of the neuromusculoskeletal system: muscle, the musculoskeletal system, and the sensorimotor behavioral system. In each layer, we choose a specific dataset and build a corresponding computational model that can realistically simulate it. Suggested models are developed by extending an existing model or by introducing a new hypothesis. Figure 1.1 provides an overview of studies presented in the thesis. In Chapter 2, we focus on muscle mechanics. We begin with a brief review of muscle mechanics and its current problems, after which we investigate the me6  chanics of muscle shortening. Focusing on the shortening dynamics of muscle, we compare two feasible forms of the Hill-type muscle model: the “force scaling model” that was originally proposed by Zajac [203] and is widely used in present studies, and the relatively obscure “f-max scaling model” proposed by Abbott and Wilkie [2] in the 1950s. We evaluate their abilities to simulate isotonic shortening profiles of two different types of muscle, the mouse soleus and the toad depressor mandibulae. Contrary to our expectation that the force scaling model would better predict the data, our results show that the f-max scaling model does a much better job of predicting the slow phase of the shortening profile with high accuracy, while the force scaling model causes a significant error in prediction. This supports the claim of the thesis that using a standard model without re-investigation can be problematic. In the latter part of Chapter 2, we investigate other mechanical aspects of muscle that are not well-explained by existing models. After a brief review of the problems of current muscle models, we propose a new hypothesis called the “winding filament model” that incorporates the role of titin in active force development. After describing the main idea of the model and its explanatory power, we develop a sarcomere-level computational model that includes mechanical winding of titin on the thin filament. As a result, we show that the proposed model can simulate residual force enhancement, a well-known problem in muscle mechanics. In Chapter 3, we address an issue in the modeling of musculoskeletal system. We focus on the suitability of the line segment-based representation of muscle geometry that is widely used in musculoskeletal models. To validate, we propose a novel “top-down” optimal design framework as an alternative to the traditional “bottom-up” framework. The suggested framework finds optimal design parame7  ters that best simulate the effective mechanical output of the muscle at the limb. Our implementation result of modeling two muscles in the rat hind limb, the posterior part of the biceps femoris and the gracilis posticus, shows that the line segmentbased muscle representation does a reasonably good job of predicting the measured six dimensional force and moment field generated in situ at the ankle, but the estimated optimal insertion point for the biceps significantly deviates proximally compared to the traditional choice that takes the anatomical center of the insertion area. Interestingly, this result led to our later anatomical examination of the insertion area of the biceps, where we found a special tendon structure that routes a large portion of the muscle fibers toward the proximal side. This exemplifies the active role of computational modeling to help understand the system. In Chapter 4, we switch our scope to the behavioral level of movement. We present a model of visuomotor coordination that, for the first time, can simulate the realistic behavior of object interception in full three dimensional space. Based on our preliminary observations on object interception behavior using motion capture and head-mounted eye tracking, we hypothesize a “shared motor program” that generates discrete motor commands for both the eye and hand, by which the coordinated behavior of the eye and hand observed in the data can be explained and can also be efficiently simulated. From the simulation results, we show that the proposed model can simulate not only the existing data, but can also simulate realistic behaviors under non-trivial interception scenarios such as intercepting bouncing or disappearing objects, which emphasizes the power of generalization of the computational modeling approach. It may be too early to propose an unified integrated model that embraces the entire hierarchy of the system, and therefore the conclusion suggested in each in8  dividual chapter cannot be completely amalgamated. Nevertheless, we believe that this multifocal view has excellent potential toward the development of a unified, globally valid model of the neuromusculoskeletal system.  9  Chapter 2  Modeling Muscle Mechanics 2.1  Background  Muscles are the actuators of the neuromusculoskeletal system. In response to the electrical stimulation delivered by motor nerves, coordinated contractions of muscles generate complex biological movement. The physiology of muscle contraction has been extensively studied over centuries, but its detailed mechanism and operating principles are still not clear. Despite our lack of understanding, computational models of muscles are used in various areas such as orthopedic surgery planning [32], development of prosthetic devices [49], and realistic computer animation [107, 187]. Most studies focus on large-scale simulation of multiple muscle systems, where the reliability of the muscle model itself is less considered. However, error originiating from incorrect muscle models should not be overlooked in higher level modeling procedure. Since the modeling error usually propagates in a bottom-up manner, the prediction error of the entire neuromusculoskele10  tal model is likely to be most sensitive to the errors in the lowest, muscle level modeling. Therefore, understanding the mechanical behavior of muscle and building a reliable model of muscle contraction is by far the most critical concern in this thesis. In this chapter, we first provide a brief review of the physiological structure of muscle and terminology that will be used in the rest of this manuscript. Next, we discuss the well-known Hill-type muscle model with two forms to represent the general dynamics of shortening muscle. After that, we present our work of quantitatively evaluating the two forms using experimental data. While the suggested model for shortening muscle falls into the category of the traditional Hill-type model, there are many other mechanical properties of muscle, such as the lengthening behavior, that are hard to explain with existing models. We first discuss the currently unsolved problems of muscle mechanics and introduce a new model, called the winding filament model, as an alternative model. Using the winding filament model, a computational simulation of the residual force enhancement is then presented.  2.1.1  Anatomy and Physiology of Skeletal Muscle  Skeletal muscles consist of bundles of striated muscle fibers. A muscle fiber is a single cell that consists of smaller fiber-like components, called myofibrils. Under polarized light, each myofibril shows a repeated stripe pattern of dark and bright regions, named A and I bands respectively, based on their different refractive indices (second row in Fig. 2.1). Based on this repeated pattern, a myofibril is divided into serially connected elements called sarcomere, which is the smallest building block of the muscle contraction mechanism. It consists of an A band and two half 11  I bands on both sides. The A-band has a lower refractive index inner band called the H-zone. The center of the H-zone has a greater refractive index and is called the M-line, which is the border between two adjacent half sarcomeres. The I-band also has a higher refractive index inner band at its center called the Z-disk, the border between two neighboring sarcomeres. Structures of the myofibril and the sarcomere are illustrated in Fig. 2.1. The observed stripe pattern originates from an interwoven structure of two myofilaments, called the thin and thick filament. Muscle force is generated by the mechanochemical interaction between the thin and thick filament; the sliding of one filament with respect to the other. This is called the “sliding filament model” [130]. The thin filament consists of a double-helical twist of F-actin with additional proteins, called troponin and tropomyosin, that are embedded along the grooves of the helices. The thick filament consists primarily of myosin whose heads sprout radially from the filament axis, as shown in Fig. 2.1 (bottom row). In the overlapped area between thin and thick filaments, cross-bridges are formed. A cross-bridge is a transient connection formed between the binding site on the thin filament and a myosin head. Repetitive attachment and detachment of the cross-bridges generates the contractile force. The detailed molecular mechanism of cross-bridge cycling is beyond the scope of this thesis, and we refer the reader to the classical book by McMahon [123]. Instead, we provide the following brief summary of the cross-bridge based force generation mechanism: 1) When calcium ions, Ca2+ , flow into the sarcomere, they bind to and deform the troponintropomyosin complex on the thin filament to reveal the binding site on the thin filament. 2) After the binding site becomes available, the myosin head on the thick filament is attached to the binding site and applies force to the thin filament by 12  Figure 2.1: Inner structure of myofibril and sarcomere (figure reproduced from [7]) releasing adenosine diphosphate (ADP) that was previously attached to the myosin head. This is called the cross-bridge powerstroke. 3) After the powerstroke, a free adenosine triphosphate (ATP) is attached to the myosin head and then hydrolyzed to ADP; this hydrosis “reloads” the myosin head back to the configuration with higher potential energy that will be used in the next cycle. 13  2.2  Simulating the Dynamics of Shortening Muscle Using Hill-type Model  Based on the cross-bridge theory described above, mathematical models of crossbridge mechanics have been proposed [78]. However, most of the cross-bridge based models generally suffer from a large number of parameters and tuning functions, due to their probabilistic nature. For this reason, phenomenological models that capture the mechanical behavior of the muscle with relatively few parameters are important for large-scale musculoskeletal simulation. Rather than focusing on the detailed molecular-level mechanism, phenomenological models view the muscle as a standard mechanical system, the dynamics of which can be represented by mechanical parameters such as force, length, velocity, and control input (activation). As a dynamical system, the contractile force produced by activated muscle depends on its length and velocity. Starting from A.V. Hill’s seminal work [65], the force-length (FL) and force-velocity (FV) relationships are well established [16, 40, 123, 203], and phenomenological models based on this mechanical realization are inclusively called the Hill-type muscle model. Usually, FL and FV are obtained under static conditions by fixing one property while varying the other: FL is obtained during isometric contractions when the shortening velocity is zero, and FV is obtained by measuring the instantaneous velocity when the muscle is at its optimal length, the length at which the muscle produces its maximum isometric force. If the FLV relationship forms a parametric surface in length-velocity-force space [17], then these methods only estimate the cross-sectional shape of FLV with respect to certain values of length or velocity.  14  Therefore, the FLV behavior of muscle may not be fully described by independent models of FL and FV. In order to simulate the mechanical behavior of a muscle in the dynamic situation when the length and velocity change simultaneously, the following questions must be addressed: How can we combine the FL and FV relationships to represent the dynamic FLV relationship? Can the FLV relationship be modeled sufficiently with as few parameters as FL and FV? These important and practical questions have received surprisingly little attention compared to the independent models of FL and FV. For this reason, the goal of our study in this section is to investigate how length and velocity evolve together under physiological conditions and to develop a feasible phenomenological model that describes the FLV dynamics of muscle during shortening. In the following sections, we first clarify and define two forms of the FLV model, force scaling and f-max scaling models, and we evaluate their abilities to predict experimental data.  2.2.1  Models of FL and FV Relationships  The classical work of Gordon et al. [55] showed that the FL relationship of an isometrically contracting muscle fiber can be well modeled as a piecewise linear curve as shown in Fig. 2.2, where the force produced by the fiber is proportional to the overlap between thin and thick filaments. On a whole muscle level, this piecewise linearity tends to be smoothened by fiber inhomogeneity as shown in Fig. 2.3(a), but Winters et al. [194] recently showed that the basic shape does not differ significantly from the single fiber FL relationship. A typical FV relationship of muscle is shown in Fig. 2.3(b). Specifically, the FV relationship of a shortening muscle, the left half of Fig. 2.3B, is known to be 15  %\begin{figure}[!htb] % \centering % \includegraphics[width=0.9\columnwidth]{Figures % \caption{The force-length and force-velocity relati % \label{fig:muscle_flfv} %\end{figure}  Figure 2.2: The force-length relationship of sarcomeres in frog muscle (from [55]) force normalized by isometric force  force active passive ascending arm  0.5 L0  plateau  1.8 descending arm  L0 = rest length  Shortening  length  (a) force-length relationship  Vmax: Maximum  1.0  Lengthening  0.0  shortening velocity  velocity  (b) force-velocity relationship  Figure 2.3: The force-length and force-velocity relationship of the skeletal muscle accurately modeled by Hill’s characteristic equation [65]. The original form of the equation is 0 (F + a)(V + b) = (Fmax + a)b,  (2.1)  where F is the contractile force that the muscle generates at shortening veloc-  16  0 ity V , Fmax is the maximum isometric force, and a and b are constants. Note  that this equation describes the FV relationship when the muscle is at its optimal length. Experimental methods, such as quick-release, after-loaded and controlled release [37, 67, 86], mainly analyze the initial part of the velocity response that occurs immediately after the activated muscle is released from its optimal length under external load. The maximum shortening velocity Vmax , also called the “unloaded shortening velocity”, is achieved when there is no external load. Vmax can 0 . be obtained by setting F to zero in Eqn. 2.1: Vmax = ab Fmax 0 By normalizing force and velocity with respect to Fmax and Vmax , we get a  dimensionless version of the characteristic equation:  f=  1−v 1 + kv  (2.2)  where f and v are normalized force and shortening velocity, respectively, and k is a b dimensionless parameter defined as k = Vmax =  a 0 . Fmax  This parameter determines the  curvature of the FV relationship and is a well regulated property; in most muscles, k lies between 0.15 and 0.25 [123].  2.2.2  Two forms of FLV relationship  0 in Eqn. 2.1 is isometric force at the muscle’s optimal As mentioned earlier, Fmax  length. In other words, the equation is not defined when the muscle is at a different length, where additional assumptions about the FLV relationship are needed. When the muscle is not at its optimal length, the isometric force F 0 will be 0 by definition. It would be reasonable to assume that the resultant smaller than Fmax  force F will also be depressed in a similar way. We can define the ratio f 0 =  17  F0 0 Fmax  as the normalized FL relationship (hence f 0 ≤ 1). Then, the simplest assumption is that the force F will be scaled by f 0 . By putting the scaled force into the original equation, we obtain (  F 0 + a)(V + b) = (Fmax + a)b, f0  (2.3)  which we refer to as the “force scaling” model. Another way of modeling the decrease in force is by assuming that the muscle’s ability to produce isometric force is scaled down by f 0 , which is analogous to 0 with F 0 in Eqn. 2.1: replacing Fmax  0 (F + a)(V + b) = ( f 0 Fmax + a)b = (F 0 + a)b,  (2.4)  which we refer to as the “f-max scaling” model. Many other forms are possible as long as they produce Eqn. 2.1 at the optimal length, but these are, as far as we know, the only two forms used in previous studies. The dimensionless versions of the two models are as follows: force scaling f=  f 0 (1 − v) 1 + kv  (2.5)  f0 −v . 1 + kv  (2.6)  f-max scaling f=  These two forms have distinct properties, and both models are supported by different experimental results. First, the unloaded shortening velocity Vmax is invariant in the force scaling model, whereas it is scaled by f 0 in the f-max scaling model. Also, the dimensionless parameter k is invariant in the force scaling model 18  whereas it is scaled by  1 f0  in the f-max scaling model. The invariance of k was first  predicted by Hill [66], but he did not provide any direct experimental evidence. The invariance of Vmax has been investigated more carefully, since it is one of the main predictions of the cross-bridge model. Edman [36] presented an important result from a single fiber experiment showing that Vmax is nearly constant over a wide range of fiber lengths, from 75% to 135% of optimal length. Based on this observation, Zajac [203] proposed a phenomenological model of muscle mechanics, which is widely considered as the generic form of the Hill type muscle model in musculoskeletal simulation studies [30, 33]. Eqn. 2.5 is usually described as a multiplication of the normalized FL ( f 0 ) and FV (Eqn. 2.2) relationship. A similar argument was provided in [40]. While Vmax and k are invariant in the force scaling model, the f-max scaling model postulates that the dimensioned parameters a and b in the characteristic equation are invariant. Experimental evidence for this was presented by Abbott and Wilkie [2], who showed that a and b remain constant in the shortening region, up to 10cm from the resting length of the frog sartorius muscle. Their result was acknowledged later by Hill [69]. Another important property of the f-max scaling model is that it predicts a linear decay of the shortening velocity under constant load when the FL relationship is linear. This linear decay of the velocity results in an exponential shaped shortening curve, which is observed in many studies [15, 92]. For example, Syyong et al. [169] reported that the length profile of the isotonic shortening of airway smooth muscle shows an exponential shape. To model the observed shortening, they used the f-max scaling model for the FLV dynamics without comparison with the force scaling model. 19  Due to the conflicting results, it is difficult to determine which model better represents the mechanical behavior of muscle. In addition, many early studies were done with a skinned fiber preparation at a low-temperature and it is not trivial to assume that properties reported in those studies can be well scaled to whole muscles at physiological temperatures. Furthermore, compared to highly sophisticated models of cross bridges, these models are considerably simplified and it would make little sense to choose one model over another based on partial observations or theoretical predictions. For this reason, we need to investigate the two models more carefully, focusing on their ability to predict the shortening pattern recorded from the whole muscles working under physiological conditions. Shortening patterns of a fully activated muscle can be obtained by load-clamp, also known as quick-release, experiments: a muscle is held at a constant length while it develops tension and released suddenly to a constant load.  2.2.3  Data Collection  Muscle Preparation for Mouse Soleus To extract muscles, mice were anesthetized with 1-2% isoflurane mixed with oxygen delivered at 1-1.5L/min. The soleus muscles were removed from both the right and left legs of each mouse. Using 4-0 silk suture, the muscles were tied off securely at the junction between the tendon and the muscle to minimize the contribution of extramuscular connective tissue to the experiments (shown in Fig. 2.4 (left)) . Following removal of the muscles, mice were euthanized by cervical dislocation. Muscles were immersed in a mammalian Krebs-Ringer bath (in mM: 137 NaCl, 5 KCl, 1 NaH2 PO4 , 24 NaHCO3 , 2 CaCl2 , 1 MgSO4 , and 11 dextrose,  20  Figure 2.4: Experimental setup for mouse soleus: (Left) Mouse soleus muscle dissected and tied off using silk suture. (Right) Dual servomotor force lever and electric stimulator pH 7.4) maintained at 23-25o C and buffered with 95% O2 and 5% CO2 . The distal end of the muscle was attached to an inflexible hook and the proximal end was attached to a dual servomotor force lever (Aurora Scientific, Inc., Series 300B, Aurora, ON Canada), which provided position and force information. Fig. 2.4 (right) the overall experimental set up. In Vitro Data Collection of Mouse Soleus Contractile and elastic properties were measured in vitro from soleus muscles. The muscle length (Lm ) was adjusted such that the soleus was taut and the passive muscle tension was barely detectable by the force lever (“slack length”). Using digital calipers, Lm was measured from the proximal origin to the distal myotendinous 21  junction. Muscles were stimulated by an electrical field generated between two platinum electrodes connected to a Grass S48 stimulator. All muscle fibers were activated simultaneously with a 1 ms square pulse at supramaximal intensity (80-90 Volts). To obtain optimal muscle length (Lo ), the muscle length was adjusted until a single pulse elicited maximum twitch force, then held constant for the remainder of the experiment. It is known that the optimal muscle length of a single twitch is usually longer than that of the tetanic contraction, but the differences are generally small: less than 5% for cat soleus muscle [144] and around 5% for mouse soleus muscle (Cinnamon Pace, personal communication). A pulse train at 800-1000 ms and a maximum stimulus frequency of 70-75 Hz was delivered to the muscles to 0 ). obtain maximum isometric force (Fmax  A series of load-clamp experiments were used to record length profile of a muscle during isotonic shortening (Fig. 2.5). Muscles were stimulated tetanically to , prior to rapid step-decreases in load ( 2-95% ), and the change in length was measured for 0.8 1.0 seconds. Data Collection of Toad DM Muscle For description of muscle preparation and data collection of toad DM muscles, we refer the reader to [105].  2.2.4  Data Processing and Estimation of FV, FLV, and FL Relationships  Mouse Soleus Muscle Data Fig. 2.5 shows raw data obtained from the load-clamp experiment on mouse soleus muscle. The length profile under constant load consists of a fast phase (Fig. 2.5D),  22  A  0.16  force (N)  0.12 0.15  B  0.08 0.14  0.04 0.13  0  length (mm)  C  400  800  1200  1600  800  1200  1600  13  12 fast phase  11  D  slow phase  10 0  400  time (ms)  Figure 2.5: Load-clamp data from mouse soleus muscle. (A) Muscle force was developed isometrically until 800msec and then quickly dropped down to a constant load. (B) Trial-by-trial variaion of the force development is observed. (C) Corresponding length profile consists of fast phase (oscillatory patterns in C) and slow phase. (D) The recoil of the series elastic element and initial velocity of the contractile element were measured by extrapolating the slow phase back to the release time (grey lines).  23  a transient oscillatory response due to recoil of the series elastic element, and a slow phase (Fig. 2.5C), a low-frequency persistent response from the contractile element that is thought to be mainly governed by the FLV dynamics [2]. Behaviors of series elastic and contractile elements during the load-clamp experiment are illustrated in Fig 2.6.  length  contractile element  series elastic element  release  A  elastic recoil  B  C  D  E  time  Figure 2.6: Illustration of the load-clamp experiment from the viewpoint of the Hill-type muscle model: (A) Muscle is activated from its resting state. (B) The isometric force from the contractile element will stretch the series elastic element. (C) When the muscle is released to a lower load, the series elastic element recoils instantaneously with a transient oscillation (fast phase). (D) The contractile element is shortened (dotted curves) with a relatively slower speed (slow phase), while the length of the recoiled series elastic element remains constant. (E) The shortening will stop when the contractile element is shortened to a length where the isometric force equals to the load. First, we used the traditional FV analysis to estimate the initial velocities of the slow phase responses. As shown in Fig. 2.5D, we linearly extrapolated the slow phase back to the release time to estimate the recoil of the series elastic element 24  and velocity of the slow phase. We applied the estimated initial velocities and 0 of Eqn. 2.1 forces to the characteristic equation to find the parameters a, b and Fmax  that minimize the fitting error. However, this may produce unexpected prediction errors since, even if we keep the total stimulation duration sufficiently short, there is trial-by-trial variation of the isometric forces as shown in Fig. 2.5B. Although the observed variation is relatively small with respect to the total force, it could produce large errors in simulating the slow phase response since small error in the initial velocity can be propagated to a large error of the slow phase simulation. Therefore, we compensated for this variation by dividing every force term by the developed isometric force at release. After that, we apply the fitting algorithm 0 . suggested by Wohlfart and Edman [195] to find the parameters a, b and Fmax  Since each load-clamp trials consists of a relatively longer duration of the stimulation (1.0 to 1.6 seconds, as shown in Fig. 2.5), it is difficult to get the FL relationship from the same muscle due to fatigue. Therefore, rather than measuring the FL relationship directly, we instead parameterize the FL relationship and find the optimal parameters together with the FV parameters. Since the FL relationship of whole muscle is known to be well modeled as a piece-wise linear curve [194] and the load-clamp experiment is mainly done on the ascending limb of the FL curve, we used a simple two segment piecewise linear curve as a FL model, in which the unknown parameters were the X-Y positions of the junction (often called the “shoulder”) and the two slopes on either side. Using the FL and the FV relationships, we simulated the slow phase using two models, force scaling and f-max scaling, to evaluate their predictions. Eqn. 2.3 and 2.4 are first order ordinary differential equations (ODEs). From the initial length, the slow phase can be simulated by integrating the ODEs using a standard adaptive 25  time-step solver (ode45 in Matlab). By comparing the simulated and the observed slow phase, we find the aforementioned shoulder parameters of the FL curve using a trust region reflective nonlinear least square solver (lsqnonlin in Matlab). For each muscle, we evaluate the fit quality by calculating the R2 value of predicting the velocity profiles in the slow phase region. Toad DM Muscle Data To analyze the FLV dynamics, we assume that the fast phase response can be ignored in the analysis since it is from the series elastic element. However, since the traditional concept of the series elastic element was refuted by the cross-bridge theory [78], it is not clear that the fast phase and slow phase are completely separable. In other words, there is no guarantee that the fast phase vibration does not have any effect to the slow phase response. For this reason, we apply the same analysis to a different type of muscle that shows a much bigger fast phase response to check whether our comparison result is altered by the fast phase vibrations. Estimating the initial velocity using extrapolation is only valid when the magnitude and duration of the fast phase are sufficiently small. For this reason, we found it hard to apply the same FV analysis to the toad data that show large and long-lasting fast phase oscillations (see Fig. 2.11). Moreover, unlike the in vitro set up for the mouse muscle, we observed larger measurement uncertainties on length and force measurements due to the in situ set up of the experiment. For this reason, we normalized all length and force data by normalization using the measured optimal length and maximum force reported in Lappin et al. [105] and estimated k and Vmax in Eqn. 2.5 and 2.6 together with the parameters of the FL curve during the optimization. Since the initial recoil length is also uncertain, we ignored the error 26  estimating the initial length by taking only the variance of the error in slow phase estimation.  2.2.5  Results 1  muscle 1 muscle 2 muscle 3 muscle 4  normalized force (Fmax)  0.8  k 0.14 0.18 0.16 0.20  R2 0.998 0.985 0.993 0.996  0.6  0.4  0.2  0 0  0.2  0.4 0.6 normalized velocity (Vmax)  0.8  1  Figure 2.7: Force-velocity (FV) relationship of four mouse soleus muscles estimated from the measured initial velocity: Circles show data points and curve shows the estimated FV relationship using Hill’s characteristic equation. Fig. 2.7 shows the FV relationship of the mouse soleus muscle along with the estimated dimensionless parameter k and the corresponding R2 . As shown, the FV relationship is well predicted by the characteristic equation and the range of k agrees well with previous studies. We found that the FV relationship of muscle 1, which was released before it reached the full isometric tension, resulted in a  27  slightly smaller k, but the difference was minimal.  force scaling 2  2  muscle 2 (R = 0.883)  muscle 1 (R = 0.804) 13  13  12  12  11 11 10 10 9 9  0  200  400  600  800  800  1000  2  1400  1600  2  muscle 3 (R = 0.860)  14  1200  muscle 4 (R = 0.803) 13  length (mm)  13 12 12 11  11  10  10  9 9  800  1000  1200  1400  800  1600  1000  1200  1400  1600  time (ms)  Figure 2.8: Estimated slow phase responses of four mouse soleus muscles using force scaling model: Black curves are length profiles recorded from the load-clamp experiment and dotted red curves are simulated length changes from the model. Fig. 2.8 and 2.9 show simulation results for the mouse soleus slow phase responses using force scaling and f-max scaling models. The results clearly show that the f-max scaling model gives a better fit than the force scaling model. Corresponding estimates of the FL relationship for the f-max scaling model are shown in Fig. 2.10 (right); the estimated FL curves qualitatively agree with the ascending arm of the classical FL relationship. The differences of FL shapes were larger  28  f-max scaling 2  2  muscle 2 (R = 0.969)  muscle 1 (R = 0.963) 13 13  12  12  11  11  10 10  9 9  200  400  600  800  800  2  1200  1400  1600  2  muscle 3 (R = 0.983)  14  1000  muscle 4 (R = 0.972) 13  length (mm)  13 12  12  11  11  10  10  9 9 800  1000  1200  1400  800  1600  1000  1200  1400  1600  time (ms)  Figure 2.9: Estimated slow phase responses of four mouse soleus muscles using f-max scaling model: Figure legends are same as Fig. 2.8. than those observed in the FV analysis, but they seem to be within the range of the differences in measured FL curves reported in the previous study [192]. Fig. 2.11 shows the simulation results for the toad DM muscle data. The result also shows that the f-max scaling is better at predicting the real data. We can observe the tendency that the force scaling model cannot predict the gradual decrease of the velocity of the shortening under lower loads. In addition, the estimated initial slow phase responses during the fast phase oscillations, which were not considered during the optimization, are highly deviated from the approximated center line of  29  force scaling muscle 1 muscle 2 muscle 3 muscle 4  1 normalized isometric force (Fmax)  f-max scaling  0.8  0.6  0.4  0.2  0  10  11  12  13  9  10  length (mm)  11  12  13  length (mm)  Figure 2.10: Estimated force-length (FL) curves for four mouse soleus muscles using force scaling and f-max scaling model: Each curve is modeled as two connected line segments. the oscillations in the force scaling model, but are described well by the f-max scaling model. Finally, as shown in Fig. 2.11 (bottom row), the estimated FL curves are different, with the ones from the f-max scaling appearing more physiological. To check whether the identified model overfits the data, we performed leaveone-out cross validation: for each muscle, we excluded a single trial and checked how well the model built by the rest of trials predicts the excluded trial. By combining individual predictions on the omitted trials, we reconstructed the entire prediction on the whole muscle data and evaluated its R2 value. If the model overfits data, we will observe a large drop of the R2 value. Table 2.1 shows the cross validation results for all muscles used in the analysis. Compared to the R2 values reported in Fig. 2.8, 2.9 and 2.11, we observed no significant difference in fit quality, which implies that our optimization methods for both mouse and toad are not likely to introduce the overfitting. Relatively larger drops in R2 values of force scaling (average 6.71%) compared to those of f-max scaling (average 1.25%) also suggest that parameters in the f-max scaling model  30  force scaling  f-max scaling  2  muscle 1 (R = 0.928)  2  muscle 1 (R = 0.941)  1  1  0.9  0.9  0.8  0.8  0.7  0.7  200  250  300  350  400  200  2  muscle 2 (R = 0.773)  0.9  0.9  0.8  0.8  0.7  0.7  300  350  400  200  2  length (L0)  400  250  300  350  400  2  muscle 3 (R = 0.891)  muscle 3 (R = 0.944)  1  1  0.9  0.9  0.8  0.8  0.7  0.7  200  350  2  1  250  300  muscle 2 (R = 0.965)  1  200  250  250  300  350  400  200  250  300  350  400  time (ms) 1  force (Fmax)  0.8  1 muscle 1 muscle 2 muscle 3 0.8  0.6  0.6  0.4  0.4  0.2  0.2  0 0.6  0.7  0.8  0.9  1  0 0.6  muscle 1 muscle 2 muscle 3  0.7  0.8  0.9  1  length (L0)  Figure 2.11: Estimated slow phase responses and corresponding FL curves of the toad DM muscle: Force scaling (left column) and f-max scaling model (right column) are compared. Figure legends are same as Fig. 2.8 and 2.10. 31  mouse  toad  muscle 1 2 3 4 1 2 3  force-scaling 0.777 0.866 0.817 0.785 0.839 0.740 0.717  f-max scaling 0.956 0.967 0.980 0.970 0.923 0.944 0.910  Table 2.1: Cross validation result of the FLV models for all data: Each entry represents the R2 values of the reconstructed prediction combined from individual leave-one-out predictions. are better constrained by the data. Taken together, both results taken from data in mice and frogs strongly suggest that f-max scaling is a better model of the FLV relationship in isotonic shortening.  2.2.6  Discussion  The results of these simulations clearly support f-max scaling as a better representation of the FLV relationship than the popular force scaling model. This is surprising, since force scaling is widely used and regarded as a textbook model of muscle biomechanics [40, 203]. In particular, prediction by the force scaling model consistently overestimates the shortening velocity at shorter lengths. This may imply that the cross-bridge theory’s prediction that the unloaded shortening velocity is independent of muscle length, also the main basis of the force scaling model, does not hold at shorter lengths. Since the operating ranges of vertebrate muscles cover almost the entire ascending arm [19], using the force scaling model in large scale musculoskeletal simulation studies, such as locomotion studies, could produce significant errors in dynamic predictions. In fact, Edman [36] reported that the unloaded shortening velocity decreases 32  at a length shorter than 75% of the optimal length. This also corresponds to Brown and Loeb’s model [16] that there is an additional spring element that generates a “negative” pushing force in the shortening region of the muscle. If there is an additional element or mechanism in the ascending region that slows down the maximum shortening speed, it could provide a possible reason why the forcescaling model does not work on the ascending limb. Another possible cause of the slow shortening speed is that the length change of the muscle may have caused a shortening-dependent or work-dependent deactivation [63, 193]. Rather than introducing additional models of the deactivation mechanism or negative force, our simulations suggest that it may be preferable to use a simple f-max scaling model in large scale simulations. From a mechanical viewpoint, it is interesting to speculate on the exponential shape of the shortening profile. The observed linear decay of the velocity shows that the muscle actually works as a linear spring-damper system when it shortens, although its viscous property is commonly described as a nonlinear, multiplicative damper. Considering the significant advantage of controlling a system with linear dynamics, this inherent mechanical property of muscle may help simplify the neural control of movement. In spite of the good prediction ability of the f-max scaling model, our result should not be over-generalized to the entire range of FLV dynamics. Evaluation of the model in broader and more realistic situations such as work-loop analysis [89] or locomotion analysis [157] including the lengthening contraction and frequency response remains as future work, as does the incorporation of innervationactivation.  33  2.3  Winding Filament Hypothesis  So far, we have shown that the shortening dynamics of muscle can be sufficiently well modeled using a Hill-type model with a small modification of the standard Hill-Zajac model. As we have discussed, the Hill-type model simplifies the complicated molecular level mechanism of cross-bridges into a simple constitutive equation with a high quality of fit. For this reason, the Hill-type model has been successfully applied to many studies in biomechanics and motor neuroscience. Some thermodynamic and transient mechanical properties that cannot be explained by Hill’s original model were explained using the mechanical model of the crossbridge cycling subsequently suggested by Huxley and colleagues [75, 76, 78]. The cross-bridge mechanism has been intensively studied and more complicated models beyond the original two-stage model have been proposed. For the state of the art studies in cross-bridge modeling, we refer the reader to a recent book by Rassier [145] or other review papers [79, 95]. Although the Hill-type model and the cross-bridge model have been enormously successful in revealing the fundamental mechanism of the muscle contraction and much experimental evidence to support these models have been suggested, researchers also have indicated many problems that cannot be explained by either model. In the following paragraphs, we briefly summarize the current problems of muscle mechanics.  34  2.3.1  Currently Unsolved Problems in Muscle Mechanics  Sarcomere Instability, History Dependence, and Residual Force Enhancement One of the most obvious and also the most studied problems is the mechanical instability of a system of serially connected sarcomeres. Hill’s model assumes that for any given length, the force is always uniquely determined in all situations by the FL relationship. In other words, for a fully activated muscle, force generated at a given length is assumed to be always determined by the FL curve. However, the FL curve has negative slope on its lengthening segment, called the descending limb, as shown in Fig.2.3(a). Since the slope of FL relationship is defined as “stiffness”, this negative stiffness of FL relationship predicts that a sarcomere on the descending limb produces less force when it is lengthened. If we consider a set of serially connected active sarcomeres whose lengths are all on the descending limb, and if we perturb a single sarcomere so that its length becomes slightly longer than the others, then it will become weaker than others due to the negative stiffness. Since this weakened sarcomere cannot sustain the force generated by neighboring sarcomeres, it will be stretched more and therefore becomes even weaker. As a result, the perturbed sarcomere will not go back to its original length. Since small differences in sarcomere lengths are amplified, the Hill-type muscle model predicts that whole muscle becomes unstable on the descending limb. It is important to notice that this instability of the sarcomere is not just a theoretical defect that scarcely happens in reality, but a very practical problem in muscle mechanics. It is commonly argued that, for some muscles, the passive elastic element with a positive stiffness compensates the negative stiffness on the  35  descending limb. However, even with the passive force, fusiform muscles such as sartorius still have an unstable region. Further more, Azizi and Roberts [5] showed that the frog plantaris muscle operates primarily on the descending limb with very small passive force. To resolve this problem, most simulation studies that include muscle as an actuator use ad hoc modifications of the FL curve or limited operation range to prevent the simulation from becoming numerically unstable [107]. Morgan [60] suggested a hypothesis called “sarcomere popping” to explain the instability problem. The hypothesis argues that the instability between sarcomeres is a real phenomenon, but the stability of the whole muscle is still achieved by non-uniformly distributed sarcomere lengths. Specifically, the hypothesis predicts that some unstable sarcomeres “pop” to extreme lengths while other sarcomeres stay at constant length or shorten. Morgan’s argument has been one of the popular hypotheses to resolve the problem and it is still under debate. However, after experiments on a single myofibril, by which the individual behavior of the sarcomere under mechanical stretch can be monitored, became available, Morgan’s argument has received less attention since no study has observed sarcomere popping [90, 171]. Based on the instability problem, the question of whether the FL curve is a real constraining curve of the muscle dynamics has arisen. After the initial observation by Abbott and Aubert [1], many researchers observed that muscle has significant amount of history dependence [38, 140, 147]. In contrast to the traditional Hilltype model in which the muscle force is determined by the current states such as length and velocity, history dependence means that the generated contractile force is affected by the past history, such as when the activation starts and how much the muscle has been stretched or shortened after the activation onset. A simple 36  comparison can be made by two experiments: 1) stretch an active muscle while maintaining activation, and 2) stretch a passive muscle and then activate it. Even if the final lengths and activations are the same, the steady state forces produced by these two experimental protocols are found to be quite different and the difference becomes larger on the descending limb. The conclusion is that the produced force does not always simply follow the FL relationship. Since stretch with constant activation produces a larger steady state force, the effect is named the “residual force enhancement”. This observed residual force enhancement of the muscle naturally resolves the instability problem: under activation, the force always increases upon strechting even on the descending limb. However, as discussed, it is clear that the Hilltype muscle model cannot explain the residual force enhancement. Investigating whether more sophisticated cross-bridge based models can account for the residual force enhancement, recent work of Walcott and Herzog [185] used mathematical modeling to show that residual force enhancement cannot be explained without violating the basic assumptions of the cross-bridge theory. Rather than adding more ad hoc assumptions to the cross-bridge models, it would be more intuitive to assume that there is an unmodeled mechanical element that keeps the positive stiffness of the system under activation and therefore stabilizes the entire muscular structure. Epstein and Herzog [40] suggested a simple mechanical model that has a constant positive stiffness after activation begins. Forcinito et al. [42] proposed a new rheological model that has an elastic rack element, a linear dashpot, and a spring to explain the residual force enhancement. However, most of these models were developed in a purely engineering sense without solid physiological bases; for instance, using ad hoc elements such as elastic 37  racks, and do not truly take into account the molecular mechanism.  Elasticity and Energetics of Active Shortening and Lengthening Studies of the FV relationship are closely related to the thermodynamics of muscle. Hill [65] first suggested that his characteristic equation can also explain the thermodynamic observation, which states that the heat production rate of a shortening muscle is proportional to the shortening speed. However, he revised his measurments later, showing that the heat production does not monotonically increase and rather that it levels off when the velocity of shortening is more than half its maximum [68]. Later, Huxley [76] proposed the two-stage model of the cross-bridge to explain the decreasing energy liberation at high speed, which is considered as a classical model of the cross-bridge mechanics. However, the model has a critical problem in that it does not specify any alternative source of the energy that is responsible for the shortening without heat production. Instantaneous elasticity of muscle during shortening is also not well explained by the previous models. As shown in the FLV modeling in the previous section, the fast phase of isotonic shortening is due to the elastic energy stored in the muscle. One of the problems of the Hill-type muscle model is the direct connection between the contractile element and series elastic element, since this connection cannot explain the instantaneous recovery of the tension after a quick shortening [77], called the “rapid elasticity” [123]. A. F. Huxley’s cross-bridge theory model did a good job of quantitatively explaining this rapid elasticity as a short-range elasticity that resides in the cross-bridges [78]. However, an important feature of the crossbridge elasticity is that it can only be valid in a very short range of the length change, approximately 6nm per half sarcomere [77]. As supporting evidence, the 38  original experiments reported that quick shortening by 6nm per half sarcomere can drop the tension to near zero. However, Galler and Hilber [47] reported, based on a later experiment, that the amount of shortening required to drop the tension to zero increases up to 18nm per half sarcomere at a higher temperature. Moreover, an experiment by Sugi et al. [168] on horseshoe crab muscle, whose sarcomere is much longer (7µm) than in typical muscle fiber but has the same cross-bridge mechanism, reported that 210nm of shortening per half sarcomeres is needed to drop the tension to zero, which is significantly longer than the stiffness range of the cross-bridges. On a whole muscle level, Lappin et al. [105] reported that the amount of fast phase shortening in toad DM muscle (shown in Fig. 2.11) is up to 21.4% of the rest length, which is also far longer than any known elasticity inside the muscle. Notwithstanding the above problems, the elasticity and energetics of the shortening part of the muscle have been relatively well studied. However, those of the lengthening part of muscle mechanics still remain surprisingly unknown, which has been even referred as a “dark side of muscle mechanics” by T. A. McMahon [114]. When the external force is greater than the muscle’s force production capability, the muscle will be actively lengthened, which is called an eccentric contraction (as opposed to the concentric contraction, the active shortening). It is well known that the muscle spends much less energy to produce the force resisting the external force during the eccentric contraction [12]. Furthermore, the muscle can store the elastic energy during lengthening and use it in the following shortening, which enables an efficient use of energy in locomotion [114, 150]. Neither the source of the force that resists high external load during lengthening nor the anatomical component that stores the elastic energy within muscle during lengthening are modeled 39  yet. In summary, there are many unsolved problems in muscle mechanics. These problems are not just ignorable side effects that inevitably arise in modeling studies; they can severely alter the stability and efficiency of the entire system. For developers of a biomechanical simulation model, it is a frustrating experience that the traditional Hill-type model cannot simulate a simple model of two serially connected sarcomeres. The absence of a reliable muscle model is undoubtedly the most serious issue of neuromusculoskeletal modeling and practically all the higher level modeling and simulation works are potentially vulnerable to error propagating bottom-up from a flawed muscle model. As a significant departure, we suggest a new model of muscle mechanics in the next section that can possibly explain the unsolved problems.  2.3.2  Winding Filament: a New Hypothesis  To build a model that is realistic from both the mechanical and molecular points of view, we propose a new hypothesis based on recent findings on the activationdependent behavior of titin and the rotation of the thin filament. This hypothesis was initiated by Dr. Nishikawa, and elaborated in collaboration with the coauthors [130]. The hypothesis can conceptually explain many aspects of muscle mechanics that have been previously described. The main idea of the winding filament hypothesis is to include titin in active force development. The mechanical behavior of titin in active muscle is hypothesized based on two new observations in molecular physiology: the calciumdependent binding of titin and the rotation of the thin filament. In the next section, we briefly introduce these ideas and discuss the potential of the hypothesis in ex40  plaining several mechanical observations. Titin as an Unmodeled Protein in the Sarcomere Titin is the largest protein found in striated muscle tissues and spans the entire half sarcomere from the Z disk to the M line. Although the existence of titin, formerly called connectin, was hypothesized in early studies [80], its real structure and properties were only revealed relatively recently, in the late 1970s [119, 188]. This is one of the reasons why titin is not included in traditional models of muscle mechanics, which were developed more than a decade before these findings. The layout of titin shown in Fig. 2.12(a) exhibits an interesting mechanical structure. In the A-band, the titin is entwined with the thick filament. At the boundary between the A-band and the I-band, six titin filaments radiate from the tip of the thick filament toward the surrounding six thin filaments and anchor to the thin filament. A portion of titin, spanning the tip of the thick filament to the anchor point on the thin filament, is considered to be a mechanically free portion of titin, and mainly consists of three distinct domains: the proximal immunoglobulin (Ig), PEVK, and distal Ig domains. PEVK is a distinct motif type that mainly consists of proline, glutamate, valine, and lysine residues [103]. The elastic properties of these different domains were clarified by Linke et al. [115]. When a single myofibril is passively stretched from the resting length, the initial stretch mainly unfolds the proximal Ig domain whose stiffness is very small at low force. At a longer length, the stiffness of the Ig domains increases rapidly and the stretch begins to elongate the PEVK domain. The elasticity of the PEVK domain follows the modified worm-like chain model of polymer elasticity: it behaves as an entropic spring at lower force and becomes an enthalpic, Hookean spring at higher force. The distal part of the Ig domain is 41  Figure 2.12: Schematic of a skeletal muscle half-sarcomere illustrating the layout of titin(yellow with a red N2A segment):(a) Each titin molecule is bound to the thin filaments (blue) in the I-band, and to the thick filaments (green) in the A-band. Note that, for simplicity, thick filaments are illustrated as double-stranded, whereas in vertebrate skeletal muscle, they are actually triple-stranded. The N2A segment is located between the proximal tandem Ig segment and the PEVK segment. (b,c) Ca2+ -dependent binding of N2A to thin filaments. (b) Resting sarcomere at slack length at low Ca2+ concentration (pCa = 9). (c) Upon Ca2+ influx (pCa = 4.5), N2A titin binds to the thin filament (blue), which shortens and stiffens the titin spring in active sarcomeres. Figure and caption are adopted from [130]. considered to be inelastic [10]. Possible Contribution of Titin in the Active Force Developement The PEVK domain of titin can be conceptualized as a mechanical spring and could be a possible candidate for the source of the spring-like property of muscle. However, as shown in the passive stretch experiment, the mechanical contribution of  42  the PEVK domain is only effective in an elongated sarcomere when the length is longer than 2.4µm, while the compliant proximal Ig domain is mainly unfolded at a shorter length. Moreover, titin has been considered as a purely passive element that does not contribute to active muscle force. This is another reason why titin is excluded in the active mechanism of muscle contraction in a normal operational range. However, many recent studies have shown that titin may have different mechanical properties in active muscle. A study by Labeit et al. [102] showed that an influx of Ca2+ in the muscle fiber increases the stiffness of titin. Since titin stiffness is changed by activation, the elastic force of stretched titin will depend on the onset time of activation. This provides a possible explanation for history-dependent behavior such as residual force enhancement. A possible mechanism for the change in stiffness is that titin binds to the thin filament by Ca2+ . An in vitro study showed that the product of titin proteolysis has high affinity for actin in the presence of Ca2+ [96], which implies that the proximal part of titin could bind to the thin filament in a calcium-dependent manner. Based on this observation, Rode et al. [153] proposed a model called the “sticky spring mechanism” that hypothesizes that the PEVK domain binds to the myosin binding sites on the thin filament and results in an increase in the stiffness and corresponding history-dependent properties. However, the reported amount of enhancement by calcium is too small to explain the observed enhancement. Joumaa et al. [91] observed the force enhancement of troponin-C depleted myofibril, by which the cross-bridge interaction is completely prevented, and reported that the influx of Ca2+ produces the force enhancement, but the amount of the enhancement is too small. 43  The above results imply that calcium increases the stiffness of titin, but without cross-bridge cycling, the amount of the stiffness change induced by Ca2+ only is too small to explain residual force enhancement. A recent result by Leonard and Herzog [112] provides definitive evidence showing that titin contributes to active force generation and its mechanical properties are regulated by both Ca2+ and cross-bridges. These authors actively stretched a single myofibril to a very long length where the thin and the thick filaments did not overlap at all, and showed that the sarcomere was still able to produce active force. Since titin is the only possible filament that can mechanically transfer force at that extreme length, the result shows that titin plays a significant role in active force generation. Another important observation of this study is that the ability of titin to produce active force disappears when strongly bound, force generating cross-bridges are inhibited by butanedione monoxime. This strongly supports the hypothesis that cross-bridge cycles regulate the active force of titin, but no previously suggested model can explain the actual mechanism of the cross-bridge based regulation of titin stiffness. Rotation of the Thin Filament and Possible Winding of the Titin The abovementioned studies show that the calcium-dependent increase in titin stiffness cannot solely explain the observed history-dependent behavior of the muscle. Here, we discuss another observation on thin filament rotation induced by crossbridge cycling and propose the winding filament hypothesis by combining the two observations. A rich group of studies has shown that the thin filament rotates in active force development. This means that the cross-bridge powerstrokes produce not only axial force but also lateral force that generates torque to twist the thin filament. This 44  rotation can be estimated by measuring the change in the helical pitch in the actin double helix structure. Nishizaka et al. [131] observed a right-handed rotation of the thin filament in vitro. Bordas et al. [14] performed an in vivo X-ray diffraction study showing that the helical pitch decreases during unloaded shortening. Based on these observations, Williams et al. [191] recently proposed a model of axial and radial force generation by cross-bridges. Taken together, by combining the two observations on the calcium dependent binding of titin to the thin filament and cross-bridge induced rotation of the thin filament, we can draw a reasonable analogy that titin may be wound around the rotating thin filament. This is the main idea of the winding filament model. Since the exact site on titin that binds to the thin filament is still unknown, we hypothesize that a segment in titin that connects the proximal Ig and PEVK domains, called the N2A segment, works as a strong site for binding to the thin filament (Fig. 2.12(b,c)). Also, we assume that the PEVK domain also binds to the thin filament by electrostatic force, similar to the sticky spring mechanism. Since the proximal Ig domain becomes un-stretchable when the N2A segment binds to the thin filament, the stiffness of titin in an active muscle will be determined by the elasticity of the PEVK domain. This explains the calcium-dependent change in titin stiffness. As the thin filament rotates, the proximal part of the PEVK domain will be wound around the thin filament and titin stiffness will be determined by the unwound portion of the PEVK domain, which possibly explains how cross-bridges regulate the stiffness of titin. Fig. 2.13 illustrates the mechanism of titin binding and the corresponding stretch of the PEVK segment. Note that the figure is only an illustration of winding and that the number of turns may be less as shown later. 45  Figure 2.13: Schematic illustrating the winding filament hypothesis(a,b): Ca2+ -dependent binding of N2A to thin filaments contributes to length-tension relationship. If N2A (red) binds non-selectively to thin filaments (blue) in the presence of Ca2+ , and if the binding site depends on sarcomere length at the time of Ca2+ influx, then a plateau is predicted in active force at sarcomere lengths between (a) 2.4 and (b) 2.6 micron in rabbit psoas muscle. (c,d) Cross-bridge cycling results in titin winding. Note that these are conceptual illustration; Actual winding pitches and total winding angles predicted from the simulation are different from these illustrations (see Fig. 2.16). (c) Cycling of the cross-bridges winds PEVK on the thin filaments (arrow indicates direction of rotation). In the model, the winding angle depends only on sarcomere geometry. (d) Stretch of an active sarcomere extends the PEVK segment and enhances the active force. Figure and caption are adopted from [130]. The winding filament hypothesis gives us a new insight into muscle mechanics by introducing the rotational degree of freedom that can store elastic energy without changing the length. Note that the hypothesis does not deny the cross-bridge model; it rather provides a supplementary mechanism that regulates the stiffness 46  of the system. The hypothesis can give potential answers to some of the unsolved problems mentioned above. The length of the PEVK domain in the active sarcomere depends on the length at which the N2A segment is bound to the thin filament. Such a mechanism can explain the observed history dependency. After titin becomes taut, the total stiffness of the sarcomere is the sum of the stiffness of the cross-bridges based on the FL relationship and the elastic stiffness of titin. Therefore, the total stiffness of the system can be kept positive and stability is preserved even when the cross-bridge force is on its descending limb. In a thermodynamic sense, the hypothesis predicts a complex exchange relationship between 1) external work, 2) ATP hydrolysis in cross-bridges, and 3) the elastic strain energy of titin. During active shortening, the strain energy of titin can either be accumulated by winding or be dissipated by shortening. When the rate of energy accumulation by winding and the rate of energy dissipation by shortening are equal, the strain energy stored in titin will remain constant. At a high shortening speed, i.e., when the shortening speed exceeds the winding speed, the work done by the muscle will mainly come from the dissipated strain energy of titin, which may be the unexplained source of energy during high speed shortening. At low speed, when the winding speed exceeds the shortening speed, the energy produced in cross-bridges will be used either for executing external work or for accumulating the strain energy of titin by winding. This mechanism can account for the fact that shortening velocity near the isometric tension increases slowly in the FV curve as tension drops (Fig. 2.3(b), the left-hand side tangential line in grey). When actively lengthening, both external work and winding will accumulate titin strain energy, which can explain the rapid increase in force in the FV relationship 47  Quality  Value (nm)  References  Thin filament radius  5  Thin filament length  1000  Burkholder and Lieber [19]  Thick filament length  800  Burkholder and Lieber [19]  21.5  Irving et al. [82]  Distance between and thick filaments  thin  Titin radius  Length of titin distal Ig  Shi et al. [162]  2  Funatsu et al. [46], Higuchi et al. [64]  100  Bennett et al. [10]  Table 2.2: Geometric parameters of sarcomere used in the modeling in eccentric contraction (Fig. 2.3(b), the right-hand side tangential line in grey).  2.3.3  Computational Model of the Winding Filament Hypothesis  Although the winding filament hypothesis provides a promising way of explaining previously unexplained observations, having a computational model and building a framework to validate its behaviors are necessary for achieving theoretical completeness. As a first step toward a fully working computational model, we present a computational model of the winding filament hypothesis in this section that explains residual force enhancement. Specifically, we develop a three dimensional mechanical model of the winding that quantifies the effects of thin filament rotation on the stiffness of titin. We first define geometries of thin, thick and titin filaments and then develop differential equations describing how the pitch of the winding and the corresponding stiffness of titin evolve as the thin filament rotates. We finally add additional models of  48  torques from the cross-bridges and Z-lattice twisting that are responsible for the rotation of the thin filament, and simulate the residual force enhancement. Three Dimensional Model of Half-sarcomere We model the sarcomere of the rabbit psoas muscle whose geometrical parameters and force-length properties have been well studied. The geometrical parameters of myofilaments per half sarcomere in the rabbit psoas muscle are summarized in Table 2.2. Using the above parameters, we can build a three dimensional geometric model of the thin filaments, thick filaments, and titins, topological relationship of which is described in the previous section. As mentioned in the previous section, the main hypothesis of the winding filament model is that the N2A segment binds to the thin filament under the presence of Calcium ions, and the torque generated by the cross-bridges rotates the thin filament and therefore winds the PEVK segment around the rotating thin filament. Therefore, to simulate the kinematics of winding, the only additional information we need is the binding location of the N2A segment to the thin filament at the onset of activation. The distance of the N2A segment from the Z-disk in a passive half sarcomere, denoted as DN2A , was reported in [115] as follows: DN2A = 134s3 − 1231s2 + 3808s − 3648,  (2.7)  where s is the sarcomere length in microns. Mechanics of Winding When N2A is attached and PEVK starts to wind around the thin filament, the geometric configuration of the winding is illustrated in Fig. 2.14. The notations defined 49  in the figure will be used in the following paragraphs. bound PEVK thin filament  fre  attachment  eP  EV  K  thick filament + distal Ig  lateral view  r  axial view  φ  xw (bound length)  h  δxw r δφ θ δd1  x (f  ree  len  gth  )  h d2-d1  d1  θ  x  d2-d1  d2  Figure 2.14: Kinematics of titin winding. As the PEVK segment is wound around the thin filament, we divide the whole PEVK segment into two parts: bound PEVK and free PEVK (Fig. 2.14, upper left). Since the model assumes that the wound part of the PEVK is rigidly bound to the thin filament due to the electrostatic force, only the free PEVK contributes to the mechanical properties of the titin. The goal of the modeling is to simulate 1) how the geometry of free PEVK, denoted as a grey triangle shown in the bottom figure of Fig. 2.14, changes as winding proceeds, and 2) how much strain is accumulated in the free PEVK.  50  Change of the free PEVK geometry can be modeled as a set of differential equations. When the thin filament rotates by a small angle δ φ , the tip of the triangle will be wrapped around the thin filament. The wrapped portion forms a similar triangle that is illustrated as a smaller triangle in the lower right corner of the Fig. 2.14. Using trigonometry and dividing the equation by δt and taking limits, the rate of change of d1 is d˙1 = (r tan θ )φ˙ =  r(d2 − d1 ) ˙ φ. h  (2.8)  Similarly, the rate of change of the bound PEVK length (xw ) is  x˙w =  rx ˙ r φ˙ = φ. cos θ h  (2.9)  It is important to note that x˙w = −x˙ because the length of x is not constant due to the stretch of the free PEVK. Instead, x˙ can be derived by trigonometry of a larger triangle in the lower right corner of the Fig. 2.14:  x2 = h2 + (d2 − d1 )2 .  (2.10)  Taking time derivatives of both side, we get  xx˙ = hh˙ + (d2 − d1 )(d˙2 − d˙1 ).  (2.11)  Since h is assumed to be constant, h˙ is zero. By replacing d˙1 using Eqn. 2.8 we get  x˙ =  d2 − d1 ˙ d2 − d1 ˙ r(d2 − d1 )2 ˙ (d2 − d˙1 ) = d2 − φ. x x hx  51  (2.12)  Note that d˙2 is given in the isokinetic experiment (especially, d˙2 = 0 for the isometric experiment). Similarly, accumulation of strain in the free PEVK can also be modeled by a differential equation. Since strain is defined by a relative change in length with respect to the initial length, we can calculate the strain by combining the length of free PEVK and its initial length. Although the initial length remains constant in most mechanical problems, we need to consider it as a variable in our case, since it is the initial length of the free PEVK, which changes as the winding proceeds. Given the length and the initial length, x and x0 , the strain ε is defined as:  ε=  x − 1. x0  (2.13)  When x and x0 are changed by δ x and δ x0 respectively, the change of the strain will be ε +δε =  x+δx − 1. x0 + δ x0  (2.14)  Given ε, x, and x0 , the change of the strain δ ε can be determined by δ x and δ x0 . Since δ x can be calculated using Eqn. 2.12, we only need to compute δ x0 . Intuitively, x0 decreases as the winding proceeds since more of the PEVK is changed from free to bound. If length δ xw is additionally bound to the thin filament by rotation, the problem is to estimate how much initial length of free PEVK will be “lost” by losing δ xw . Using current strain ε, it can be computed as follows:  δ x0 = −  52  δ xw . ε +1  (2.15)  By applying Eqn. 2.8, 2.9, 2.12 and 2.14 to 2.13, ε˙ can be expressed as  ε˙ =  1+ε {(d2 − d1 )d˙2 + rhφ˙ }. x2  (2.16)  Taken together, we can define a nonlinear ODE to simulate the winding for the given thin filament rotation profile φ (t) and the sarcomere length profile d2 (t): r(d2 − d1 ) ˙ d˙1 = φ h rx x˙w = φ˙ h d2 − d1 ˙ r(d2 − d1 )2 ˙ d2 − φ x˙ = x hx 1+ε ε˙ = 2 {(d2 − d1 )d˙2 + rhφ˙ }. x  (2.17)  Initial conditions of xw and ε are 0 and those of d1 and x are determined by the geometry at initial activation. For binding position of N2A, the initial value of d1 is given by DN2A in Eqn. 2.7.  Figure 2.15: Elasticity of PEVK in single titin (figure reproduced from [115])  53  The computed strain, ε, can be used to calculate the force generated by titin using Linke et al.’s model shown in Fig. 2.15. The only change is that the force is multiplied by 0.5 since Linke et al. assumes three titins per thick filament, whereas recent finding from the same group has modified their assumption that there are actually six titins per thick filament [100]. Since the elasticity of the titin is represented as a relationship between fractional extension and force per titin in pico-Newton in their paper, the strain ε in our equation should be converted to the fractional extension:  f ractional extension =  x0 (ε + 1), xmax  (2.18)  where x0 is the initial length of the PEVK before winding and xmax is the contour length used in [115] (476nm). Given the force of titin ft obtained from Fig. 2.15, the titin’s axial and radial components can be calculated by its geometric relationship: d2 − d1 ft , x h = r ft , x  fta =  (2.19)  τt  (2.20)  where fta is the axial force generated by titin and τt is the torque applied to the thin filament by multiplying the radial force of the titin by a constant moment arm, r, the radius of the thin filament (2nm). Combining Models of Cross-bridges and Z-lattice Twisting By combining this geometric model of the winding and the cross-bridge mechanics in three dimension, we can simulate the active force generation in sarcomere. 54  Unlike previous models that consider the cross-bridge cycling in a single axial dimension, we instead simplify and lump the cross-bridges as a pure force generator that applies a force to the thin filament in both axial and radial directions. Based on the model proposed by Morgan [125], the force applied from thick filament to thin filament is tilted by 10.038 degrees with respect to the axial direction. We can also convert the radial force to the torque by multiplying r. The magnitude of the generated force from cross-bridges is tuned based on the single myofibril study by Pavlov et al. [135], which is 300pN at plateau (2400nm). In the radial direction we need to consider the balance of torques generated by cross-bridges and titin. Here we assume that another source of the torque is from the lattice structure in the Z-disk. Consistent with the winding hypothesis, it has been observed that the lattice of alpha actinin forms a small square pattern in the cross-sectional view during the resting state, and then changes to a basket-weave pattern under force development [53, 85]. Therefore, it is reasonable to assume that change of the lattice structure will apply a resisting torque against the rotation of the thin filament. We model the relationship of twisting angle and torque of the z-lattice as a simple exponential spring. Since we have not seen any previous study that measured this relationship, the parameters were tuned by hand so that the final result was realistic. The equation used in the simulation is: τz = 10−5 e4.5φ pN · nm. To summarize, there are two types of axial forces, cross-bridge and titin, which balance the external force, and three types of torques, titin, Z-lattice and crossbridge, which are balanced on the thin filament. Note that each thin filament is surrounded by three thick filaments and three thin filaments and each thick filament is surrounded by six thin filaments. Therefore, for a single thin filament, the  55  constitutive equation can be written as follows:  fa =  3 fca + 3 fta , 6  3 τc = 3τt + τz + bφ˙ , 6  (2.21)  (2.22)  where fa is total axial (contractile) force, fca is cross-bridges axial force, and fta is titin axial force. τc , τt , τz are torques from cross-bridges, titin and z-lattice respectively. The last term bφ˙ is a viscosity term, by which we assume that φ˙ is related to the torque difference. b determines how fast the torque difference is compensated by the rotation of the thin filament, which is set to 30 in our simulation. Note that this is not a physical constant and the model assumes the quasi-static situation, which is a valid assumption since the residual force enhancement is known to have no velocity dependence [38]. Given τc , τt and τz , φ˙ is determined and used in Eqn. 2.17. Simulation Results Fig. 2.16 shows the simulation result of isometric contraction at three different lengths. Only one thin and one thick filament, and one titin molecule are shown for clarity. The result shows that the suggested model can simulate the winding behavior as a combined mechanical consequence of N2A attachment to the thin filament and thin filament rotation. Fig. 2.17 shows the simulation result of the activation-lengthening scenario. Since the portion of PEVK bound to the thin filament is not stretched, the strain accumulated in the free portion of PEVK due to the stretch well exceeds that of the isometric contraction at the stretched length (Fig. 2.16, second and third columns). 56  fractional extension  1.2 1.0  1000  1000  900  900  800  800  700  700  600  600  500  500  500  400  400  400  300  300  200  200  100  100  1000  0.8 0.6  900  0.4 0.2  800  0.0  Axial distance (nm)  700  600  0 −10  0  10  0 10 20 30  300  200  100  0 −10  0  10  0 10 20 30  0 −10  10 0  10  30  0  20  Radial distance (nm)  Figure 2.16: Simulation of winding in isometric contraction for different lengths (2400, 2800, 3200nm): Lower and higher grey cylinders represents thin and thick filaments respectively. Thinner black cylinders at the tip of the thick filament show the distal Ig region of the titin. Red dot on the thin filament represents the binding site of the N2A region, the axial distance of which is determined by Eqn. 2.7. Colored strands are PEVK region, where the color represents the fractional extension, which can be converted to the strain using Eqn. 2.18  57  These results suggest that the history dependence of the muscle can be explained by the suggested model of the titin winding. Using the proposed model, we can also simulate the residual force enhancement. We use the same activation-lengthening protocol used by Edman et al. [38]. The goal is to qualitatively simulate the enhancement pattern observed in their study (Fig. 2.18). By simulating the activation-lengthening for different lengths on the descending limb, an overall pattern of the residual force enhancement can be obtained as shown in Fig. 2.19. From the simulation, the pattern of the enhancement reported by Edman et al. can be qualitatively explained (Fig. 2.20). It would be useful to speculate what feature of the mechanism explains the length-dependent variation of the residual force enhancement. According to our model, the axial strain of free PEVK portion depends on two features: 1) the winding pitch angle (θ ) and 2) cross-bridge torque. As the length of muscle increases on the descending limb, the former increases and the latter decreases. Therefore, the sum of these two effects may explain the rise and fall of the residual force enhancement pattern.  2.3.4  Conclusion  In conclusion, we have proposed a new model of muscle mechanics based on the hypothesis of titin winding. The model introduces a new rotational degree of freedom to the contraction mechanics of muscle, which may explain many previously unexplained mechanical properties. Based on the proposed model, a simulation study of residual force enhancement was conducted. From the simulation of an activation-lengthening experiment with the rabbit psoas muscle, the pattern of residual force enhancement can be explained by the proposed model with a mini58  Sarcomere length(nm) 2400 to 2800 fractional extension  1.2 1.0  2800 to 3200  1000  1000  900  900  800  800  0.8 0.6 0.4 0.2 0.0 700  Axial distance (nm)  700  600  600  500  500  400  400  300  300  200  200  100  100  0  0 -10  0  10  30  20  10  0  -10  0 10 0  20 10  30  Radial distance (nm) Figure 2.17: Simulation of winding in isometric contraction followed by stretch: (Left) Isometric contraction at 2400nm and stretched to 2800nm. (Right) Isometric contraction at 2800nm and stretched to 3200nm. Figure legends are the same as Fig. 2.16. 59  Figure 2.18: Residual force enhancement of single fibers from the frog tibialis anterior muscle (from Edman et al. [38]) 0.16 0.14  stress (Nmm-2 )  0.12 0.1 0.08 0.06 0.04 0.02 0 2400  2600  2800  3000  3200  3400  3600  3800  4000  sarcomere length (nm)  Figure 2.19: Simulation of activation-lengthening of the rabbit psoas muscle by winding filament hypothesis: Green baseline represents the active force-length relationship of the descending limb, whereas red and blue baselines represent passive (titin) and total force-length relationship. Each blue branch shows the enhanced force by active lengthening for two different stretches (60nm and 150nm per sarcomere). If blue branches are subtracted by red and green baselines respectively, the enhancement by titin (red branches) and resultant enhancement of the active force (black branches) are obtained.  60  -2  residual force enhancement per cross−section area (N/mm )  0.04  0.03  0.02  0.01  0  −0.01 2400  2600  2800  3000 3200 3400 sarcomere length (nm)  3600  3800  4000  Figure 2.20: Simulated residual force enhancement plotted in the same format as Figure 2.18. mal number of free parameters. Although the result of this simulation is promising, the current model is still at a very early stage of development. First of all, we need direct evidence of titin winding. Even if many experimental observations indirectly support the hypothesis of titin winding, the nano-scale phenomenon of actual titin winding in an active muscle is still not technically observable. On the modeling side, the velocity-dependent properties of winding and its interaction with cross-bridges should be added since, as discussed earlier, the current model assumes a quasi-static situation and therefore cannot simulate transient properties such as the FV relationship and thermodynamics. Another crude assumption of the current model that needs to be revised is the rigid binding of N2A and PEVK to the thin filament. In other words, the current model cannot explain how the bound N2A and PEVK segments become unbound from the thin filament depending on the depletion of Ca2+ or a change in tension. In order to fully simulate calcium-dependent winding and unwinding, as well as possible changes in 61  winding pitch with length change, a more realistic description of binding is required. Finally, the simulation result shows that the contribution of PEVK to the instantaneous elasticity [78] cannot be solely explained by the winding filament model, since the simulated axial strain of free PEVK is too small to explain the observed elastic recoil of the muscle [105]. For this reason, an additional mechanism such as half-sarcomere non-uniformity [135, 168, 170], needs to be integrated to the current model. After building a realistic simulation model that reproduces those basic observations, the next step will be focused on a more detailed analysis of muscle mechanics, such as the rate of energy liberation [68], stiffness variation in isokinetic release and stretch [157], and the storage of elastic energy [105]. It is hoped that all these issues will be explained by the hypothesis of titin winding.  2.4  Conclusion  We have investigated two computational models of muscle mechanics. First, we have shown that the force-length-velocity dynamics of shortening muscle can be well simulated by the Hill-type muscle model, but the formulation of which should be different from the traditional form. Secondly, after discussing several problems of the Hill-type model or cross-bridge model, such as instability and stiffness, we have proposed a winding filament model. The model is based on a new hypothesis that titin winds around the thin filament by cross-bridges generated torque, and therefore is actively engaged in muscle contraction. As a promising example, we have showed that the proposed computational model of titin winding can simulate the residual force enhancement pattern. Since the presented models are focused on different mechanical aspects of the 62  muscle, combining them to a unified model requires substantial future work. The first question will be how the complex interaction between titin and cross-bridges in shortening muscle is related to the relatively simple Hill’s characteristic equation and how it is related to our first result that f-max scaling model better represents the shortening dynamics. Viscous unwinding of titin during active shortening may provide a possible mechanism for the depression of the shortening velocity, which is the main mechanical insight of the f-max scaling model, but more detailed modeling study beyond the conceptual explanation is needed. Although the winding filament model qualitatively presents a possible mechanism responsible for the FV characteristics of shortening muscle, the simulation of this mechanism cannot be done directly since it also requires an alternative computational model of cross-bridge mechanics. Nevertheless, it is hoped that incorporating winding filament model, i.e. adding a new active mechanical element in muscle contraction, can help us to find a better and more succinct model of the cross-bridges mechanics that is not clearly explained even with highly complicated models. Other important directions for future work are to investigate the effects of largescale mechanical behaviors such as sarcomere heterogeneity [21, 170] and force transmission in fiber architecture [6, 73], which are not covered in this thesis.  63  Chapter 3  Modeling Musculoskeletal Systems 3.1  Background  In the previous chapter, we focused on building a realistic computational model of muscle. On the larger scale, the body is made up of muscles, bones and connective tissues such as cartilage, tendons, and ligaments. These parts form a complicated mechanical network of the musculoskeletal system. Unlike man-made machinery, most of these components are highly deformable with large interaction forces, which makes interpreting and simulating their coupled dynamics a challenging task. In this chapter, we develop a framework of building a realistic model for the musculoskeletal system. Similar to muscle mechanics, a simulation-based computational model can be a promising tool for understanding the musculoskeletal system. In previous studies, computational models of the musculoskeletal system [31, 167, 172] have been 64  used in many different contexts and have been helpful in gaining insights into biological motor control. Although there are a wide variety of musculoskeletal models used in the literature, in general these models are created from “the bottom up”: investigators make a large number of independent measurements of physiological parameters (e.g., optimal muscle length, physiological cross-sectional area, moment arms, etc.) which are then used to specify a model of muscle action [16, 20, 71, 88, 93, 94, 121, 175]. Such reductionist approaches are very appealing since they provide a standardized way to construct musculoskeletal models across a wide range of species and measured parameters are usually directly related to well defined muscle properties. Although this approach to developing musculoskeletal models has been generally successful, there are several issues concerning this approach which are often difficult to evaluate. At a basic level, the validity of such models can be difficult to establish: the measurement of muscle parameters is subject to error and this error can result in a range of estimated muscle actions [158, 159, 182]. One critical property of muscles included in such models is their moment arm. In muscles that are represented as line segments, the moment arm of a muscle depends on its origin and insertion. However, many muscles have a broad insertion and/or origin. If a muscle is represented with a single origin and insertion to define its line of action, it is not always obvious based on anatomical dissection where to best place these points. Further, it is not always clear what degree of model complexity is necessary to capture muscle action. For instance, to model especially complex muscles with very broad insertions or origins, investigators often divide the muscle into several distinct regions, modeling each region with an independent insertion and origin which together span the observed muscle anatomy [24, 71]. How to divide these 65  complex muscles into distinct regions and the most appropriate location of their respective insertion and origin points is not usually obvious. The present study develops and evaluates an approach which can potentially help guide and validate the development of musculoskeletal models. We take an integrative approach to modeling the actions of individual muscles, using “top down” empirical characterizations of measured muscle actions to create parameterized muscle models. This approach is an extension of that developed in Loeb et al. [117] and is related to other work attempting to estimate muscle moment arms by direct measurement of muscle action [109, 136, 181]. In particular, we measure in situ the forces and moments produced by stimulation of individual muscles in the rat hindlimb. These forces and moments are a direct consequence of the intrinsic force generating properties of a muscle working through its attachment to the skeleton [117]. Any muscle model attempting to characterize the action of a muscle should therefore be able to predict these forces and moments well. Further, because these measurements are closely linked to the properties of the stimulated muscle, it should be possible to estimate the parameters for a muscle model (i.e. its anatomical origin and insertion and its intrinsic force generating properties) based on the forces and moments observed experimentally. The present study uses this approach, using the experimentally observed forces and moments to create a parameterized muscle model to capture the action of muscles. We apply this approach to two different muscles in the rat hindlimb, gracilis posticus and biceps femoris [56]. Gracilis posticus is a relatively simple, bandlike muscle with well defined origin and insertion. We expected that this muscle should be well described with a model consisting of a single insertion and origin. In contrast, biceps femoris is a complex muscle with an insertion distributed across 66  the tibia from the knee to the ankle [51]. Given this complex insertion, it is not obvious where one should place the insertion of this muscle on the tibia, or whether multiple insertions might be required to describe the action of this muscle. We expected that this muscle should be poorly described by a single insertion and origin and would require a more complex model to capture its action. These two muscles, representing a range from simple to complex muscle architectures, thus allow us to evaluate the utility of our approach for the development and validation of muscle models.  3.2 3.2.1  Methods Experimental Preparation  We report the data obtained from five adult female rats (250-300 g). All procedures were approved by the IACUC at Northwestern University. Rats were anesthetized with ketamine/xylazine (80/10 g/kg, i.p.). Metal posts were attached to three points on the sacrum using bone screws and cement: one each on the rostral and caudal ends of the contralateral (right) pelvis and another on the rostral end of the ipsilateral pelvis (Fig. 3.1A). A threaded attachment was then cemented to bone screws attached at the distal tibia so that it was approximately perpendicular to the tibia. The anterior head of gracilis was severed and gently peeled back to expose the nerve innervating the posterior head of gracilis. The nerve was gently freed from the surrounding tissue and a small cuff electrode with bipolar exposed leads was placed around the nerve, insulated with Vaseline, and secured to surrounding tissue using cyanoacrylic glue. To prepare biceps femoris (BF) for stimulation, the anterior head of biceps femoris (BFa) was first severed close to its origin on the  67  A  Y  Y X  Z  posts in pelvis  transducer attachment  Z  B  Y  sagittal plane (X-Y)  X  hip  hip  origin (X,Y,Z)  fem ura  th pa  ur  mu sc le  fem  W  la  xis  pori  pk  V  pin  s  pa  fk  fm  knee  pm  U tibia insertion (U,V,W)  -fibu  la  ankle  Figure 3.1: (A) Experimental setup. The pelvis was immobilized by three posts; two in the rostral and the caudal ends of the contralateral pelvis and one in the rostral end at the ipsilateral pelvis. A six dimensional force and moment transducer was secured to an attachment cemented to the distal tibia on the left hindlimb. In this configuration, the X-YZ axes are aligned with the rostral-caudal, dorsal-ventral and mediallateral direction respectively. (B) Kinematic structure and coordinate frames defined for the hindlimb. The origin of the muscle is defined as a point (X-Y-Z) in a coordinate frame with an origin at the hip and the insertion is defined as a point (U-V-W) in a coordinate frame fixed to the tibia-fibula with an origin at the knee. Skeletal geometry is modeled as a two-link serial chain, where the hip joint is modeled as a ball joint (X-YZ rotation) fixed to the ground (pelvis) and the knee joint is modeled as a universal joint (U-W rotation). The overlaid skeleton is shown to help visualize the coordinate frames. Fig. 1A is modified with permission from Tresch and Bizzi [176] . 68  tibia -fib  ula  fa ma  vertebrae and partially removed by gently cutting its adhesion to the posterior portion of biceps (BFp). BFp here refers to the entire region of BF which has an origin on the ischium and inserts along the extent of the tibia (see Discussion). We then exposed the branch of the tibial nerve innervating BFp by cutting gluteus maximus and caudofemoralis, then peeling back caudofemoralis to expose the nerve. The branch of the nerve to semimembranosus was cut, the remaining nerve freed from surrounding tissue, placed in a cuff electrode, insulated with Vaseline, and secured to the surrounding tissue. We then separated semitendinosus from biceps femoris, cutting through the connections between the two muscles and severing semitendinosus distal to the fusion of its medial and lateral heads. Exposed muscles were covered with moist gauze. Finally, we performed a dorsal laminectomy to expose the lumbar spinal cord and dorsal roots. The L4 and L5 ipsilateral dorsal roots could be visually identified by their size and, after opening the dura and applying a few drops of lidocaine to the cord surface, the ipsilateral dorsal roots from L1 to L6 were severed. We usually also cut a dorsal root rostral and caudal to the identified roots in order to ensure full deafferentation of the limb. The deafferentation was confirmed during the experiments by observing a lack of response to strong pinch of the ipsilateral foot while observing a response to the same pinch applied to the contralateral foot used to assess anesthetic state. The animal was then transferred and positioned in the experimental setup (Fig. 3.1A). The posts cemented to the pelvis were attached to magnetic stands anchored to the table surface so that the pelvis was approximately level in medial-lateral and anterior-posterior directions. The threaded rod in the distal tibia was attached to a 6 axis force transducer (ATI; Mini40). The transducer was mounted on a series of stages which allowed it to be moved along the dorsoventral, mediolateral, and 69  rostrocaudal axes of the rat. Once the limb was attached to the transducer the only possible movement at the ankle was a rotation around the mediolateral axis. This rotation could be prevented or permitted by set screws. The ankle was positioned so that it was approximately in the same mediolateral plane as the hip center. The temperature of the rat was monitored with a rectal thermometer and maintained between 37 and 39◦ C using a heat lamp placed above the animal. Fluids were given to the animal immediately after the initial surgical preparation (Dextrose 3cc subcutaneous) and in between collections of force field data sets (approximately every 3 hours).  3.2.2  Experimental Protocol and Data Collection  We first found recruitment curves for individual muscles, varying the current level and measuring the change in evoked force. Stimulation trains (75Hz, 0.2ms pulses, 500-1000ms duration) were applied through the implanted cuff electrodes in order to evoke a muscle contraction with the limb fixed at a specific configuration in the workspace. The forces and moments measured by the transducer were collected (2000Hz) using custom software in Labview for a period of 1s before and 3s after the stimulation train onset. Stimulation trains were applied at a fixed interval of 1 minute. Current was increased gradually from the threshold response until the evoked peak magnitude reached a plateau level. We measured the response for several current levels above the plateau level as well, making sure that the evoked force was not altered, in order to ensure lack of current spread to adjacent muscles. From these recruitment curves, we chose a current level which evoked a maximal activation of the muscle (between 0.06 and 0.45mA). Once the current level was established, we then measured how the evoked 70  forces and moments varied across the hindlimb workspace. The limb was fixed in a configuration, the muscle stimulated and the evoked forces and moments measured at the ankle with the attached force transducer. Two types of trials were collected for each muscle at a single limb configuration. In one trial, all degrees of freedom at the ankle were restrained (by securing the set screws; referred to as “restrained” trials). In the other trial, rotation about the mediolateral axis of the ankle was permitted (by releasing the set screws; referred to as “released” trials). By visual inspection, any such rotations were generally very small during evoked responses. Stimulation trains to any muscle were applied at an interval of 1 minute. After collecting the evoked responses at a single limb configuration, the ankle was moved to a new position. We took care not to move the limb to the extremes of the workspace to avoid kinematic singularities or overstretching muscles. Stimulation was then applied to each muscle at the new position, for both restrained and released conditions. 19-24 limb configurations were measured for each animal in order to span a range of limb configurations throughout the hindlimb workspace. Periodically through the collection of these responses the limb was returned to the original configuration and stimulation repeated to evaluate the stability of the evoked responses. At the end of the experiment, the position of the hip center of rotation was measured and the limb was dissected to find the distance from the hip to the knee center of rotation and the distance between the knee center of rotation and the rod attached to the foot. The mediolateral distance from the origin of the force transducer to the center of the foot was also measured for each animal in order to set the tool transform for the collected forces and moments.  71  3.2.3  Data Analysis and Modelling  We examined only the active force produced by each muscle, obtained by subtracting the averaged forces and moments measured prior to the onset of each stimulation train. We then found the peak force produced by each muscle and extracted the 6 dimensional vector of peak forces and moments at each position. The set of these 6 dimensional vectors measured at each limb configuration comprises a “force/moment field” similar to the force fields which have been used previously [52, 117]. As described above, these fields are the result of the intrinsic force generating properties of each muscle acting through the attachment of the muscle to the skeleton. In the approach described here we attempt to find a model of muscle action which is capable of explaining the observed muscle actions. Parameterized Muscle Model We characterized the muscle by two parameters: maximal force (peak tetanic tension) and optimal muscle length. Note that we did not include an extra parameter for tendon length in these models. Although this parameter is in general necessary to account for the difference between musculotendon length and muscle fascicle length, including this extra parameter did not substantially improve the quality of the fits described here. Moreover, as long as the muscles acts predominantly upon the ascending portion of the force length functions (see Fig. 3.8), the parameters of maximal force, optimum length, and tendon slack length have redundant effects on the predicted muscle force for the isometric measurements made here. For fully activated muscle, we approximated the relationship between muscle length and force magnitude with a Gaussian curve with s.d. equal to 0.2 times muscle length:  72  fm = fmax exp −  (l/lo − 1)2 2(0.2)2  (3.1)  Where l is muscle length, fmax is the maximal force and l0 is the optimal muscle length. The force produced by the muscle acts through its attachment points on the skeleton. In this study we considered each muscle to have a point origin and insertion. Each point could be an arbitrary position in three dimensional space, with the insertion described in a local frame attached to the tibia-fibula. Thus, the action of each muscle was characterized by 8 parameters, three for the origin, three for the insertion, the optimal length, and the maximal force. Kinematic Model of the Rat Hindlimb As shown in Fig. 3.1B, we modeled the rat hindlimb as a two-link serial chain consisting of the femur and tibia-fibula. The foot was omitted from analyses since we were concerned with muscles acting on the proximal joints. The hip was modeled as a ball joint and the knee was modeled as a universal joint with motion in flexion/extension in the sagittal plane and motion in pronation/supination of the tibia-fibula. There were thus 5 degrees of freedom in the limb. Constant centers of rotation were used for both the hip and knee joints. Relating Muscle Force to Forces and Moments at the Ankle To predict forces and moments measured at the ankle for a given muscle force, we used standard methods from robotics based on Jacobians to translate between coordinate frames [126]. For six dimensional space, the structure of Jacobian for n  73  degrees of freesdom limb is represented as follows:     v1 · · · J= w1 · · ·  vn  6×n ∈R , wn  (3.2)  where vi and wi are three dimensional column vectors representing the screw parameters of the ith degree of freedom. For rotary joint, wi is the axis of rotation and vi is represented as qi × wi , where qi is an arbitrary point on the axis of rotation. Using Jacobian transpose, we can relate the force produced by the muscle to the ankle forces and moments according to the following equation:  τ = JaT ψa = JmT fm ,  (3.3)  where fm is a (3 × 1) vector of forces produced by the muscle, ψa is the (6 × 1) vector of forces and moments measured at the ankle, Ja is the (6 × 5) Jacobian relating muscle force/moments measured at the ankle to joint torques, and Jm is the (3 × 5) Jacobian relating muscle force to joint torques (in this case, we only use the upper three rows of the Eqn. 3.2), τ (a 5 × 1 vector).In this equation, both sides of the equation are equal to joint torques. Ja and Jm were calculated using standard formulations of the Jacobian based on a rigid body analysis. The direction of fm was calculated from the orign and insertion of each muscle and its magnitude was calculated by Eqn. 3.1. Ja , Jm and fm vary with the configuration of the limb. If one considers the configuration of these experiments, it can be seen that there is a mismatch in the number of degrees of freedom within the hindlimb (5) and the measured forces and moments (6). Because of this, when the ankle was completely restrained by fixing the rotation around the mediolateral axis, the limb 74  was overconstrained by the attachment to the force transducer. In such a situation, translating between intrinsic forces produced by the muscle and the forces and moments measured by the force transducer is not a well posed problem: i.e. in general, there is not a unique mapping between the two sets of forces although the physical system has a unique mapping, dictated by the intrinsic stiffnesses of the force transducer and hindlimb [127]. Mathematically, this problem can be described as finding ψa in Eqn. 3.3 for a given Ja , Jm , and fm . Note that because Ja is not square it cannot be inverted, so that although transforming ankle forces and moments into joint torques is well defined, the mapping from joint torques to ankle forces and moments is not. In order to translate between intrinsic muscle force and endpoint forces/moments, we therefore formed the Moore-Penrose (M-P) pseudoinverse [54] of Ja , and used it to predict the forces and moments at the ankle which result from contraction of the muscle.  ψa = Ja (JaT Ja )−1 τ = Ja (JaT Ja )−1 JmT fm .  (3.4)  The M-P pseudoinverse is similar to the standard matrix inverse in that JaT ψa gives τ back, as can be easily seen from Eqn. 3.4. Of the many possible values of ψa , Eqn. 3.4 picks the one with the smallest magnitude. However, other pseudoinverses are also possible [127]. We therefore compared this approach to a complementary one in which a pseudoinverse was not necessary. In the second approach, referred to here as the “released” method, we used the trials in which the rotational degree of freedom around the mediolateral axis was released (the “released” trials described above). If this degree of freedom is  75  released and we assume that as a consequence there are no moments about this mediolateral axis, then the 5 degrees of freedom of the limb are matched with the 5 degrees of freedom measured by the force transducer. In this situation, ψa is reduced to a 5 dimensional vector ψ˜ a (3 forces and 2 moments) and we can define a reduced (5 × 5) Jacobian J˜a by removing the row from Ja corresponding to the rotation along the released axis. Since J˜a is square and thus invertible, ψ˜ a can be uniquely determined as follows;  ψ˜ a = (J˜aT )−1 τ = (J˜aT )−1 JmT fm .  (3.5)  Note that by releasing this degree of freedom the forces and moments acting at the ankle can, in principle, differ from those measured with the limb fully restrained. However, we found that the forces and moments were only minimally altered when this degree of freedom was released (see Results). Parameteric Optimization We found the set of muscle model parameters which best predicted the force/moment fields measured at the ankle, minimizing the error between observed and predicted force/moment fields. If forces and moments were measured at N limb configurations, the optimization problem can be defined as: N  min ∑ ψi − ψ¯ i (d) d i=1  2  ,  (3.6)  where d are the 8 muscle model parameters and ψ¯ i (d) is the ankle force and moments predicted for the ith configuration. ψ¯ i (d) was the predicted action of the muscle from either the Jacobian or static analyses described in the previous section,  76  and ψi was the data taken from either the restrained or released trials. Units of force and moment were adjusted to have similar scales. The set of parameters which minimized this error was found using lsqnonlin in Matlab. The optimization routine was started with reasonable guesses for the initial parameters, with the origin placed caudal to the hip center and the insertion placed near the tibia-fibula. We also tried different initial guesses, and found that the result was robust to the selection of initial guesses. A concern in fitting any model is whether the identified model overfits the data, explaining noise in the data rather than the structural variation. We therefore performed an n-fold cross validation. Each force field was divided into 5 sets of positions, chosen randomly. One set of positions was set aside and the remaining 4 sets were used to find model parameters, using the optimizations described above. We then evaluated the prediction of the identified model over the set of positions which were set aside. This procedure was repeated for each of the 5 sets of positions. Evaluating Fit Quality To quantify the quality of the fit to the observed data, we used a centered R2 statistic, computed as 1 minus the ratio between the residual error variance in the model and the total variance of the data. An R2 value near 1 indicates that there is minimal residual error, whereas an R2 near 0 indicates that the model performs the same as the mean of the data. Note that this is a conservative measure of fit quality in the present case, since predicting the mean force produced by the muscle is non-trivial but is a result of the specific origins and insertions found by the optimization.  77  Comparison to Existing Muscle Models Johnson et al. [88] report a comprehensive musculo-skeletal model for the rat hindlimb based on detailed anatomical measurements. We used these data plus the muscle architecture data of Eng et al. [39] as the basis of a forward model to compute the force/moment fields in the present study. Johnson et al.’s variable joint center data, femur and tibia lengths, and muscle origin and insertion locations were incorporated in the model. This model will be referred to in this manuscript as the Johnson parameter model. Note that Johnson et al. have updated the model in OpenSim [33] to include wrapping surfaces, made slight changes in some parameters, and the most recent model version considers the knee to move in only flexion extension. In the predictions shown here we used the model parameters and structure as described in Johnson et al. As with our method above, the model has 5 degrees of freedom: X,Y,Z rotations for the hip joint and U,W rotations for the knee joint (see Fig. 3.1). To approximate the experimental conditions, we assumed that the knee was in the same mediolateral plane as the ankle (see Figs). Using the femur and tibia lengths from Johnson et al., we found the corresponding joint angles, allowing for the variable joint centers described in Johnson et al. These 5 angles were used to position the skeleton and determine the origin and insertion of the BFp and GRp. Ankle force and moments were calculated using Eqn. 3.4. Ja and Jm were determined numerically from the Johnson parameters for each set of joint angles. Because in the Johnson et al. manuscript joint centers are a function of joint angle, Ja and Jm as calculated here are slightly different compared to the analytic calculations used above. The magnitude of fm was determined using a force-length relationship based on the data of Eng et al. [39]. The shape of the  78  function was approximated using a idealized sarcomere force-length function [55] corrected for actin length in the rat [146]. The width of the function was scaled by muscle fiber length. The optimal length was set to muscle length. The maximal force is based on PCSA multiplied by 22 N/cm2 .  3.3 3.3.1  Results Forces and Moments Evoked by In Situ Stimulation of Rat Hindlimb Muscles  Fig. 3.2A shows an example of the forces and moments evoked from stimulation of biceps femoris at a single limb configuration. As seen in the trial, the evoked response quickly reached a plateau level which was sustained for the duration of the stimulation train then returned to the prestimulation baseline levels. Fig. 3.2B shows an example of the variation in the evoked force magnitude with the current level used to activate the nerve. In this example, the magnitude reached a plateau level at 0.13 mA. Further increases of current past this level did not result in increases in the evoked force magnitude, suggesting that the muscle was maximally active. In Fig. 3.2C we show the force and moment vectors in the XY plane taken at the time of peak response for the trial shown in Fig. 3.2A. We separate the plots of forces and moments in order to aid their visualization, although in the data analyses and development of muscle models described below they were combined in a single 6 dimensional vector. We also show the force and moment vectors evoked on repeated trials during the collection of a single force field in Fig. 3.2C. The vectors were collected at intervals of 4, 6, 6, and 7 minutes. The similarity of these vectors demonstrates the consistency of muscle stimulation in these experiments, although there is also fatigue over time. Taken together, the results shown in Fig. 3.2 (see 79  A  X  500 msec  FORCE  3 N cm  MOMENT  1N  X Y  Y  Z  Z STIMULATION  STIMULATION  B  C  3  MOMENT  FORCE (N)  2.5 2 1.5  FORCE  1 0.5 0.06  0.08  0.1  0.12  0.14  CURRENT (mA)  Figure 3.2: Force and moments evoked by in situ muscle stimulation in the rat hindlimb: (A) Recorded forces and moments evoked by BFp stimulation for one trial. Force and moment quickly reached plateau levels at the onset of the stimulation and returned to baseline levels after the end of the stimulation. (B) Magnitudes of the evoked force of biceps femoris across a range of current levels. As the current level increases, the force reached a plateau level (around 0.13 mA) and did not increase further, suggesting that the muscle was maximally activated after 0.13 mA. (C) Variation of evoked force and moment by biceps femoris in the same limb configuration across 5 repeated trials, separated by 4, 6, 6 and 7 minutes. 80  also Figs 3.3,3.4) demonstrate the reliability of the experimental preparation and collected data used in the present experiments to characterize muscle actions.  3.3.2  Force and Moment Fields Produced by Stimulation of Gracilis Posticus and Biceps Femoris  Fig. 3.3 shows the force and moment fields evoked by stimulation of GRp for each of the 5 animals, with the ankle movement fully restrained. The equivalent fields for BFp are shown in Fig. 3.4. GRp produced a set of forces which drove the limb caudally, dorsally, and medially. Fig. 3.3 shows that the force and moment fields measured for GRp were very consistent across animals. Fig. 3.5A shows an example of two force and moment fields measured for one animal either with the ankle fully restrained or with rotation around the mediolateral axis released. As can be seen in the figure, the moments around the mediolateral axis were minimal when the rotation was permitted while other degrees of freedom were only minimally affected. Although the forces and moments measured for GRp were generally consistent with what might be expected based on the anatomy of this muscle (see below), the responses measured for BFp were less expected. Although Fig. 3.4 shows that the responses evoked by BFp stimulation were similar across animals, the evoked forces were directed caudally in extension with only a minimal dorsal component. Such a caudally directed force suggests that BFp acted predominantly as a hip extensor with minimal knee flexion. Fig. 3.5B shows an example of force and moment fields collected for BFp with mediolateral rotation either released or restrained, showing that the evoked forces and moments were minimally affected by releasing this degree of freedom.  81  seen from caudal direction (YZ plane)  seen from lateral direction (XY plane)  Y  Y Z  X  ANIMAL 5  ANIMAL 4  ANIMAL 3  ANIMAL 2  ANIMAL 1  hip center of rotation  MOMENT (1  FORCE (1 N)  Figure 3.3: Force and moment fields evoked by in situ stimulation of GRp for each animal: First two columns show the recorded force (filled arrows) and moment (open arrows) field as seen from the lateral direction (XY plane, see Fig. 3.1A), while the last two columns show the field as seen from the caudal direction (YZ plane, see Fig. 3.1A). The position of the hip in each plot is indicated by the filled circle.  82  seen from caudal direction (YZ plane)  seen from lateral direction (XY plane)  Y  Y Z  X  ANIMAL 5  ANIMAL 4  ANIMAL 3  ANIMAL 2  ANIMAL 1  hip center of rotation  MOMENT (5  FORCE (5 N)  Figure 3.4: Force and moment fields evoked by in situ stimulation of BFp. Conventions are the same as in Fig. 3.3.  83  A  hip center of rotation  GRp  GRp  force in restrained trial force in released trial  moment in restrained trial moment in released trial  B BFp  BFp  XY plane  XY plane  YZ plane  YZ plane  Figure 3.5: Comparison of force and moment fields measured in “restrained” trials, in which all movement at the ankle was constrained, and “released” trials in which rotation about the z axis was allowed: (A) Comparison of GRp data collected from animal 3. (Gray arrows) Measured force and moment fields in “restrained” trial, which is same as data shown in Fig. 3.3. (Black arrows) Measured force and moment fields in “released” trial. Note that Z moments in the “released” trial are close to zero, since the ankle is allowed to rotate around Z axis. The other forces/moments were only minimally altered between the two restraint conditions. (B) Comparision of BFp data collected from animal 3. Conventions are the same as A. Note that the force/moment fields for BFp were nearly identical for released and restrained trials. That is, although both trials are plotted in each figure, at most positions they are indistinguishable.  84  3.3.3  Using Measured Force and Moment Fields to Create Muscle Models  The force and moment fields illustrated in Figs 3.3 and 3.4 reflect the intrinsic force generating properties of each muscle and the attachment of the muscle to the skeleton. We therefore attempted to use the observed muscle actions to develop a parameterized muscle model for GRp and BFp. Each model was described by 8 parameters: the optimal muscle length, the maximal force, and the origin and insertion (each a point in 3 dimensional space) of each muscle. The parameters which minimized the error between the observed and predicted forces/moment fields were found by optimization. We applied this optimization to the data obtained with the ankle either fully restrained or released, either using the Jacobian with M-P pseudoinverse (restrained analysis) or using the invertible Jacobian without z-moments (released analysis). Fig. 3.6 shows the predictions of the model found for GRp using the “restrained” method for one animal. In each plot, the thick dark line indicates the predicted path of the muscle from its origin to its insertion. As can be seen in the plots, the identified model does a very good job of capturing both the forces and moments produced by stimulation of GRp. The one exception to this good fit is for the z-moments: the observed response produced a negative z-moment, whereas the model predicted a positive z-moment (see Discussion). In addition to this generally good fit to the data, the origins and insertions identified by the model are also consistent with the expected anatomy of GRp, with an origin near the midline of the pelvis and caudal to the hip, and with an insertion on the tibia partially towards the ankle. Fig. 3.8A (results from the restrained method are shown in empty symbols) shows the parameters for the origins and insertions found for GRp in each animal. 85  GRp  measured moment predicted moment  GR p  FE MU R  measured force predicted force  TIB  IA-  FIB U  LA  XY plane  YZ plane  YZ plane  XY plane  Figure 3.6: Example prediction of the force and moment fields evoked by GRp stimulation: (Gray arrows) measured force and moments collected from animal 1 (as shown in Fig. 3.3). (Black arrows) force and moments predicted from the model identified using the “restrained” method. The model prediction shows good fit in both forces and moments, but note discrepancy of moments in YZ plane. The estimated muscle path connecting origin and insertion point are shown in a thick black line in all plots for a single limb configuration. Detailed parameter values are shown in Fig. 3.8.  measured force predicted force  measured moment predicted moment  BF p  FE  MU  R  BFp  TIB  IA-  FIB  UL  A  XY plane  XY plane  YZ plane  YZ plane  Figure 3.7: Example prediction of force and moment fields evoked by BFp stimulation. Conventions are the same as Fig. 3.6.  86  A  2  empty - ‘restrained’ method filled - ‘released’ method  animal 1 animal 2 animal 3 animal 4 animal 5  1.5  0.5 0  1  −0.5 −1  C  0.5 0 −0.5  x  y  z  u  v  −1  w  D  2.5  2  2  1.5  POSITION (cm)  POSITION (cm)  1  B  animal 1 animal 2 animal 3 animal 4 animal 5  x  y  z  u  v  w  20 18 16  FORCE (N)  FORCE (N)  14  1.5  1  12 10 8 6  0.5  4 2  0 1.5  2  2.5  3  3.5  4  4.5  LENGTH (cm)  0 1.5  2  2.5  3  3.5  4  4.5  LENGTH (cm)  Figure 3.8: Estimated parameters for GRp (A) and BFp (B) models for each animal: Empty and filled symbols compare parameter values estimated by “restrained” and “released” method, respectively. First three parameters are position of the muscle origin represented in X-Y-Z coordinate and last three represent muscle insertion in U-V-W coordinate (coordinates are defined in Fig. 3.1). Force-length relationships of GRp (C) and BFp (D) estimated by the model defined in Eqn. 3.1. Line shows the parameterized force length function and each symbol shows a point on the curve where the force was actually evaluated.  87  As can be seen in the figure, these parameters were very similar across animals, as might be expected based on the similarity of the force/moment fields shown in Fig. 3.3. In Fig. 3.8C we show the predicted force length functions for each animal, suggesting that in each animal GRp was acting on the ascending portion of its force length function. Fig. 3.7 shows the results of the model found for BFp using the “restrained” method for the same animal shown in Fig. 3.3. As was the case for GRp, the identified model did a very good job of capturing the observed actions of the muscle for both forces and moments. The origin of the model is also consistent with the expected anatomy of the muscle which arises from the ischium caudal to the hip. Note the more lateral and dorsal origin of BFp as compared to GRp. The insertion identified in the model was very close to the knee, consistent with the mainly hip extensor action suggested by the plots of Fig. 3.4. Fig. 3.8B (empty symbols) shows the origins and insertions found for the BFp model in each animal, again showing the consistency of the estimated parameters across animals. Fig. 3.8D shows that, similar to GRp, the model predicts that BFp acts along the ascending portion of its force length function. As described above, we also identified models based on the data sets with “released” trials, in which the rotation around the mediolateral axis was permitted. These models were found using the invertible form of the Jacobian. In Fig. 3.8A,B (filled symbols) we show the parameters found using these models from the released analysis for GRp and BFp, respectively. Comparison of empty and filled symbols in plots shows that, in general, the model parameters found using the two approaches were very similar to one another. The models using the “released” method consistently had high R2 values for 88  both BFp (mean=0.88, s.d.=0.02) and GRp (mean=0.89, s.d.=0.04), consistent with the good fit quality seen in the plots. The models found using the “restrained” method for BFp also had high R2 values (mean=0.84, s.d.=0.06). For the GRp restrained analysis, the R2 values were considerably lower (mean=0.27, s.d.=0.29), primarily because of the consistent offset in the z-moment as illustrated in Fig. 3.3. The R2 values for GRp increased considerably if this z-moment was ignored (mean=0.85, s.d.=0.06). Note that this is a centered R2 measure which considers only variations around the mean value in each force/moment dimension. Finally, we examined whether the models identified above were overfitting the observed data. We performed an n-fold cross validation procedure, fitting the data on 80% of the data then evaluating the fit on the remaining 20% of the data. We found that the quality of the fits was only minimally reduced when evaluating such cross validation fits, considering either the force/moment fields predicted using the “released” method for GRp (mean=0.85, s.d.=0.03) or BFp (mean=0.85, s.d.=0.06) or using the “restrained” method for GRp (mean=0.81, s.d.=0.06, ignoring the zmoments) and BFp (mean=0.85, s.d.=0.03, ignoring the z-moments). This lack of overfitting suggests that the 6 dimensional data vectors measured at multiple positions across the workspace place strong constraints on the 8 parameters used in the model. Taken together, the results show that simple muscle models consisting of point origins and point insertions are able to capture the actions of individual muscles in the rat hindlimb very well using either the “restrained” or the “released” procedures. Further, the models identified by the two procedures and identified across animals were very similar to one another.  89  3.3.4  Comparison to Existing Musculoskeletal Models of the Rat Hindlimb  A  GRp  moment (1 N cm)  GR p  GR p  FE  FE  MU R  MU R  force (1 N)  TIB I  TIB  A-F I  BU  IA-  LA  B  FIB U  LA  BFp  BF p  FE  MU  R  moment (5 N cm)  BF p  FE  MU  R  force (5 N)  TIB  TIB  IAFI  IA-  BU  LA  XY plane  FIB U  LA  XY plane  YZ plane  YZ plane  Figure 3.9: Prediction of force and moment fields of GRp (A) and BFp (B) using anatomical parameters taken Johnson et al. [88]. Conventions are the same as those in Fig. 3.6. The results are shown both in lateral (XY plane) and caudal (YZ plane) direction. Fig. 3.9 shows the force and moment fields predicted from the Johnson parameter model. As can be seen from the figures, the fields predicted for GRp from the published parameters were very similar to those found in the present study by direct stimulation of GRp (Fig. 3.3). On the other hand, the force/moment field predicted for BFp is very different than the one observed experimentally from di90  rect stimulation (Fig. 3.4). In particular, there is a much more prominent knee flexion action predicted for BFp using the published parameters. This knee flexion action is a consequence of the insertion of BFp being placed distally along the tibia, reflecting the anatomically distributed insertion of BFp across the tibia. The results of this comparison suggest that the mechanical action of BFp as measured here is better predicted by a more proximal insertion, such as those shown in Fig. 3.4 based on direct in situ measurements of BFp action, than by a more distal insertion, based on measurements of BFp anatomy.  3.4 3.4.1  Discussion Characterization of Muscle Actions by In Situ Stimulation  As has been suggested in previous work [52, 117], the force and moment fields measured here give a succinct representation of muscle actions and their potential contributions to motor control. For instance, although standard characterizations might describe both GRp and BFp as hip extensor and knee flexors [56], the force and moment fields shown in Figs 3.3, 3.4 illustrate that the actions of these muscles are clearly distinct, both in terms of their actions in the sagittal plane but also along other degrees of freedom. The description here involving 6 dimensional forces and moments provides an especially rich description of muscle actions. It is likely that similar characterizations for other hindlimb muscles [106, 201] can help to highlight the unique contributions of individual muscles and provide insights into the coordination strategies implemented by the nervous system. The force and moment fields measured in response to GRp stimulation were consistent with the actions expected from anatomical inspection of this muscle and  91  were similar to those predicted based on parameters taken from a previous study using anatomical dissections [88]. This consistency helps validate the force and moment fields found here and the general approach taken in this study, in addition to helping validate the findings of that previous study. The fields measured in response to BFp stimulation, on the other hand, were more surprising. They reflected a predominantly hip extension action, with a relatively small knee flexion contribution. This dominant hip extension action was also confirmed by the proximal insertion of BFp estimated by our optimization procedures. Straightforward anatomical inspection of BFp suggests that it should produce a more significant knee flexion action than that measured here [88]. In fact, the fields predicted by parameters estimated in a previous study using anatomical inspection were clearly different from the ones measured here, reflecting the more distal insertion of BFp estimated in that study. It is unlikely that the mainly hip extension action of BFp observed here reflected an incomplete activation of all regions of BFp since with increasing stimulation strength the evoked force magnitude reached a constant level which did not change with further increases in strength (Fig. 3.2B). This difference between the insertion of BFp estimated from anatomical dissection and the insertion estimated from measuring its action illustrates the utility of the approach described here in helping validate and refine musculoskeletal models. Closer inspection of the structure of BFp might provide an explanation for how this predominantly hip extensor action is consistent with the anatomy of BFp. Although the insertion of BFp appears uniformly distributed along the tibia from the knee towards the ankle, reflecting the muscle to expose the medial surface reveals a well defined tendon near the knee into which a large portion of the muscle fibers in BFp insert. BFp is also thickest proximally and tapers more distally. Fig. 3.10 shows the details. Both 92  Figure 3.10: Tendon structure in the insertion area of BFp, seen from the lateral (left) and the medial surface (right). of these features make it likely that this proximal insertion dominates the action of BFp when the entire muscle is stimulated. This division between a proximal region of BFp with a well defined tendon and a more distal region of the muscle with a distributed insertion along the tibia has apparently been recognized in Greene [56], who labeled the more proximal region biceps femoris posterior and the more diffuse distal region biceps femoris accessory. However, the exact division between these two regions is not obvious from simple inspection [128, 134]. In the cat, BFp has been suggested to be divided into three compartments defined by separate innervation territories [23, 26] and the divisions described here in the rat might similarly correspond to different compartments. Further experiments stimulating individual nerve branches to BFp or partially cutting its insertion might help clarify the actions of different regions of BFp. It should also be noted that in order to isolate the muscles stimulated here we had to perform extensive dissections of the hindlimb. These dissections obviously interrupted connective tissue interactions between muscles, for instance between BFp and semitendinosus, BFa and BFp, and GRa and GRp, all of which have  93  close adhesions. Similarly, in order to place bone screws on the distal tibia we interrupted the fascial sheath of the lateral limb. As has been described in previous work, these connective tissue interactions can in some circumstances significantly affect the forces produced by stimulation of individual muscles [74, 81], but these interactions may not be significant in physiological conditions [118, 156]. Future experiments might measure the force and moment fields with and without such connective tissue interactions to gain a better understanding of the consequences of these nontraditional pathways of force transmission.  3.4.2  Developing Musculoskeletal Muscle Models Based on In Situ Measurements of Muscle Actions  Beyond its use in validating and refining musculoskeletal models, we also demonstrated that the approach described here can be used directly to specify the parameters of a muscle model. This approach is related to previous work in which moment arms were estimated at single postures based on forces applied to tendons in cadavers [109, 136, 158, 181, 182] and in which planar force fields were used to find functional approximations of moment arm variations across the workspace [117]. The present study extends that work to use the forces and moments to specify a simple anatomical model of muscle action. We showed that a simple model of muscle action, based on a single origin and insertion along with the optimal length and force of the muscle, was able to describe the actions of both GRp and BFp very well. The good fits for GRp with this model were consistent with the relatively circumscribed origin and insertion of this muscle, suggesting the validity of the approach used here. However, we were surprised by the ability of this simple model with a single  94  origin and insertion to explain the actions of the complex muscle BFp. We expected this muscle to show the limitations of this simple model structure, showing how the measurements from direct muscle stimulation could be used to justify the creation of distinct subregions and to specify their placement. As discussed above, it might be that this good fit reflected the dominant action of the more proximal tendon of BFp, allowing it to be well described with a single insertion. Whether other complex muscles can be as well described with a single origin and insertion will require further experiments. Ultimately, whether a complex muscle such as BFp should be represented by a single insertion or multiple insertions and where those insertions should be placed depends not only on being able to explain the mechanical action of the muscle but also on being able to explain the activation of the muscle by the nervous system. For example, even though the present experiments suggest that we could explain the mechanical action of the entirety of BFp by a single insertion point, it is likely the nervous system does not activate BFp homogenously. Previous work in the adult cat [25, 143] and our own work in the neonatal rat [98] have shown that regions of BFp are activated differentially in behaviors such as locomotion and posture. Obviously, the mechanical consequences of this differential activation cannot be explained using a musculoskeletal model of BFp with a single origin and insertion. Such considerations emphasize the iterative process of refinement and elaboration involved in creating musculoskeletal models, requiring a combination of both biomechanical and neurophysiological studies. The results of this study were mainly focused on estimating the origins and insertions of GRp and BFp based on the structure of the measured force/moment fields. These origins and insertions determine how the direction of the 6 dimen95  sional vector measured at the ankle varies across limb configurations, which is a critical aspect of all musculoskeletal models. Our models also included estimates of the force length properties of each muscle, as defined by the maximal force and optimal muscle length parameters. These properties influence the magnitude of the vector at each limb configuration but do not affect the direction of the predicted vector. The model of force generation used in this study only included parameters defining the maximal force and optimal length of each muscle and did not include a parameter for the resting tendon length [203]. Because we only made isometric measurements in these experiments and because these muscles are likely acting mainly on the ascending portion of their force/length functions (see Fig. 3.9 and from predictions from Johnson model parameters), including this extra parameter is not expected to increase the fit quality of the identified models, consistent with the very good fits shown here. For other muscles, such as distal muscles with long tendons or muscles operating on both the ascending and descending portion of the force length function, it should be possible to estimate the tendon length and increase the quality of the observed fits by including this parameter. Note, however, that because of the independence between force/moment direction and intrinsic force generation described above, this simplification would not be expected to alter the estimated origin and insertion of each muscle. Similarly, although equations other than the one used here could be used to characterize the force length relationship of muscle, it is unlikely that their use would affect the estimated anatomical parameters. Further experiments in which muscle actions are characterized dynamically or across a wider range of conditions might make it possible to estimate these force generating properties of muscles more completely. Note also that although the experiments described here involved extensive dis96  sections to isolate individual muscles, other less invasive methods of estimating muscle actions have been developed [101]. Such less invasive methods might be used in combination with the computational analyses described in this study to construct musculoskeletal models.  3.4.3  Conclusion  The present study demonstrates the potential utility of in situ measurements of muscle actions to validate and develop musculoskeletal muscle models. However, it is important to point out several limitations to the present study. First, the measurements made here are obviously subject to experimental error. This error arises for a number of reasons, including imperfect mechanical isolation of the hip, muscle fatigue both within and between trials, and idealized kinematic models of joints. For instance, the z-moments observed for GRp could not be explained with the models considered here. In fact, the negative z-moments observed from stimulation should be impossible for any muscle if the hip were completely mechanically fixed. Obtaining ideal mechanical fixation is clearly difficult because of the considerable internal flexibility in the sacrum and pelvis: despite placing bone screws on the contralateral pelvis and on the rostral ipsilateral pelvis, there could still be some movement of the hip during muscle stimulation because of internal flexion of the pelvis. We were unable to restrain the skeleton further without interrupting the insertion of BFp on the pelvis. Although this imperfect fixation can apparently alter the measured actions of muscles, it is important to note that these effects did not substantially affect the models identified in these experiments: if we excluded the z moment from the optimizations, the anatomical parameters estimated were affected only minimally. Using more accurate joint models and tracking the exact 97  skeletal configuration at each limb configuration would also potentially increase the accuracy of models developed using the methods described here. It will be interesting in future work to evaluate the sensitivity of the models identified using the “top down” approach of this study to these sources of error. Whether or not the consequences of these measurement errors are greater in this top down approach than they are in more traditional “bottom up” approaches based on anatomical measurements [159], it seems clear that combining the information from both approaches can be a useful way to create and refine complex musculoskeletal models.  98  Chapter 4  Modeling Sensorimotor Behavior 4.1  Background  In the previous two chapters, we have focused on the mechanical aspects of the movement: how the movement is generated by mechanical interaction between myofilaments and how individual muscle contractions are related to body movement through the musculoskeletal structure. In this chapter, we change our focus to the highest level, the behavioral aspect of the movement: how movement is planned and controlled in response to sensory stimuli. At this level, movement is not just a mechanical phenomenon. It should be viewed as a product of behaviorallevel dynamics, the coordination between sensation and movement. Specifically, we focus on the eye-hand component of sensorimotor coordination, called visuomotor coordination, in object interception. Visuomotor coordination of eye and hand in object interception is one of the most important parts of human movement and has been extensively studied in the movement sciences. However, even though some one dimensional models have been proposed (e.g., [34]), 99  much of this work focuses on finding a descriptive model. Unlike previous studies, we propose a fully generative, computational model that can simulate realistic interception behavior in three dimensional space. The model is grounded in experimental observations of human behavior, both in our own experiment and reported in the literature. A key feature of the model is that it is based on simultaneous measurements of eye movements along with traditional motion capture of hand and body movements. We combine a model of active vision with a model of movement generation using short duration discrete submovements, based on experimental observations. More importantly, we propose that submovements of the eyes, head, hand, and body are tightly synchronized, which provides an exceptionally convenient way to produce natural looking coordination. Since the proposed model is not based on raw motion capture data but on the underlying principles, it generalizes very well to novel scenarios, such as catching with poor visibility and sudden changes in trajectory. The proposed model is based on our preliminary data. Therefore, our goal in this chapter is to develop a computational model that can be practically used for computer animation. We do not focus on rigorously verifying the model via a statistical investigation. Nevertheless, such a model is still valuable for applications in computer graphics because it captures significant features of real human behavior.  4.2  Observations  To understand how normal, untrained people actually perform interception tasks, it is helpful to look at some data. We collected pilot data of catching behavior of an untrained subject, simultaneously measuring the motion of the ball, hand, head, and body using a motion capture system, and eye movements using a head mounted 100  eye tracker. See Fig. 4.1; details of the measurement system are given in Sec. 4.6. First, the ball trajectory may be simple but is not known to the human catcher who has to quickly estimate it, primarily using vision. This immediately suggests that eye and head movements that enable clear vision of the moving ball are an important part of real catching behavior. This can be seen even by informal observation and in Fig. 4.1(a). This is why we pay special attention to eye movements, and not just to body movements as is the usual practice in full body animation (i.e., other than in face animation). Two types of eye movements are observed (see Fig. 4.2). (1) smooth pursuit eye movements which attempt to continuously track the ball, especially prominent after the ball has passed the apex of its trajectory. (2) catch-up saccades, which are very fast movements, especially at the beginning when the quality of pursuit is not good. Second, it is not uncommon for the flight duration to be very short, about a second, and human vision is relatively slow, so there is very little time to obtain an accurate estimate and plan the catch before moving the hand. Instead the hand begins to move very early, when only a crude estimate is known, and appears to have an initial “open loop” phase, followed by a “closed loop” or “homing” phase. This was observed a century ago by Woodworth [197], who proposed this “two component” model. Thus hand movements are not straight, preplanned trajectories, but are curved, with a stereotypical initial phase, followed by individual corrections to the target. See Fig. 4.1(b). Third, much experimental evidence suggests that human movements are not planned in their entirety but generated by blending together discrete, short-duration submovements, which is a long-standing hypothesis in movement science [184]. There has been considerable experimental support for this hypothesis, starting from 101  the seminal work of Soechting and Terzuolo [163] and Vallbo and Wessberg [184]. Although there are criticisms [165] that this may not apply to rhythmic movements, the concept of submovements continues to have wide support and is used in recent analyses on hand movement [132] and clinical treatment after stroke [35]. The example from pilot data shown in Figs. 4.1(b) and 4.5 show that real hand trajectories appear to have discontinuous higher derivatives, consistent with the hypothesis. This hypothesis is particularly appealing for computer animation, as we can leverage standard tools for constructing smooth curves from basis functions. Finally, somewhat surprisingly, we observed close synchronization between the submovements detected in the hand movements and catch-up saccades of the eyes. This phenomenon can be seen in Fig. 4.3 and consistently in other trials. Saccades and the peak hand velocity appear to be very well synchronized. Furthermore, if we decompose the hand velocity profile into submovements, we find that the start times of the submovements are also well synchronized to saccades. We can also qualitatively observe from decomposition results in Fig. 4.5 that the three dimensional hand trajectory is well represented by piecewise linear segments that correspond to the submovement directions. This analysis suggests an appealing hypothesis about movement generation, at least for visually guided interceptive movements: the eye and hand share the same “motor program” that is triggered by sensory events. Analogous hypotheses have been proposed by others, in slightly different contexts, e.g., [61, 87, 164]. This primacy of gaze in locomotion is neatly summarized by Berthoz [11] in the phrase: “ ‘Go where I’m looking’, not ‘Look where I’m going.’ ” It also underlies the coach’s admonition to “keep your eyes on the ball” to improve performance. This hypothesis deserves more detailed study in the movement sciences, presumably us102  (b)  (a)  Figure 4.1: Captured catching behavior: (a) Temporal change of the upper body configuration. (b) Captured hand trajectories for eleven trials (red) and corresponding ball trajectories (gray). ing a minimalistic highly controlled experiment design, as in the work cited above, but that is not our goal here. Rather, our goal is to show that we can exploit it to efficiently generate realistic animations of goal directed behavior of a 3D character for computer graphics.  4.3  Related Work  Due to its importance, object interception and its related topics have been investigated. We briefly review the most relevant work here. Computer Graphics. There have been very few papers written specifically on the topic of interception 103  vertical 0.1 rad 0.1 rad  horizontal  Figure 4.2: Eye movements for three different ball trajectories, as seen in a reference frame fixed at the initial position of the eye. The ball (gray) goes upward initially and then downward after passing the apex (green circle); note that due to perspective projection, the true apex differs from that of the projected trajectory. Gaze position (red) are captured every 10ms (dots). movements. Gillies and Dodgson [50] proposed a method to solve the classical “baseball outfielder problem” of how to run to catch a ball, but no implementation results were reported. More recently, several character animation papers included object interception as an example [3, 28, 43], but did not consider gaze. Papers on object manipulation have either not included gaze (e.g., [116, 142]) or included gaze as a post-process for realism (e.g, [177, 198]). Gaze has been an important part of virtual characters [48, 59, 83, 108, 137, 138, 198], and crowd simulation [57, 161]. Perhaps the closest to our view of coupling perception and movement are the seminal papers by Terzopoulos and co-workers (e.g., [107, 173, 178]). Movement Neuroscience. Visuomotor coordination and object interception have been extensively studied in  104  Hand Tangential Velocity (m/s) 60  70  80  90  16  1.2  12  0.8  8  0.4  4  0.0  0.1  0.2  0.3  0.4  0.5 Time (s)  0.6  0.7  0.8  Saccade Velocity (rad/s)  0  1.6  0  Figure 4.3: Tangential velocity of captured hand (green) and eye movements (red): Red arrows are observed saccades of which the estimated trigger times are drawn as dotted lines. Gray curves shows the decomposed submovements. the movement sciences. We refer the reader to a review by Zago et al. [202]. Submovements. A long standing hypothesis is that human (and animal) movement is not generated continuously, but produced by combining discrete building blocks, called submovements. There has been considerable experimental support for this hypothesis, starting from the seminal work of Soechting and Terzuolo [163] and Vallbo and Wessberg [184]. Although there are criticisms [165] that this may not apply to rhythmic movements, the concept of submovements continues to have wide support and is used in recent analyses on hand movement [132] and clinical treatment after stroke [35]. Robotics. Several impressive visuomotor control models for catching have been proposed in robotics [8, 13, 72, 99, 148, 149, 160]. Since their main goal is to generate robust high performance catching, the proposed models are optimized for the robot system itself, rather than for producing human-like movement. 105  4.4  Visual Estimation of Target Movement  Human vision is remarkably complex and the subject of intense ongoing study. For computer animation we need a relatively simple model that captures the relevant features of human vision, but by no means all the complexities. It is commonly hypothesized that the brain is able to use the internal models of object dynamics [196] to estimate and predict the object’s state over time. We use a Kalman filter for predicting the state of ball, using plausible models of noise introduced by visual sensing.  Y Z  u  X  q E  d  p  O Figure 4.4: Eye coordinate frame and the corresponding uncertainty of an object: Uncertainty is represented as a multivariate normal distribution. The error covariance ellipsoid around the true object location p is calculated by its deviation θ from the gaze direction (X-axis) and the axis of rotation ω. Fig. 4.4 defines a spatial coordinate frame attached to the eye, with its origin at the center of the globe, X-axis aligned with the visual axis, and Z-axis vertical in a reference position looking straight ahead. For ball catching, we define the perceived position and velocity of the ball as a six dimensional state vector in the p  eye coordinate frame, x = ( p˙ ).  106  Table 4.1: Parameters for error standard deviations of foveal vision: units for last three rows are identical to their original unit.  Retinal Position Retinal Velocity Depth Depth Velocity  4.4.1  Parameter σqz , σqy σq˙z , σq˙y σx σx˙  Value 2.9 × 10−4 rad 0.05|q| ˙ 0.03d ˙ 0.5|d|d  Vision  We approximate the probability of the state of a moving point seen by the eye as a multivariate normal distribution in the eye coordinate frame. As shown in Fig. 4.4, the probability of the perceived position p and velocity p˙ of the target in the eye-fixed frame is:      p p¯  2 T   = N   , AΣ A  p˙ p˙¯ Σ=  1 diag([σx , σy , σz , σx˙ , σy˙ , σz˙]), α(θ )   [ω]θ O3×3  e A= , O3×3 e[ω]θ  (4.1)  (4.2)  (4.3)  where p¯ and p˙¯ are the true position and the velocity of the ball in 3 dimensional space, Σ is the error covariance matrix in eye frame and A the corresponding transformation matrix from the target to the eye center. The block diagonal elements of A represent the rotation needed for the eye to foveate the target and [ω] is a skew symmetric matrix of the axis of this rotation. The noise parameters in Eqn. 4.2 are chosen as follows and are summarized in Table 4.1. The location of a point on the retina is represented using spheri-  107  cal coordinates, sequential rotation around Z and Y axes of the eye frame, called retinal coordinates: q = (qz , qy )T . We set the standard deviations σqz and σqy to correspond to the spatial resolution of standard 20/20 vision. Using the small angle approximation, we can convert the error from retinal to eye coordinates: (σz , σy )T = d(σqz , σqy )T , where d is the distance from eye to the object. The ability to discriminate the difference in retinal velocity is roughly proportional to the actual retinal velocity. The Weber fraction of velocity detection is known to be around 5% [122]. Therefore, if the probability of perceived velocity is also assumed to be normally distributed over its true velocity, the standard deviation of the error, (σq˙z , σq˙y )T , can be chosen to be proportional to the actual velocity. Similar to the position error, the velocity error represented in eye coordinate will be (σz˙, σy˙ )T = d(σq˙z , σq˙y )T . Our ability to measure depth or change in depth is known to be much worse than detecting the retinal location and is perceived using binocular vision and other cues. Simulating stereo vision is beyond the scope of this study; we set the depth variance, σx and σx˙ , to be much higher and also proportional to the depth to approximate the vergence angle. Finally, the factor α accounts for the drop in acuity with distance to fovea, the central region of highest acuity.  4.4.2  Internal Model and Bayesian State Estimation  Using the vision model described above, a noisy observation of the target states is obtained. For simulation, we discretize the state update with a time step ∆ t = 20ms.  108  The internal model of the ball’s dynamics is represented as  xk+1 = Fxk + b,      I3×3 F= O3×3    I3×3 ∆t  O3×1  , b =  , I3×3 g∆t  (4.4)  (4.5)  where g is gravity vector, which we assume to be constant. In other words, we assume that the brain has prior knowledge about gravity (see, e.g., [120]). Using the vision model described above, we can then model the observer. If Rk is the orientation and rk the position of the eye frame at time k, the observed position and velocity of the ball zk is:  zk = Hk xk + hk + vk ,  (4.6)  where vk is the observation noise defined in Eqn. 4.3 while Hk and hk are the transformation from the spatial to the eye frame:   T    O3×3   Rk Hk =  , T O3×3 Rk    T −Rk rk  hk =  . O3×1  (4.7)  We use a standard Kalman filter as a model of Bayesian inference with internalized dynamics in the brain. For completeness, the predict-update filtering algorithm is given below.  109  • Predict:  xk|k−1 = Fxk−1|k−1 + b  (4.8)  Pk|k−1 = FPk−1|k−1 FT  (4.9)  • Update:  yk = zk − Hk xk|k−1 − hk  (4.10)  Sk = Hk Pk|k−1 Hk T + Ak Σk 2 Ak T  (4.11)  Kk = Pk|k−1 Hk T Sk −1  (4.12)  xk|k = xk|k−1 + Kk yk  (4.13)  Pk|k = (I − Kk Hk )Pk|k−1 ,  (4.14)  where Ak and Σk are previously defined in Eqn. 4.3. Since we assume that the subject initially has no prior knowledge about the states, the initial value of the prior estimate covariance, P0|0 , is made sufficiently large and x0|0 is chosen to be the zero vector. The performance of the Kalman filter in object tracking is known to deteriorate seriously if the tracked object changes its trajectory suddenly, especially when the error covariance is small, such as the case when the ball we are tracking bounces unexpectedly. Therefore, in order to implement human vision’s ability to track an object with a sudden trajectory change, we should endow an adaptive behavior to the filter. This is typically achieved by resetting the prior error covariance matrix Pk|k−1 to a sufficiently high value when an abnormal tracking error is recognized.  110  We choose the reset decision variable as a Mahalanobis measure of the innovation: yk T Sk −1 yk , which indicates the abnormality of the error between prediction and observation.  4.5  Movement Generation  Based on the visually estimated target movement, we can then generate the gaze movement to the target (Sec. 4.5.1), followed by synchronized movements of the head (Sec. 4.5.2), the hand (Sec. 4.5.3), and the body (Sec. 4.5.4). To make this complex whole body movement tractable, the movements are generated kinematically, rather than using a dynamic model of human biomechanics. This is an approximation we make primarily for efficiency, but one that is well supported by human movement research. Even though dynamic properties of the body are important for the ultimate control of the movement, the brain appears to represent external movement goals kinematically at a high level. This can be seen by the fact that one’s handwriting looks the same regardless of posture or writing tool used; this phenomenon was termed “motor equivalence” by Donald Hebb [62]. Since then, several kinematic invariants and laws have been proposed for movement planning, though the topic remains controversial, as are most topics in neuroscience. The dynamics must influence the movement at some level, at least in the lower level controllers. We take this evidence pragmatically to mean that a kinematic approximation is plausible, and helps to improve the efficiency of the animation system. As discussed in Sec. 4.1 we assume that movements are produced by blending short duration submovements together. Each submovement is smooth, with a bell-  111  shaped velocity1 profile. Different smooth shapes have been proposed, including those that minimize jerk (time derivative of acceleration) [41] and the more flexible delta lognormal shape [141]. We chose the minimum jerk profile for simplicity, as it has fewer parameters. The tangential velocity v(t 0 ,t f ,t) with unit displacement is defined in a time interval t 0 < t < t f as:  v(t 0 ,t f ,t) =  (t f  30 (t − t 0 )2 (t − t f )2 . − t 0 )5  (4.15)  ˙ is represented as a superposition of submovement The velocity of a movement u(t) velocity profiles u˙ i : n  n  i=1  i=1  ˙ = ∑ u˙ i (t) = ∑ bi v(ti0 ,tif ,t), u(t)  (4.16)  where bi is the basis vector representing the direction and magnitude of the submovement. Very fast movements, such as saccades, can consist of single submovement, and slower movements of head and hand can consist of multiple submovements. Care should be taken in defining superimposed submovements: when the destination of the new submovement is given, bi is not the difference between the destination and current position, but rather, it should be the difference between the new destination and the previous submovement destination. Fig. 4.5 shows the velocity profile and decomposed submovement of some measured hand trajectories. Decompositions are done by nonlinear least squares optimization with pre-chosen number of submovements: finding a set of submovement parameters, bi ,t 0 ,t f , that minimizes the error between captured and com1 More  accurately, speed, but we follow the usage in the human movement sciences.  112  posed trajectories. Note that this submovement decomposition may not be unique and several different combination can exist depending on the initial condition and selection of the basis [154]. Using this framework, animation of body movement, including head movement and hand movement becomes very clear: for each update of visual information, corresponding submovements of each body part is determined and superimposed on the current motion.  4.5.1  Gaze Movements  Given the estimated ball state, we first determine gaze, that is, where to look. We define (binocular) gaze as a point in space where the eyes are looking. This determines the orientations of the eyes in space. As described earlier, a distinguishing feature of our approach to animation is that the gaze is determined prior to determining body movement. In humans and other animals, gaze is not significantly affected by body movement unless the eyes reach a biomechanical limit. Indeed, several fast low-level reflexes, including the vestibulo-ocular and vestibulo-collic reflexes, stabilize gaze in space and cancel disturbances due to body movements. To foveate an object of interest, two types of gaze movements are used, called saccade and pursuit. The saccade is a very fast, impulsive eye movement (up to 800◦ /s in humans) [111], but the high speed comes at a cost: vision is very poor during a saccade. Humans and other primates can use smooth pursuit eye movements to continuously track moving objects, but only at relatively low speeds (an order of magnitude slower than saccades). The properties of saccades and pursuit are well studied [152, 200], in particular [133] provide a recent perspective. Our goal here is not to precisely represent 113  5 4 3  5  4  2  3  1  1  2  0  2  4  6  0  8  2  4  6  8  10  5 4  5  4 3  3  2  2  1  Tangential Velocity  1  0  2  4  6  8  10  12  0  2  4  6  8  10  Time (100 msec)  Figure 4.5: Decomposed submovement of hand trajectory: [first and third row] Real (green) and synthesized (orange) submovements shown in three dimensional space. Dotted lines show the direction of each submovement. [second and bottom row] Decomposed submovements shown in tangential velocity  114  Perceived target position  Gaze  200ms +  -  Ax+b  Saccade velocity  Position error  +  Sampler  Saccade generator  Pursuit generator  Limiter  Integrator  Pursuit velocity  Prediction  Internalized dynamics  +  Figure 4.6: Overall control structure of gaze the decision to saccade or pursue, but rather to have a simple model that captures the statistical regularities observed in the pilot data. The typical time required for two consecutive saccades is known to be around 200ms [151, 200], corresponding to the 200ms recharge duration of the saliency map in the model proposed by Itti et al. [84]. Tracking fast moving objects is challenging, and is likely to require a rapid sequence of catch-up saccades. As shown in Fig. 4.3, we also observe this in our pilot data, with catch-up saccades triggered at approximately 200ms interval. Therefore, we assume 200ms as the visuomotor update time that regulates the general movement including head and hand movement. The overall block diagram of the eye movement velocity generator, including saccade and pursuit, is described in Fig. 4.6. The saccade amplitude is determined by current error between gaze and target angles. Given its amplitude, the likely duration of a saccade can be well estimated from the “main sequence” relationship [111]. The velocity of saccade is chosen as in Eqn. 4.15; real saccades have a more asymmetric velocity profile, but the difference is imperceptible for computer animation.  115  Since pursuit is a smooth movement, we set its velocity to be updated at every ¯ time step ∆t with v = (p−p)/∆t. Here, p is the current gaze and p¯ is the position of the target at the next step as predicted by the internal model. The maximum pursuit velocity is limited to 100◦ /s as per human oculomotor limits [124], restricting it to follow only low speed targets. When an object is detected, the estimated positions in the early stages will be very noisy. Therefore, if the saccade is simply triggered to the estimated position, the initial eye movement will be very inaccurate and unrealistic. The brain avoids this problem by suppressing the initial saccade until the likelihood of position of the object exceeds some significance criterion [22]. We relate this decision criterion to the innovation observed by the perceptual system, using the error covariance of innovation that is defined in Eqn. 4.11: the initial saccade is triggered when the estimated error of the object position in the retinal coordinate, which is the Y-Z coordinate of the eye frame, is reduced below a certain threshold.  4.5.2  Head Movements  We assume that the main goal of head movement in target tracking is to assist the eye in tracking the target more easily. Therefore, the movement for the head is defined similarly to the gaze: we define head movement by setting the “head gaze”, a spatial point that the head is facing. Since the rotation of the head follows Donders’ law [29], which states that three dimensional head rotation can sufficiently be described by two parameters, the longitudinal rotation followed by latitudinal rotation, the rotation matrix of the head, Rh , is determined by the head gaze. If the head gaze is represented with spherical coordinates, [qz , qy ], with respect to the head frame in the home configuration, Rh0 , the corresponding head orientation is 116  R = Rh0 Rz (qz )Ry (qy ). It is known that the same brain areas are involved in eye and head movement in both saccades [45] and pursuit [104]. Therefore, we assume that the head also generates saccades and pursuit in the same way as the eye, whereas the head saccade is a slower, longer lasting bell shaped velocity profile compared to the eye. Given a gaze shift, the corresponding saccade amplitude of the head is linearly related with a 20◦ dead zone [111]. There is a more sophisticated model proposed by Freedman [44] for the kinematic relationship between the head and eye movement, but this model is not directly applicable to 3-dimensional, large angle movements since it is focused on horizontal movements with a limited range. Once the change of head gaze is determined, the peak velocity is also determined by its well known linear relationship to the amplitude [111], with slope varying between 4s−1 to 8s−1 . If we apply this relationship to the submovement shape function in Eqn. 4.15, we get a constant submovment duration around 400ms. Taken together, we can define the submovement parameters in Eqn. 4.16 for saccade velocity of the head gaze, g˙ :  ξ  = max(0, θ − 0.349)  bi = (e[ωi ]ξ − I)(gi−1 − p)  (4.17) (4.18)  where [ωi , θi ] is the axis-angle representation of the displacement from the last head gaze destination gi−1 to the current eye gaze destination with respect to the current head position p, and ξ is the required angular displacement of the head saccade. Here θ is the non-negative magnitude of the rotation, with the sign of the rotation absorbed in ω. Note that, this definition will not generate a minimum jerk 117  angular velocity profile, but it is sufficiently close when the gaze is not very near. The pursuit velocity of the head, which is identical to that of the eye, is then added.  4.5.3  Hand Movements for Interception  With predetermined gaze behavior, the strategy of manual interception can be modeled as a simple algorithm. When the gaze is successfully tracking the ball, we can simply move the hand towards the gaze. In other words, if the hand is “latched” to the eye frame while the gaze is latched to the target, the hand will be driven sufficiently close to the ball when it reaches the body. This strategy can be observed when we try to catch a fly. Before we trigger the final snatch, we move our hand in the same way that we move our eyes and head. When the hand is synchronized to the gaze, the remaining decision will be to determine the distance of the catch. The pilot data shown in Fig. 4.7 strongly suggests that humans have a preferred interception distance from the head frame. This distance may be affected by considerations such as the manipulability of the arm or the effective compliance produced by muscle properties, but for simplicity we choose a fixed distance as a decision variable for catching. If we draw a sphere whose radius is the preferred distance around head, the catch position is determined as an early intersection point between the sphere and the predicted ball trajectory. Note that this interception planning is only valid when we have a sufficiently accurate prediction of the object trajectory. For this decision, we use the determinant of the current posterior error covariance of the state, defined in Eqn. 4.14, as a decision variable. Note that a decision made too early can cause a serious mistake, whereas one that is too late gives too little time to move. Therefore, we empirically searched 118  Distance from head to the interception point (mm)  1000  800  600  400  200  0  5  10  15  20  25  30  35  40  45  50  55  Trial Number  Figure 4.7: Distance from head to the interception point summarized from 59 ball catching trials for one subject for the best decision timing and corresponding value of the decision variable. Interestingly, from the simulation result, we see that in most cases a good decision is made around the apex of the ball trajectory as illustrated in Fig. 4.8. This agrees with the transition from catch-up saccades to pursuit that is previously shown in Fig. 4.2, and also corresponds to the large change of the hand trajectory observed in the pilot data. Fig. 4.9 shows the captured hand trajectory for different catching trials and we can see that a clear branch of the trajectory toward the target occurs around the apex of the ball. The reason why the decision is made around the apex can be intuitively understood: since the velocity of the ball reaches its minimum at the apex and the gaze usually catches up with the ball by saccades and start to pursue, the quality of the position and velocity estimation is considerably improved around the apex. This corresponds to the fact skilled jugglers only look at the juggled balls at their apex [9]. Based on this hypothesis, we divide the interception strategy into two parts: 119  60  Decision variable = log( det[ Covariance ] )  40  20  0  Decision criterion −20  −40  0  20  40  60  80  100  120  Time (10 msec)  Figure 4.8: Decision variables and ball apex: Lines show the decision variable, the determinant of the posterior error covariance matrix, for various catching simulations. For each trial, the time when the ball reaches its apex is marked as a green dot. Gray dotted line shows the criterion used for the prediction. a reactive phase and a proactive phase, similar to Woodworth’s two component model. The strategies for composing submovements for each phase are as follows: Reactive Phase. A submovement to drive the hand to the gaze direction with preferred catching distance is generated. The destination of the submovement is not exactly set to the gaze center to avoid the occlusion of the target due to the hand. For a right-handed subject, the destination is set to be the 45 degrees lower right point with respect to the gaze that corresponds to the captured configuration in Fig. 4.9. We name this position as the ready pose, pready , that is fixed with respect to the eye frame. If Re and pe are the orientation and position of the eye frame, and p˜ m is the destination of the previous submovement, direction vector b in Eqn. 4.16 can be  120  Figure 4.9: Transition of the hand trajectory around the apex (from pilot data): The subject’s pose corresponds to the ball apex. Red curves show the entire hand trajectories. It can clearly be seen that significant modification of the trajectory occurs around the apex. determined as follows: b = Re pready + pe − p˜ m  (4.19)  We limit the maximum velocity to 2m/s as observed in the pilot data. We also assume from the pilot data that the duration of the submovement is approximately 400ms so that it reaches the maximum velocity when the next submovement is triggered. This is consistent with the observation of Novak et al. [132]. Proactive Phase. When the subject has a sufficiently accurate prediction of the ball trajectory, the observed ball position and velocity can be extrapolated into the future to determine the interception point. Here we limit the prediction to a specific time window, 400ms, that reflects the prediction limit of the internal model. 121  When the ball penetrates the sphere of the preferred distance, the interception point and remaining time is determined and a corresponding submovement is generated. Otherwise, the hand goes in the direction of the closest point on the predicted trajectory with maximum 2m/s velocity. The magnitude of the submovement is also adjusted so that the final hand position does not exceed the biomechanical limit of reach. Using this strategy, we can generate the distinctive curved path of the hand trajectory observed from the pilot data. Fig. 4.10 shows the resultant hand trajectory and submovements. The orientation and grasping angle of the hand can also be determined in this phase by applying a relationship observed in the pilot data. From the data, we observe that the angle of attack between the ball and the palm remains at 90 degrees and grasping starts with the last submovement, when the decision to catch is made.  4.5.4  Body Movements  As previously described, we focus on generating coupled motion between gaze, head, and hand. Body movements are not the focus of this work as they have been extensively studied in computer animation. For a given trajectory of the hand and head, we solve for the configuration of the rest of the body using inverse kinematics (IK). The character model used consists of a torso that is connected to the ground by a 3 DOF rotational joint and a 7 DOF (right) arm model. The position of the head is fixed on the torso. The input provided to the IK algorithm consists of the position and the normal vector of the palm that is adjusted to satisfy the preferred angle of attack between the ball and the hand. Although general treatment for IK is not our main focus, its solution is not trivial due to kinematic redundancy, with biomechanical constraints and learned 122  4 3 5  6 5  2  4 6 3  2  1 1  Figure 4.10: Submovement composition of the hand trajectory for a given ball trajectory (black curve): Left: Circles show the perceived position of the ball over time. The simulated hand trajectory (orange curve) is compared to the captured pilot data (green curve). As the catch phase changes from reactive (gray circles) to proactive (red circles), corresponding submovements are generated. Right: dotted lines (gray for reactive and red for proactive phase) are directions and amplitudes of submovements triggered at each visuomotor update (numbers match to the estimated ball position on the trajectory). styles contributing to realistic solutions. We first implemented the IK solution proposed by [58], but it did not generalize well when the movement approaches or crosses the boundary of the training set where the learned movements need to be extrapolated. More recent latent variable models [179, 186] may fare better, but we instead used a simpler approach of combining example configurations with weighted pseudo-inverse of the Jacobian to obtain our IK solution (for all 3+7  123  DOF). Our approach first solves for the torso and arm joint angles according to the specified hand position, we then solve for the relative hand orientation according to the ball velocity. To solve for the torso and arm configuration, we first search in our training set (derived from pilot motion capture data) for configurations with hand position most similar to the target position. The average of these configurations is then blended with the previous timestep configuration and a weighted pseudo-inverse of the Jacobian is used to correct the resulting configuration so that its hand position matches the target. During the reactive phase, the rotation of the hand is interpolated from its home orientation to a predetermined ready pose orientation with a set velocity. Subsequently, in the proactive phase, we find the rotation that gives the nearest interpolation from the previous palm normal vector to the estimated ball velocity vector and apply it to the hand. Finally, we use a canned animation to generate the finger grasping motion during the final moments of interception.  4.6  Results  Using the proposed methods, we developed a fully generative 3D character animation model of the ball catching, shown in Fig. 4.11. Validation. We compared the results of our simulation to measurements of human catching, to see if there is a qualitative agreement. Note that the measurements were not used for estimating the parameters of the model, so qualitative agreement is the best one can expect. We recorded the movement of the upper body using an 8-camera Vicon MX motion capture system (Vicon, Los Angeles), and eye movements using a head mounted C-ETD eyetracker (Chronos Vision, 124  Figure 4.11: Simulation of three dimensional ball catching. The full eye and body movements of the character are automatically generated by given ball trajectory as a visual stimuli. Berlin) to capture the head-unrestrained eye movements. Both systems recorded at 100 Hz, which was sufficient for our needs. Two subjects were seated on a comfortable stool and instructed to catch a 5cm diameter ball (a ping-pong ball filled with water) thrown in the subject’s direction. The gaze behavior for different ball trajectories was previously shown in Fig. 4.2. The simulated data generated using our system is shown in Fig. 4.12 and looks qualitatively similar. Comparison of measured and simulated hand trajectories in Fig. 4.13 shows that the proposed algorithm is sufficient to generate qualitative details of real hand movement as well.  125  Figure 4.12: Simulated gaze behavior for a given ball trajectory: Colors and data rates are the same as Fig. 4.2.  Generalization. We simulated ball catching for many different trajectories. The character is able to catch the ball successfully in 98 percent of the trials. To test the model capability to generalize beyond simple trials, we introduce several variations in the catching scenario. In one scenario, we assume that the ball unexpectedly bounces off a transparent wall. We found that the Kalman filter with adaptive resetting successfully tracks the ball after the bounce. The character often fails to catch if the time to recover the prediction is too short, just as humans would. In another scenario, we simulated poor vision by increasing the standard deviation of the error in Table 1 by a factor of 10. As expected, the resultant motion of character showed a slow start-up followed by an abrupt catching with an 81.2 percent success ratio. Performance. By design, the model is very efficient. We implemented the entire algorithm in MATLAB (version 7.9.0, R2009b, MathWorks Inc.) without using externally compiled subroutines. The software ran on a 3.40 GHz Intel Core i7 computer. The computation for generating the entire catching movement, in-  126  Figure 4.13: Comparison between real (green) and simulated (red) hand trajectories for the same ball catching tasks: Submovements vectors are plotted as a gray dotted lines. Note that in the bottom right trials, the simulated character failed to catch the ball since it is subjected to a very extreme condition. cluding IK, for 100 frames, with a 1 second ball trajectory, takes 5.93 seconds (16.9 fps).  4.7  Conclusions  In this chapter, we have proposed a framework of animating visually guided interception based on the neurophysiology of the visuomotor system. The framework includes novel features such as gaze based motor coordination and submovement composition, that have not been previously introduced in computer animation.  127  Since the model is grounded in human behavior, the simulated catching movements look very human-like. Interestingly, we found that the use of submovements produces very subtle behaviors such as discontinuities and hesitation, which are clearly different from continuous, robotic motions. Our framework currently has several limitations. These are mostly due to this study’s focus on constructing a complete and practical generative model for computer animation, which required simplifying many of the known complexities of the human sensorimotor system. First, our vision model is a highly simplified approximation of the human visual system. Our linear Kalman filter could be improved using more sophisticated (and complex) Bayesian filters. Although a visuomotor update interval of about 200ms is observed in both our pilot data and in previous studies, it is likely that there is some variability as a result of a probabilistic decision process. A more accurate model of this decision process may be able to account for the variability. In general, it may also be desirable to model context-dependent and individual variability in the coordination pattern. Finally, we currently rely on kinematic control, but it is clear that a dynamic model, based on realistic biomechanics and neural control (e.g., [110, 167]), may be able to better account for movements at the limits of human performance. Notwithstanding the limitations listed above, we believe our framework provides a new approach to computer animation, based on principles that may be used by the human brain to control movement, such as visual motion estimation, use of gaze to control body movements, and generating movements using overlapping submovements. Since we base our animations on these principles, rather than on a corpus of raw data, the method is very efficient and generalizes in plausible ways to novel scenarios. 128  Chapter 5  Conclusions In this thesis, we have presented a set of studies that model different layers of the neuromusculoskeletal system. We have especially focused on building computational models that can simulate real data in a sufficiently wide range of space and time. There are three distinctive levels of modeling, and conclusions for each part are presented in Chapter 2, 3, and 4. Here, we provide a brief overview and directions for future work.  5.1  Overview of Contributions  In muscle mechanics, the thesis has made two major contributions to the current research literature. First, we identified and named two different forms of the Hill type muscle model presented in previous studies which were not well distinguished. Based on the definition, we have shown, for the first time, that the isotonic shortening profile of skeletal muscle can be quantitatively simulated well using the rather obscure f-max scaling model, but not with the standard force-scaling model. The winding filament model presented in the latter part of the analysis of mus129  cle mechanics gives us a new insight into muscle mechanics and provides a possible explanation for many unsolved problems. Together with the explanatory value of the hypothesis, we have proposed a computational model of residual force enhancement that is based on a succinct protein level structure and is able to produce a realistic pattern of the enhancement that qualitatively matches previous observations. An important contribution of this thesis to musculoskeletal modeling is that we have presented a general method to build a physically realistic musculoskeletal model. Two points have been raised: 1) we have shown that the action of even some complex muscles can be simulated by the traditional line segment based model. The successful results from two geometrically different muscles, GRp and BFp, have been presented. 2) However, the geometric modeling parameters for the model should be chosen through careful consideration, as they sometimes deviate quite a bit from the intuitive choice based on anatomy. As a convincing example, we have shown that the model of the BFp muscle with a conventional insertion point can produce a significant amount of error. This finding has been explained later by closer anatomical inspection on the insertion area of BFp, by which we have revealed a special tendinous structure that proximally distributes muscle fibers. The last part of the thesis contributes to the development of a behavioral-level simulation model. First, we suggest for the first time a fully generative model of visuomotor coordination in three dimensional space. Compared to previous models that describe behavior [87] or simulate data in a limited dimension [133], this model can simulate observed data in three dimensional space. Furthermore, although it needs further investigation, we have suggested a novel hypothesis con130  necting submovements of the limbs and eye movements, which may explain how submovements are composed in coordination with vision. The second contribution of the work is in the area of computer graphics; our study, for the first time, introduces the concept of visuomotor coordination to human character animation. The ability to automatically generate realistic human-like interception behavior from a given visual stimulus can be widely used in various applications of character animation such as in computer games and virtual agents. Overall, this thesis has applied a computational approach to various modeling problems of the neuromusculoskeletal system, spanning from filaments to behavior. The results of each chapter have shown that the proposed computational models provide insights into the working principles of this complex system, guide future experiments, and are also effective in application areas such as computer animation of human movement.  5.2  Future Work  Computational modeling of the neuromusculoskeletal system is a relatively new area that is still at an initial stage of development. Therefore, the work presented here leaves plenty of room for improvement. For the muscle mechanics model, as discussed, a direct validation of the winding filament hypothesis is currently the most urgent step. A new technique should be developed that enables nanometer-scale observation of titin winding in an active single myofibril. On the modeling side, future work should focus on integrating the model of titin winding with the traditional cross-bridge theory and to eventually develop a fully computational model that accounts for the viscoelastic behavior of muscle. 131  Future work in musculoskeletal modeling should develop a comprehensive, physically valid model of the whole rat hind limb that includes all major muscles. Together with collecting more data on other muscles, analyzing parameter sensitivity using the Monte Carlo method [159] and introducing a free-form curve-based model of muscle geometry [167] are currently in progress [189]. For the model of visuomotor coordination, the proposed model originally focused on developing a generative model that can be used in computer animation; therefore, the scientific validation of the hypothesis was relatively less emphasized. We are currently planning to validate the hypothesis with more subjects in a more rigorously controlled environment. In a longer-term perspective, the ultimate goal of the study is to build a unified, multi-layered model of the neuromusculoskeletal system from muscle filaments to behavior. Although it may look overambitious to build such model, we see that our results from the different chapters of the thesis can possibly be connected together as follows: From muscle mechanics to musculoskeletal modeling. The proposed optimal design framework to simulate realistic static force and moment field can be extended to the stiffness and velocity dependent properties that are expected to be realistically simulated by the winding filament model. It would be an interesting future work to see how well these properties of individual muscle build up to those of the end effector by the point-to-point musculoskeletal model. From muscle mechanics to sensorimotor behavior. Our visuomotor coordination study assumes that the submovement of hand is planned as a kinematic quan132  tity. However, in reality, we also need to plan the dynamics along the trajectory for robust control of the movement, which is the main idea of impedance control [70]. From this viewpoint, the winding filament model provides a feasible mechanism by which the mechanical stiffness of muscle can be regulated by a history dependent winding mechanism. This connection can be elucidated and verified by realistic musculoskeletal model. It is commonly apprehended that this integrative approach of building a comprehensive, multi-layered model will result a system with huge number of parameters that are almost impossible to handle and to understand. However, it is important to note that the integration procedure does not always increase the complexity, but it also helps to reduce the complexity; the structural understanding of the inter-layer relationships provides a systematic way to evaluate the mechanical significance of each parameter to the entire system [180] and could eventually help us to find the minimal parametic representation of each layer to maintain the overall integrity of the system. For this reason, we expect that our integrative, computational modeling approach will contribute significantly to understanding the control of movement [129, 183].  133  Bibliography [1] B. Abbott and X. Aubert. The force exerted by active striated muscle during and after change of length. The Journal of physiology, 117(1): 77–86, 1952. → pages 36 [2] B. Abbott and D. Wilkie. The relation between velocity of shortening and the tension-length curve of skeletal muscle. The Journal of physiology, 120 (1-2):214, 1953. → pages 7, 19, 24 [3] Y. Abe and J. Popovi´c. Interactive animation of dynamic manipulation. In Proceedings of the 2006 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 195–204. Eurographics Association, 2006. → pages 104 [4] A. Arnold, F. Anderson, M. Pandy, S. Delp, et al. Muscular contributions to hip and knee extension during the single limb stance phase of normal gait: a framework for investigating the causes of crouch gait. Journal of biomechanics, 38(11):2181, 2005. → pages 4 [5] E. Azizi and T. Roberts. Muscle performance during frog jumping: influence of elasticity on muscle operating lengths. Proceedings of the Royal Society B: Biological Sciences, 277(1687):1523–1530, 2010. → pages 36 [6] E. Azizi, E. Brainerd, and T. Roberts. Variable gearing in pennate muscles. Proceedings of the National Academy of Sciences, 105(5):1745–1750, 2008. → pages 63 [7] T. Baechle and R. Earle. Essentials of strength training and conditioning. Human Kinetics Publishers, 2008. → pages 13 [8] B. Bauml, F. Schmidt, T. Wimbock, O. Birbach, A. Dietrich, M. Fuchs, W. Friedl, U. Frese, C. Borst, M. Grebenstein, et al. Catching flying balls 134  and preparing coffee: Humanoid rollin’justin performs dynamic and sensitive tasks. In Robotics and Automation (ICRA), 2011 IEEE International Conference on, pages 3443–3444. IEEE, 2011. → pages 105 [9] P. Beek and A. Lewbel. The science of juggling. Scientific American, 273 (5):92–97, 1995. → pages 119 [10] P. Bennett, T. Hodkin, and C. Hawkins. Evidence that the tandem ig domains near the end of the muscle thick filament form an inelastic part of the i-band titin. Journal of structural biology, 120(1):93–104, 1997. → pages 42, 48 [11] A. Berthoz. The brain’s sense of movement. Harvard Univ Pr, 2000. → pages 102 [12] B. Bigland-Ritchie and J. Woods. Integrated electromyogram and oxygen uptake during positive and negative work. The Journal of physiology, 260 (2):267–277, 1976. → pages 39 [13] O. Birbach, U. Frese, and B. Bauml. Realtime perception for catching a flying ball with a mobile humanoid. In Robotics and Automation (ICRA), 2011 IEEE International Conference on, pages 5955–5962. IEEE, 2011. → pages 105 [14] J. Bordas, A. Svensson, M. Rothery, J. Lowy, G. Diakun, and P. Boesecke. Extensibility and symmetry of actin filaments in contracting muscles. Biophysical journal, 77(6):3197–3207, 1999. → pages 45 [15] B. Brenner. The necessity of using two parameters to describe isotonic shortening velocity of muscle tissues: the effect of various interventions upon initial shortening velocity (v i) and curvature (b). Basic research in cardiology, 81(1):54–69, 1986. → pages 19 [16] I. Brown and G. Loeb. A reductionist approach to creating and using neuromusculoskeletal models. Biomechanics and neural control of posture and movement, pages 148–163, 2000. → pages 5, 14, 33, 65 [17] I. Brown, S. Scott, and G. Loeb. Mechanics of feline soleus: II Design and validation of a mathematical model. Journal of muscle research and cell motility, 17(2):221–233, 1996. → pages 14 [18] T. Buchanan, D. Lloyd, K. Manal, and T. Besier. Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements 135  from measurements of neural command. Journal of applied biomechanics, 20(4):367, 2004. → pages 6 [19] T. Burkholder and R. Lieber. Sarcomere length operating range of vertebrate muscles during movement. Journal of Experimental Biology, 204(9):1529–1536, 2001. → pages 32, 48 [20] T. Burkholder and T. Nichols. Three-dimensional model of the feline hindlimb. Journal of morphology, 261(1):118–129, 2004. → pages 65 [21] K. Campbell. Interactions between connected half-sarcomeres produce emergent mechanical behavior in a mathematical model of muscle. PLoS Computational Biology, 5(11):e1000560, 2009. → pages 63 [22] R. Carpenter and M. Williams. Neural computation of log likelihood in control of saccadic eye movements. Nature, 377(6544):59–62, 1995. → pages 116 [23] D. Carrasco and A. English. Mechanical actions of compartments of the cat hamstring muscle, biceps femoris. Peripheral and spinal mechanisms in the neural control of movement, 123:397, 1999. → pages 93 [24] E. Chadwick, D. Blana, A. van den Bogert, and R. Kirsch. A Real-Time, 3D Musculoskeletal Model for Dynamic Simulation of Arm Movements. IEEE transactions on bio-medical engineering, 56(4), 2008. → pages 65 [25] C. Chanaud and J. Macpherson. Functionally complex muscles of the cat hindlimb. III. Differential activation within biceps femoris during postural perturbations. Exp. Brain Res, 85:271–280, 1991. → pages 95 [26] C. Chanaud, C. Pratt, and G. Loeb. Functionally complex muscles of the cat hindlimb. II. Mechanical and architectural heterogenity within the biceps femoris. Experimental brain research. Experimentelle Hirnforschung. Experimentation cerebrale, 85(2):257, 1991. → pages 93 [27] J. Clark, J. Cham, S. Bailey, E. Froehlich, P. Nahata, R. Full, and M. Cutkosky. Biomimetic design and fabrication of a hexapedal running robot. In Robotics and Automation, 2001. Proceedings 2001 ICRA. IEEE International Conference on, volume 4, pages 3643–3649. IEEE, 2001. → pages 4 [28] S. Cooper, A. Hertzmann, and Z. Popovi´c. Active learning for real-time motion controllers. In ACM Transactions on Graphics (TOG), volume 26, page 5. ACM, 2007. → pages 104 136  [29] J. Crawford, J. Martinez-Trujillo, and E. Klier. Neural control of three-dimensional eye and head movements. Current opinion in neurobiology, 13(6):655–662, 2003. → pages 116 [30] M. Damsgaard, J. Rasmussen, S. Christensen, E. Surma, and M. de Zee. Analysis of musculoskeletal systems in the AnyBody Modeling System. Simulation Modelling Practice and Theory, 14(8):1100–1111, 2006. → pages 19 [31] S. Delp and J. Loan. A graphics-based software system to develop and analyze models of musculoskeletal structures. Computers in Biology and Medicine, 25(1):21–34, 1995. → pages 64 [32] S. Delp, J. Loan, M. Hoy, F. Zajac, E. Topp, and J. Rosen. An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures. Biomedical Engineering, IEEE Transactions on, 37(8): 757–767, 1990. → pages 4, 10 [33] S. Delp, F. Anderson, A. Arnold, P. Loan, A. Habib, C. John, E. Guendelman, and D. Thelen. OpenSim: Open-source software to create and analyze dynamic simulations of movement. IEEE transactions on biomedical engineering, 54(11):1940, 2007. → pages 19, 78 [34] J. Dessing, C. Peper, D. Bullock, and P. Beek. How position, velocity, and temporal information combine in the prospective control of catching: Data and model. Journal of cognitive neuroscience, 17(4):668–686, 2005. → pages 99 [35] L. Dipietro, H. Krebs, S. Fasoli, B. Volpe, and N. Hogan. Submovement changes characterize generalization of motor recovery after stroke. Cortex, 45(3):318–324, 2009. → pages 102, 105 [36] K. Edman. The velocity of unloaded shortening and its relation to sarcomere length and isometric force in vertebrate muscle fibres. The Journal of Physiology, 291(1):143, 1979. → pages 19, 32 [37] K. Edman and R. Josephson. Determinants of force rise time during isometric contraction of frog muscle fibres. The Journal of physiology, 580 (3):1007–1019, 2007. → pages 17 [38] K. Edman, G. Elzinga, and M. Noble. Residual force enhancement after stretch of contracting frog single muscle fibers. Journal of General Physiology, 80(5):769, 1982. → pages 36, 56, 58, 60 137  [39] C. Eng, L. Smallwood, M. Rainiero, M. Lahey, S. Ward, and R. Lieber. Scaling of muscle architecture and fiber types in the rat hindlimb. Journal of Experimental Biology, 211(14):2336, 2008. → pages 78 [40] M. Epstein and W. Herzog. Theoretical models of skeletal muscle: biological and mathematical considerations. Wiley, 1998. → pages 14, 19, 32, 37 [41] T. Flash and N. Hogan. The coordination of arm movements: an experimentally confirmed mathematical model. The journal of Neuroscience, 5(7):1688–1703, 1985. → pages 112 [42] M. Forcinito, M. Epstein, and W. Herzog. Can a rheological muscle model predict force depression/enhancement? Journal of biomechanics, 31(12): 1093–1099, 1998. → pages 37 [43] J. Francik and A. Szarowicz. Character animation with decoupled behaviour and smart objects. In 6th International Conference on Computer Games CGAIMS, Louisville, Kentucky, USA, 2005. → pages 104 [44] E. Freedman. Interactions between eye and head control signals can account for movement kinematics. Biological cybernetics, 84(6):453–462, 2001. → pages 117 [45] E. Freedman, T. Stanford, and D. Sparks. Combined eye-head gaze shifts produced by electrical stimulation of the superior colliculus in rhesus monkeys. Journal of neurophysiology, 76(2):927–952, 1996. → pages 117 [46] T. Funatsu, E. Kono, H. Higuchi, S. Kimura, S. Ishiwata, T. Yoshioka, K. Maruyama, and S. Tsukita. Elastic filaments in situ in cardiac muscle: deep-etch replica analysis in combination with selective removal of actin and myosin filaments. The Journal of cell biology, 120(3):711–724, 1993. → pages 48 [47] S. Galler, K. Hilber, et al. Tension/stiffness ratio of skinned rat skeletal muscle fibre types at various temperatures. Acta physiologica scandinavica, 162(2):119, 1998. → pages 39 [48] M. Garau, M. Slater, V. Vinayagamoorthy, A. Brogni, A. Steed, and M. Sasse. The impact of avatar realism and eye gaze control on perceived quality of communication in a shared immersive virtual environment. In Proceedings of the SIGCHI conference on Human factors in computing systems, pages 529–536. ACM, 2003. → pages 104 138  [49] H. Geyer and H. Herr. A muscle-reflex model that encodes principles of legged mechanics produces human walking dynamics and muscle activities. Neural Systems and Rehabilitation Engineering, IEEE Transactions on, 18(3):263–273, 2010. → pages 10 [50] M. Gillies and N. Dodgson. Ball catching: An example of psychologically-based behavioural animation. Eurographics UK, 1999. → pages 104 [51] G. Gillis and A. Biewener. Hindlimb muscle function in relation to speed and gait: in vivo patterns of strain and activation in a hip and knee extensor of the rat (Rattus norvegicus). Journal of Experimental Biology, 204(15): 2717, 2001. → pages 67 [52] S. Giszter, F. Mussa-Ivaldi, and E. Bizzi. Convergent force fields organized in the frog’s spinal cord. Journal of Neuroscience, 13(2):467, 1993. → pages 72, 91 [53] M. Goldstein, L. Michael, J. Schroeter, and R. Sass. Two structural states of z-bands in cardiac muscle. American Journal of Physiology-Heart and Circulatory Physiology, 256(2):H552–H559, 1989. → pages 55 [54] G. Golub and C. van Loan. Matrix Computations. Johns Hopkins University Press, 1989. → pages 75 [55] A. Gordon, A. Huxley, and F. Julian. The variation in isometric tension with sarcomere length in vertebrate muscle fibres. The Journal of physiology, 184(1):170, 1966. → pages 15, 16, 79 [56] E. Greene. The anatomy of the rat. Hafner New York, 1935. → pages 66, 91, 93 [57] H. Grillon and D. Thalmann. Simulating gaze attention behaviors for crowds. Computer Animation and Virtual Worlds, 20(2-3):111–119, 2009. → pages 104 [58] K. Grochow, S. Martin, A. Hertzmann, and Z. Popovi´c. Style-based inverse kinematics. In ACM Transactions on Graphics (TOG), volume 23, pages 522–531. ACM, 2004. → pages 123 [59] E. Gu and N. Badler. Visual attention and eye gaze during multiparty conversations with distractions. In Intelligent Virtual Agents, pages 193–204. Springer, 2006. → pages 104 139  [60] D. Gunzel and W. Rathmayer. Non-uniformity of sarcomere lengths can explain the catch-likeeffect of arthropod muscle. Journal of Muscle Research and Cell Motility, 15(5):535–546, 1994. → pages 36 [61] M. Hayhoe and D. Ballard. Eye movements in natural behavior. Trends in cognitive sciences, 9(4):188–194, 2005. → pages 102 [62] D. Hebb. The organization of behavior: A neuropsychological theory. Lawrence Erlbaum, 1949. → pages 111 [63] W. Herzog, T. Leonard, and J. Wu. Force depression following skeletal muscle shortening is long lasting. Journal of biomechanics, 31(12): 1163–1168, 1998. → pages 33 [64] H. Higuchi, Y. Nakauchi, K. Maruyama, and S. Fujime. Characterization of beta-connectin (titin 2) from striated muscle by dynamic light scattering. Biophysical journal, 65(5):1906–1915, 1993. → pages 48 [65] A. V. Hill. The heat of shortening and the dynamic constants of muscle. Proceedings of the Royal Society of London. Series B, Biological Sciences, 126(843):136–195, 1938. → pages 14, 16, 38 [66] A. V. Hill. The mechanical efficiency of frog’s muscle. Proceedings of the Royal Society of London. Series B, Biological Sciences, 127(849):434–451, 1939. → pages 19 [67] A. V. Hill. The series elastic component of muscle. Proceedings of the Royal Society of London. Series B, Biological Sciences, 137(887):273–280, 1950. → pages 17 [68] A. V. Hill. The effect of load on the heat of shortening of muscle. Proceedings of the Royal Society of London. Series B, Biological Sciences, 159(975):297–318, 1964. → pages 38, 62 [69] A. V. Hill. First and last experiments in muscle mechanics. Cambridge Univ Pr, 1970. → pages 19 [70] N. Hogan. Impedance control: An approach to manipulation. In American Control Conference, 1984, pages 304–313. IEEE, 1984. → pages 133 [71] K. Holzbaur, W. Murray, and S. Delp. A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control. Annals of biomedical engineering, 33(6):829–840, 2005. → pages 65  140  [72] B. Hove and J. Slotine. Experiments in robotic catching. In American Control Conference, 1991, pages 380–386. IEEE, 1991. → pages 105 [73] P. Huijing. Muscle as a collagen fiber reinforced composite: a review of force transmission in muscle and whole limb. Journal of biomechanics, 32 (4):329–345, 1999. → pages 63 [74] P. Huijing and G. Baan. Myofascial force transmission causes interaction between adjacent muscles and connective tissue: effects of blunt dissection and compartmental fasciotomy on length force characteristics of rat extensor digitorum longus muscle. Archives of physiology and biochemistry, 109(2):97–109, 2001. → pages 94 [75] A. Huxley. Muscle structure and theories of contraction. Progress in biophysics and biophysical chemistry, 7:255, 1957. → pages 34 [76] A. Huxley. A note suggesting that the cross-bridge attachment during muscle contraction may take place in two stages. Proceedings of the Royal Society of London. Series B, Biological Sciences, 183(1070):83–86, 1973. → pages 34, 38 [77] A. Huxley. Muscular contraction. The Journal of physiology, 243(1):1, 1974. → pages 38 [78] A. Huxley and R. Simmons. Proposed mechanism of force generation in striated muscle. Nature, 233:533–538, 1971. → pages 14, 26, 34, 38, 62 [79] H. Huxley. Fifty years of muscle and the sliding filament hypothesis. European journal of biochemistry, 271(8):1403–1415, 2004. → pages 34 [80] H. Huxley, J. Hanson, et al. Changes in the cross-striations of muscle during contraction and stretch and their structural interpretation. Nature, 173(4412):973, 1954. → pages 41 [81] L. Hyde, T. Burkholder, and T. Nichols. Reflex action of the hamstrings muscles at the feline ankle mediated by the crural fascia. Soc. Neurosc. Abst, 25:1151, 1999. → pages 94 [82] T. Irving, Y. Wu, T. Bekyarova, G. Farman, N. Fukuda, and H. Granzier. Thick-filament strain and interfilament spacing in passive muscle: effect of titin-based passive tension. Biophysical Journal, 100(6):1499–1508, 2011. → pages 48  141  [83] L. Itti. Realistic avatar eye and head animation using a neurobiological model of visual attention. Technical report, DTIC Document, 2003. → pages 104 [84] L. Itti. Quantitative modelling of perceptual salience at human eye position. Visual cognition, 14(4-8):959–984, 2006. → pages 115 [85] R. Jarosch. Muscle force arises by actin filament rotation and torque in the Z-filaments. Biochemical and Biophysical Research Communications, 270 (3):677–682, 2000. → pages 55 [86] B. Jewell and D. Wilkie. An analysis of the mechanical components in frog’s striated muscle. J Physiol, 143(3):515–540, 1958. → pages 17 [87] R. Johansson, G. Westling, A. B¨ackstr¨om, and J. Flanagan. Eye–hand coordination in object manipulation. the Journal of Neuroscience, 21(17): 6917–6932, 2001. → pages 102, 130 [88] W. Johnson, D. Jindrich, R. Roy, and V. Reggie Edgerton. A three-dimensional model of the rat hindlimb: Musculoskeletal geometry and muscle moment arms. Journal of biomechanics, 41(3):610–619, 2008. → pages v, 65, 78, 90, 92 [89] R. Josephson. Mechanical power output from striated muscle during cyclic contraction. Journal of Experimental Biology, 114(1):493–512, 1985. → pages 33 [90] V. Joumaa, T. Leonard, and W. Herzog. Residual force enhancement in myofibrils and sarcomeres. Proceedings of the Royal Society B: Biological Sciences, 275(1641):1411–1419, 2008. → pages 36 [91] V. Joumaa, D. Rassier, T. Leonard, and W. Herzog. The origin of passive force enhancement in skeletal muscle. American Journal of Physiology-Cell Physiology, 294(1):C74–C78, 2008. → pages 43 [92] G. Joyce and P. Rack. Isotonic lengthening and shortening movements of cat soleus muscle. The Journal of physiology, 204(2):475, 1969. → pages 19 [93] W. Kargo and L. Rome. Functional morphology of proximal hindlimb muscles in the frog Rana pipiens. Journal of Experimental Biology, 205 (14):1987, 2002. → pages 65  142  [94] W. Kargo, F. Nelson, and L. Rome. Jumping in frogs: assessing the design of the skeletal system by anatomically realistic modeling and forward dynamic simulation. Journal of Experimental Biology, 205(12):1683, 2002. → pages 65 [95] M. Kawai and H. Halvorson. Force transients and minimum cross-bridge models in muscular contraction. Journal of muscle research and cell motility, 28(7):371–395, 2007. → pages 34 [96] M. Kellermayer and H. Granzier. Calcium-dependent inhibition of in vitro thin-filament motility by native titin. FEBS letters, 380(3):281–286, 1996. → pages 43 [97] S. Kim, M. Spenko, S. Trujillo, B. Heyneman, D. Santos, and M. Cutkosky. Smooth vertical surface climbing with directional adhesion. IEEE Transactions on Robotics, 24(1):65–74, 2008. → pages 4 [98] D. Klein, A. Patino, and M. Tresch. Flexibility of motor pattern generation across stimulation conditions by the neonatal rat spinal cord. Journal of neurophysiology, 103(3):1580, 2010. → pages 95 [99] J. Kober, M. Glisson, and M. Mistry. Playing catch and juggling with a humanoid robot. In 2006 6th IEEE-RAS International Conference on Humanoid Robots. IEEE, 2012. → pages 105 [100] M. Kruger and W. Linke. The giant protein titin: a regulatory node that integrates myocyte signaling pathways. Journal of Biological Chemistry, 2011. → pages 54 [101] J. Kutch, A. Kuo, and W. Rymer. Extraction of individual muscle mechanical action from endpoint force. Journal of Neurophysiology, 103 (6):3535, 2010. → pages 97 [102] D. Labeit, K. Watanabe, C. Witt, H. Fujita, Y. Wu, S. Lahmers, T. Funck, S. Labeit, and H. Granzier. Calcium-dependent molecular spring elements in the giant protein titin. Proceedings of the National Academy of Sciences, 100(23):13716, 2003. → pages 43 [103] S. Labeit, B. Kolmerer, and W. Linke. The giant protein titin emerging roles in physiology and pathophysiology. Circulation Research, 80(2): 290–294, 1997. → pages 41  143  [104] J. Lanman, E. Bizzi, and J. Allum. The coordination of eye and head movement during smooth pursuit. Brain Research, 153(1):39–53, 1978. → pages 117 [105] A. Lappin, J. Monroy, J. Pilarski, E. Zepnewski, D. Pierotti, and K. Nishikawa. Storage and recovery of elastic potential energy powers ballistic prey capture in toads. Journal of Experimental Biology, 209(13): 2535, 2006. → pages 22, 26, 39, 62 [106] J. Lawrence 3rd, T. Nichols, and A. English. Cat hindlimb muscles exert substantial torques outside the sagittal plane. Journal of neurophysiology, 69(1):282, 1993. → pages 91 [107] S. Lee and D. Terzopoulos. Heads up!: biomechanical modeling and neuromuscular control of the neck. In ACM Transactions on Graphics (TOG), volume 25, pages 1188–1198. ACM, 2006. → pages 10, 36, 104 [108] S. Lee, J. Badler, and N. Badler. Eyes alive. In ACM Transactions on Graphics (TOG), volume 21, pages 637–644. ACM, 2002. → pages 104 [109] S. Lee, H. Chen, J. Towles, and D. Kamper. Estimation of the effective static moment arms of the tendons in the index finger extensor mechanism. Journal of biomechanics, 41(7):1567–1573, 2008. → pages 66, 94 [110] S. Lee, E. Sifakis, and D. Terzopoulos. Comprehensive biomechanical modeling and simulation of the upper body. ACM Transactions on Graphics (TOG), 28(4):99, 2009. → pages 128 [111] R. Leigh and D. Zee. The neurology of eye movements. Number 55. Oxford Univ Pr, 1999. → pages 113, 115, 117 [112] T. Leonard and W. Herzog. Regulation of muscle force in the absence of actin-myosin-based cross-bridge interaction. American Journal of Physiology-Cell Physiology, 299(1):C14–C20, 2010. → pages 44 [113] O. Lilienthal. Der vogelflug als grundlage der fliegekunst: Ein beitrag zur systematik der flugtechnik. R. Gaertner, 1889. → pages 3 [114] S. Lindstedt, P. LaStayo, T. Reich, et al. When active muscles lengthen: properties and consequences of eccentric contractions. News in physiological sciences, 16:256–261, 2001. → pages 39 [115] W. Linke, M. Ivemeyer, P. Mundel, M. Stockmeier, and B. Kolmerer. Nature of PEVK-titin elasticity in skeletal muscle. Proceedings of the 144  National Academy of Sciences of the United States of America, 95(14): 8052, 1998. → pages 41, 49, 53, 54 [116] C. Liu. Dextrous manipulation from a grasping pose. In ACM Transactions on Graphics (TOG), volume 28, page 59. ACM, 2009. → pages 104 [117] E. Loeb, S. Giszter, P. Saltiel, E. Bizzi, and F. Mussa-Ivaldi. Output units of motor behavior: an experimental and modeling study. Journal of Cognitive Neuroscience, 12(1):78–97, 2000. → pages 66, 72, 91, 94 [118] H. Maas and T. Sandercock. Are skeletal muscles independent actuators? Force transmission from soleus muscle in the cat. Journal of Applied Physiology, 104(6):1557, 2008. → pages 94 [119] K. Maruyama. Connectin, an elastic protein from myofibrils. Journal of Biochemistry, 80(2):405–407, 1976. → pages 41 [120] J. McIntyre, M. Zago, A. Berthoz, F. Lacquaniti, et al. Does the brain model newton’s laws? Nature Neuroscience, 4(7):693–694, 2001. → pages 109 [121] J. McKay, T. Burkholder, and L. Ting. Biomechanical capabilities influence postural control strategies in the cat hindlimb. Journal of biomechanics, 40(10):2254–2260, 2007. → pages 65 [122] S. McKee. A local mechanism for differential velocity detection. Vision Research, 21(4):491–500, 1981. → pages 108 [123] T. McMahon. Muscles, reflexes, and locomotion. Princeton Univ Pr, 1984. → pages 12, 14, 17, 38 [124] C. Meyer, A. Lasker, and D. Robinson. The upper limit of human smooth pursuit velocity. Vision Research, 25(4):561–563, 1985. → pages 116 [125] R. Morgan. Actin rotates as myosin translates. Journal of theoretical biology, 67(4):769, 1977. → pages 55 [126] R. Murray, Z. Li, and S. Sastry. A mathematical introduction to robotic manipulation. CRC, 1994. → pages 73 [127] F. Mussa-Ivaldi and N. Hogan. Integrable solutions of kinematic redundancy via impedance control. The International Journal of Robotics Research, 10(5):481, 1991. → pages 75  145  [128] S. Nicolopoulos-Stournaras and J. Iles. Hindlimb muscle activity during locomotion in the rat (Rattus norvegicus)(Rodentia: Muridae). Journal of Zoology, 203(3):427–440, 1984. → pages 93 [129] K. Nishikawa, A. Biewener, P. Aerts, A. Ahn, H. Chiel, M. Daley, T. Daniel, M. Hale, T. Hedrick, A. Lappin, et al. Neuromechanics: an integrative approach for understanding motor control. Integrative and Comparative Biology, 47(1):16–54, 2007. → pages 133 [130] K. Nishikawa, J. Monroy, T. Uyeno, S. Yeo, D. Pai, and S. Lindstedt. Is titin a winding filament? a new twist on muscle contraction. Proceedings of the Royal Society B: Biological Sciences, 2011. → pages 12, 40, 42, 46 [131] T. Nishizaka, T. Yagi, Y. Tanaka, and S. Ishiwata. Right-handed rotation of an actin filament in an in vitro motile system. Nature, 1993. → pages 45 [132] K. Novak, L. Miller, and J. Houk. The use of overlapping submovements in the control of rapid hand movements. Experimental Brain Research, 144 (3):351–364, 2002. → pages 102, 105, 121 [133] J. Orban de Xivry and P. Lef`evre. Saccades and pursuit: two outcomes of a single sensorimotor process. The Journal of Physiology, 584(1):11–23, 2007. → pages 113, 130 [134] F. Parsons. 4. Myology of Rodents. Part II. An Account of the Myology of the Myomorpha, together with a Com parison of the Muscles of the various Suborders of Rodents. In Proceedings of the Zoological Society of London, volume 64, pages 159–192. John Wiley & Sons, 1896. → pages 93 [135] I. Pavlov, R. Novinger, and D. Rassier. The mechanical behavior of individual sarcomeres of myofibrils isolated from rabbit psoas muscle. American Journal of Physiology-Cell Physiology, 297(5):C1211–C1219, 2009. → pages 55, 62 [136] J. Pearlman, S. Roach, and F. Valero-Cuevas. The fundamental thumb-tip force vectors produced by the muscles of the thumb. Journal of Orthopaedic Research, 22(2):306–312, 2004. → pages 66, 94 [137] C. Pelachaud and M. Bilvi. Modelling gaze behavior for conversational agents. In Intelligent Virtual Agents, pages 93–100. Springer, 2003. → pages 104  146  [138] C. Peters and A. Qureshi. A head movement propensity model for animating gaze shifts and blinks of virtual characters. Computers & Graphics, 34(6):677–687, 2010. → pages 104 [139] R. Pfeifer, M. Lungarella, and F. Iida. Self-organization, embodiment, and biologically inspired robotics. Science, 318(5853):1088, 2007. → pages 4 [140] G. Pinniger and A. Cresswell. Residual force enhancement after lengthening is present during submaximal plantar flexion and dorsiflexion actions in humans. Journal of Applied Physiology, 102(1):18–25, 2007. → pages 36 [141] R. Plamondon. A kinematic theory of rapid human movements. Biological Cybernetics, 72(4):295–307, 1995. → pages 112 [142] N. Pollard and V. Zordan. Physically based grasping control from example. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 311–318. ACM, 2005. → pages 104 [143] C. Pratt, C. Chanaud, and G. Loeb. Functionally complex muscles of the cat hindlimb. IV. Intramuscular distribution of movement command signals and cutaneous reflexes in broad, bifunctional thigh muscles. Exp. Brain Res, 85(2):281, 1991. → pages 95 [144] P. Rack and D. Westbury. The effects of length and stimulus rate on tension in the isometric cat soleus muscle. The Journal of physiology, 204(2):443, 1969. → pages 22 [145] D. Rassier. Muscle Biophysics: From Molecules to Cells, volume 682. Springer Verlag, 2010. → pages 34 [146] D. Rassier, B. MacIntosh, and W. Herzog. Length dependence of active force production in skeletal muscle. Journal of Applied Physiology, 86(5): 1445, 1999. → pages 79 [147] D. Rassier, W. Herzog, J. Wakeling, and D. Syme. Stretch-induced, steady-state force enhancement in single skeletal muscle fibers exceeds the isometric force at optimum fiber length. Journal of biomechanics, 36(9): 1309–1316, 2003. → pages 36 [148] M. Riley and C. Atkeson. Robot catching: Towards engaging human-humanoid interaction. Autonomous Robots, 12(1):119–128, 2002. → pages 105 147  [149] A. Rizzi and D. Koditschek. Progress in spatial robot juggling. In Robotics and Automation, 1992. Proceedings., 1992 IEEE International Conference on, pages 775–780. IEEE, 1992. → pages 105 [150] T. Roberts and E. Azizi. Flexible mechanisms: the diverse roles of biological springs in vertebrate movement. Journal of Experimental Biology, 214(3):353–361, 2011. → pages 39 [151] D. Robinson. The mechanics of human smooth pursuit eye movement. The Journal of Physiology, 180(3):569, 1965. → pages 115 [152] D. Robinson, J. Gordon, and S. Gordon. A model of the smooth pursuit eye movement system. Biological Cybernetics, 55(1):43–57, 1986. → pages 113 [153] C. Rode, T. Siebert, and R. Blickhan. Titin-induced force enhancement and force depression: A ‘sticky-spring’mechanism in muscle contractions? Journal of theoretical biology, 259(2):350–360, 2009. → pages 43 [154] B. Rohrer and N. Hogan. Avoiding spurious submovement decompositions: a globally optimal algorithm. Biological cybernetics, 89 (3):190–199, 2003. → pages 113 [155] M. Sagar, D. Bullivant, G. Mallinson, and P. Hunter. A virtual environment and model of the eye for surgical simulation. In Proceedings of the 21st annual conference on Computer graphics and interactive techniques, pages 205–212. ACM, 1994. → pages 4 [156] T. Sandercock. Summation of motor unit force in passive and active muscle. Exercise and sport sciences reviews, 33(2):76, 2005. → pages 94 [157] T. Sandercock and C. Heckman. Force from cat soleus muscle during imposed locomotor-like movements: experimental data versus Hill-type model predictions. Journal of neurophysiology, 77(3):1538, 1997. → pages 33, 62 [158] V. Santos and F. Valero-Cuevas. Anatomical variability naturally leads to multimodal distributions of Denavit-Hartenberg parameters for the human thumb. In Proc. 25th Annual Intl Conf. IEEE EMBS, 2003. → pages 65, 94 [159] V. Santos, C. Bustamante, and F. Valero-Cuevas. Improving the Fitness of High-Dimensional Biomechanical Models via Data-Driven Stochastic Exploration. IEEE transactions on biomedical engineering, 56(3), 2009. → pages 65, 98, 132 148  [160] S. Schaal and C. Atkeson. Robot juggling: Implementation of memory-based learning. Control Systems Magazine, IEEE, 14(1):57–71, 1994. → pages 105 [161] W. Shao and D. Terzopoulos. Autonomous pedestrians. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 19–28. ACM, 2005. → pages 104 [162] D. Shi, A. Somlyo, A. Somlyo, and Z. Shao. Visualizing filamentous actin on lipid bilayers by atomic force microscopy in solution. Journal of Microscopy, 201(3):377–382, 2001. → pages 48 [163] J. Soechting and C. Terzuolo. Organization of arm movements. motion is segmented. Neuroscience, 23(1):39–51, 1987. → pages 102, 105 [164] J. Starkes, W. Helsen, and D. Elliott. A menage a trois: the eye, the hand and on-line processing. Journal of sports sciences, 20(3):217–224, 2002. → pages 102 [165] D. Sternad and S. Schaal. Segmentation of endpoint trajectories does not imply segmented control. Experimental Brain Research, 124(1):118–136, 1999. → pages 102, 105 [166] S. Sueda and D. Pai. Dynamic hand simulation with strands. unpublished. → pages 4 [167] S. Sueda, A. Kaufman, and D. K. Pai. Musculotendon simulation for hand animation. ACM Trans. Graph. (Proc. SIGGRAPH), 27(3):83:1–83:8, 2008. → pages 64, 128, 132 [168] H. Sugi, T. Akimoto, T. Kobayashi, S. Suzuki, and M. Shimada. Possible contribution of titin filaments to the compliant series elastic component in horseshoe crab skeletal muscle fibers. Advances in experimental medicine and biology, 481:371–382, 2000. → pages 39, 62 [169] H. Syyong, A. Raqeeb, P. Par´e, and C. Seow. Time course of isotonic shortening and the underlying contraction mechanism in airway smooth muscle. Journal of Applied Physiology, 111(3):642–656, 2011. → pages 19 [170] I. Telley, J. Denoth, E. St¨ussi, G. Pfitzer, and R. Stehle. Half-sarcomere dynamics in myofibrils during activation and relaxation studied by tracking fluorescent markers. Biophysical journal, 90(2):514–530, 2006. → pages 62, 63 149  [171] I. Telley, R. Stehle, K. Ranatunga, G. Pfitzer, E. St¨ussi, and J. Denoth. Dynamic behaviour of half-sarcomeres during and after stretch in activated rabbit psoas myofibrils: sarcomere asymmetry but no sarcomere popping. The Journal of Physiology, 573(1):173, 2006. → pages 36 [172] J. Teran, S. Blemker, V. Hing, and R. Fedkiw. Finite volume methods for the simulation of skeletal muscle. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 68–74. Eurographics Association Aire-la-Ville, Switzerland, Switzerland, 2003. → pages 64 [173] D. Terzopoulos and T. Rabie. Animat vision: Active vision in artificial animals. Videre. Journal of Computer Vision Research, 1(1):2–19, 1997. → pages 104 [174] J. Tobin. To conquer the air: The Wright brothers and the great race for flight. Free Press, 2003. → pages 3 [175] E. Topp and J. Rosen. An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures. IEEE Transactions on Biomedical Engineering, 37(8):157, 1990. → pages 65 [176] M. Tresch and E. Bizzi. Responses to spinal microstimulation in the chronically spinalized rat and their relationship to spinal systems activated by low threshold cutaneous stimulation. Experimental Brain Research, 129 (3):401–416, 1999. → pages 68 [177] W. Tsang, K. Singh, and E. Fiume. Helping hand: an anatomically accurate inverse dynamics solution for unconstrained hand motion. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 319–328. ACM, 2005. → pages 104 [178] X. Tu and D. Terzopoulos. Artificial fishes: Physics, locomotion, perception, behavior. In Proceedings of the 21st annual conference on Computer graphics and interactive techniques, pages 43–50. ACM, 1994. → pages 4, 104 [179] R. Urtasun, D. Fleet, A. Geiger, J. Popovi´c, T. Darrell, and N. Lawrence. Topologically-constrained latent variable models. In Proceedings of the 25th international conference on Machine learning, pages 1080–1087. ACM, 2008. → pages 123  150  [180] F. Valero-Cuevas, F. Zajac, C. Burgar, et al. Large index-fingertip forces are produced by subject-independent patterns of muscle excitation. Journal of biomechanics, 31(8):693–704, 1998. → pages 133 [181] F. Valero-Cuevas, J. Towles, and V. Hentz. Quantification of fingertip force reduction in the forefinger following simulated paralysis of extensor and intrinsic muscles. Journal of Biomechanics, 33(12):1601–1609, 2000. → pages 66, 94 [182] F. Valero-Cuevas, M. Johanson, and J. Towles. Towards a realistic biomechanical model of the thumb: the choice of kinematic description may be more critical than the solution method or the variability/uncertainty of musculoskeletal parameters. Journal of Biomechanics, 36(7): 1019–1030, 2003. → pages 65, 94 [183] F. Valero-Cuevas, H. Hoffmann, M. Kurse, J. Kutch, and E. Theodorou. Computational models for neuromuscular function. Biomedical Engineering, IEEE Reviews in, 2:110–135, 2009. → pages 133 [184] A. Vallbo and J. Wessberg. Organization of motor output in slow finger movements in man. The Journal of physiology, 469(1):673, 1993. → pages 101, 102, 105 [185] S. Walcott and W. Herzog. Modeling residual force enhancement with generic cross-bridge models. Mathematical biosciences, 216(2):172–186, 2008. → pages 37 [186] J. Wang, D. Fleet, and A. Hertzmann. Multifactor gaussian process models for style-content separation. In Proceedings of the 24th international conference on Machine learning, pages 975–982. ACM, 2007. → pages 123 [187] J. Wang, S. Hamner, S. Delp, and V. Koltun. Optimizing locomotion controllers using biologically-based actuators and objectives. ACM Transactions on Graphics (TOG), 31(4):25, 2012. → pages 4, 10 [188] K. Wang, J. McClure, and A. Tu. Titin: major myofibrillar components of striated muscle. Proceedings of the National Academy of Sciences, 76(8): 3698, 1979. → pages 41 [189] Q. Wei, A. Jarc, S. Yeo, T. Sandercock, M. Tresch, and D. Pai. A biomechanicsl model of the rat hindlimb: modeling with strands and subject-specific registration. → pages 132 151  [190] Q. Wei, S. Sueda, and D. Pai. Physically-based modeling and simulation of extraocular muscles. Progress in biophysics and molecular biology, 103 (2):273–283, 2010. → pages 4 [191] C. Williams, M. Regnier, and T. Daniel. Axial and radial forces of cross-bridges depend on lattice spacing. PLoS computational biology, 6 (12):e1001018, 2010. → pages 45 [192] P. Williams and G. Goldspink. Changes in sarcomere length and physiological properties in immobilized muscle. Journal of Anatomy, 127 (Pt 3):459, 1978. → pages 29 [193] T. Williams. A new model for force generation by skeletal muscle, incorporating work-dependent deactivation. The Journal of Experimental Biology, 213(4):643–650, 2010. → pages 33 [194] T. Winters, M. Takahashi, R. Lieber, and S. Ward. Whole muscle length-tension relationships are accurately modeled as scaled sarcomeres in rabbit hindlimb muscles. Journal of biomechanics, 44(1):109–115, 2011. → pages 15, 25 [195] B. Wohlfart and K. Edman. Rectangular hyperbola fitted to muscle force-velocity data using three-dimensional regression analysis. Experimental physiology, 79(2):235–239, 1994. → pages 25 [196] D. Wolpert, Z. Ghahramani, and M. Jordan. An internal model for sensorimotor integration. Science, 269(5232):1880, 1995. → pages 106 [197] R. Woodworth. Accuracy of voluntary movement. The Psychological Review: Monograph Supplements, 3(3):i, 1899. → pages 101 [198] K. Yamane, J. Kuffner, and J. Hodgins. Synthesizing animations of human manipulation tasks. In ACM Transactions on Graphics (TOG), volume 23, pages 532–539. ACM, 2004. → pages 104 [199] K. Yin, M. Cline, and D. Pai. Motion perturbation based on simple neuromotor control models. In Computer Graphics and Applications, 2003. Proceedings. 11th Pacific Conference on, pages 445–449. IEEE, 2003. → pages 4 [200] L. Young and L. Stark. Variable feedback experiments testing a sampled data model for eye tracking movements. Human Factors in Electronics, IEEE Transactions on, (1):38–51, 1963. → pages 113, 115 152  [201] R. Young, S. Scott, and G. Loeb. The distal hindlimb musculature of the cat: multiaxis moment arms at the ankle joint. Experimental Brain Research, 96(1):141–151, 1983. → pages 91 [202] M. Zago, J. McIntyre, P. Senot, and F. Lacquaniti. Visuo-motor coordination and internal models for object interception. Experimental Brain Research, 192(4):571–604, 2009. → pages 105 [203] F. Zajac. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical reviews in biomedical engineering, 17(4):359, 1989. → pages 7, 14, 19, 32, 96  153  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items