4.9 SIMULATION VAN DEN BOGERT, A.J. NIGG, B.M. 4.9.1 INTRODUCTION Research on human movement does not always use human subjects. In most areas of biomedical research, animal models are used when rigorous control of experimental condi- tions is required, measurements require invasive techniques or when the use of human sub- jects is not ethically acceptable. In the context of biomechanics of human movement, animal models are typically used to study basic biological mechanisms. A good example is researcg on muscle mechanics (see chapter 2.7), where full experimental control over mus- cle length and activation was required and direct measurements of variables such as mus- cle force, muscle fibre length, and even heat production are done. This approach would not be possible in human subjects. For more applied questions, related to specific human inju- ries or specific human movement tasks, animal models are not appropriate because the dynamics and anatomy are not comparable. For instance, the study of knee injuries in dis- tance runners has a basic biology component and a mechanical component. Basic biology or in this case reaction of cartilage to mechanical loads, can be studied in animal models. It is, however, obviously impossible to use animal models to learn how to reduce forces in the knee joint with appropriate orthotics or other mechanical interventions. In human sub- jects, the force in the joint can not be measured, other factors contributing to injury can not be controlled, and ethical problems arise in the study of more acute, serious injuries. Anal- ogous to the use of animal models, it is therefore desirable to study such movement-related questions with methods that allow full experimental control as well as access to all mechanical variables for measurement. In the case of human movement, one possible option currently available is to build a model of the system under investigation, and then perform experiments on the model. Although physical/mechanical models have been used, computer models (also numerical models) are almost exclusively used for this purpose. Simulation is defined as follows: Simulation is the process of performing experiments on a numerical model. It is of crucial importance that the model reflects the causality of the system that is investigated: the input of the model should precede the output in a causative sense. Once such a numerical model has been developed, experiments can be performed quickly and at low cost. It should, however, be kept in mind that the experiment is always performed on the model only, and conclusions may not apply to the human system when the model is not valid. This is a central problem in modeling and simulation. It is important that a model replicates those features of the system that are essential for the question that is studied. This requires careful design of the model and appropriate validation tests once the model is finished. In biomechanics, simulation is applied in four different areas: dynamics of movement, tissue mechanics, fluid flow, and measuring techniques. In dynamics of movement, simu- lations allow the experimenter to apply controlled perturbations to the anatomy, muscle coordination, and control system, and observe the resulting changes in movement and forces. Tissue mechanics can be studied using finite element models (FEM), which are the numerical equivalent of partial differential equations (PDE) that describe the local rela- tionships between stress and strain in tissue (see chapter 3.7). These methods have been applied to bone, cartilage, tendon, ligament, muscle, and even brain tissue during impact. Fluid flow models, used in cardiovascular mechanics, also use finite element techniques, but use hydrodynamics rather than properties of solids. Biomechanical measuring tech- niques (chapter 3) often suffer from errors, and simulation may be used to determine how errors in the raw data are propagated to the final result of the calculation. The term "sensi- tivity analysis" is sometimes used, but when correctly done with realistically simulated statistical error distributions, the technique is known as Monte-Carlo simulation. The focus of this chapter will be the simulation of movement. 4.9.2 DIFFERENTIAL EQUATIONS A mathematical model for the simulation of movement always takes the form of a sec- ond-order system of ordinary differential equations (ODE). This term indicates that the equation contains not only the unknown movement variables (e.g., joint angles), but also their time-derivatives up to the second order (i.e., accelerations). Such equations have been derived for simple systems in sections 4.4.4 and 4.4.5 (***FIRST ED.). A one- dimensional mass-spring damper system (Fig. 4.4.9 ***FIRST ED.) was described by an equation, which is a single second-order linear ODE, with the unknown . A two-seg- ment planar pendulum (Fig. 4.4.33) was described by two coupled second-order non-lin- ear ODEs (equations (4.4.57) and (4.4.58)). These equations are non-linear because they involve non-linear functions (sin and cos) of the unknowns and . The only reason for making a distinction between linear and non-linear ODEs is that linear equa- tions can be solved analytically, which means that general solutions can be found in the form of equations giving movement explicitly as a mathematical function of time. The solutions for second-order linear ODEs take the form of combinations of exponential and/ or harmonic (sine and cosine) functions. Even moderately complex models of the human musculo-skeletal system, however, lead to non-linear equations. For non-linear equations, analytical solutions can only be found in special cases. But even without solving the equation completely, analytical methods can be used to study certain aspects of the behaviour of the system. As an example, we will consider a cyclist going downhill without pedaling. The cyclist will be considered as a particle and rolling friction will be ignored. A free body diagram (Fig. 4.9.1) reveals that there are three forces acting on the cyclist: gravity, ground reaction force, and air resistance. Considering only the component of force and movement parallel to the inclined plane, the equation of motion can be formulated as: (*equation*) From wind tunnel tests, it is known that the air drag is proportional to the square of the velocity: (*equation*) A typical value for is 0.15 Ns m-1 (de Groot, 1995). Combining the two equations, and substituting the derivative of velocity for acceleration, we obtain: (*equation*) Note that this is a first-order ODE, where the velocity is the unknown. Although we will not attempt an analytical solution, it is easy to determine a steady-state solution, if there is one. In a steady state, nothing changes and all derivatives are zero. Setting in equation (*equation*) and solving for gives the steady-state velocity: (*equation*) For a typical cyclist, with =75 kg, on a slope of 10, the steady-state velocity will, therefore, be just over 29 m s-1. Another aspect of the system behaviour that can be studied analytically is stability; by looking at small perturbations from equilibrium, and determining whether the forces respond to push the system back to equilibrium or away from it. A major advantage of ana- lytical methods is that they allow a general understanding of how the behaviour of the sys- tem depends on system parameters, such as mass or stiffness. However, during movement, we are usually interested in simulating the entire time course of the movement. Only in rare cases can this be done by analytical methods. Usually, a numerical solution of the ODE is required. This is where mathematical modelling becomes computer simulation. 4.9.3 NUMERICAL SOLUTION METHODS Every first order ODE can be formulated as: (*equation*) The function may be a simple mathematical function, as in equation (*equation*), or a function that requires a substantial computation, as we will see later for biomechanical systems. The variable is the state variable of the model. A large class of numerical solution methods for ordinary differential equations is based on Taylor's theorem, which relates derivatives of a function to a change in function value over a finite time interval. If is a sufficiently smooth function, Taylor's theorem states: (*equation*) where is the first derivative of , is the th derivative, and is an error term proportional to . A th order ODE solver is obtained by truncating the right hand side (the Taylor series) after the th term, and choosing small enough that the error term is small enough to be neglected. The simplest method is the forward Euler method, which is a first order method: (*equation*) If the state, e.g., position or velocity, of the system at time is known, the right hand side of the differential equation (*equation*) can be calculated, multiplied by a time step , and added to the previous state to obtain the new state at time . This is repeated again and again until reaches the end of the time interval that we are inter- ested in. This solution method can also be applied to a system of simultaneous differential equations, as shown later. Note that this solution process requires that the initial state at =0 is known. Results of a simulation always depend on this initial condition. The steady- state solution, however, usually does not depend on the initial condition, because damping in the system damps out those movements that depend on the initial state. The Euler method is the simplest method, but has a larger error for the same step size than higher order methods, where the error is proportional to a higher power of the step size (remember that is very small). This means that to obtain the same accuracy, the Euler method must use a smaller step size than some other methods, and the function needs to be calculated more frequently. However, the forward Euler method is adequate for work where computer speed is not a limiting factor and will be used in the examples in this chapter. Details of other methods may be found in textbooks on numerical methods (e.g., Shampine and Gordon, 1975; Press et al., 1992). Matlab provides a convenient pro- gramming environment for mathematical modelling and is available at low cost. Although Matlab has built-in ODE solvers, we will use the forward Euler method in all examples to make it easier to translate programs to other programming languages. As an example, the Matlab program in Listing 1 applies the forward Euler method to the differential equation (*equation*) describing the downhill coasting cyclist. The output is a graph of velocity versus time (Fig. 4.9.2). As an example of an experiment, a second simu- lation is done to simulate the effect of a 1000 m altitude where air density is approximately 10% lower. The choice of step size has an effect on accuracy, as can be seen from Taylor's theo- rem. This is demonstrated in Fig. 4.9.3, showing the effect of step size on the results. Accuracy obviously suffers with increasing step size. The solution at =20 s no longer converges to the steady state, which indicates that the numerical method has become unstable. This occurs for the forward Euler method when becomes larger than twice the smallest time-constant of the system. A time-constant can be roughly defined, for a single ODE, as the ratio between the simulated variable and its derivative. When the solution explodes, as in the example with =20 s, this is a clear indication that should be reduced. An even lower value of is necessary to achieve a good accuracy. Commercial software as well as some of the computer implementations in Shampine and Gordon (1975) automatically adjust the step size to achieve a certain, user-specified accuracy. When using a method with fixed step size, as in Listing 1, one can test the accuracy of the solution simply be running the simulation again with half the step size. If the difference between the two sets of results is smaller than the desired accuracy, the step size was suffi- ciently small. Otherwise, reduce again by a factor two and repeat the test. For this model, the step size of 0.01 s, which was chosen initially, is clearly small enough to guar- antee sufficient accuracy. Complex systems typically require a set of simultaneous ODEs: (*equation*) Once the state, which now consists of all state variables , is known, all right- hand sides can be computed, and a forward time step can be made by applying a numerical solu- tion method sequentially to each variable. Equations of motion for mechanical systems generally lead to a set of second-order ODEs, but can be transformed into a set of first-order ODEs so that standard numerical methods can be applied. If the original set of equations is: (*equation*) we can transform this to a new set of variables , with: (*equation*) which transforms equation (*equation*) to: (*equation*) This is an equation of the form (*equation*) that can be solved using the existing numerical methods. In systems of simultaneous ODEs, there are as many time constants as state variables, and it is the smallest time constant that determines whether the numerical method is stable. The largest time constant, however, indicates what the duration of the simulation should be in order to examine the complete system dynamics. When the ratio between the smallest and largest time constant is very small, we consider the ODE to be stiff. An example of a stiff system is a simulation of a walking horse, with a largest time constant representing the stride cycle of about 1 s, while the smallest time constant is about 0.1 ms, which is due to the interaction between hoof and ground (Bogert et al., 1989). Using a simple numerical method, such as described above, time steps would have to be extremely small compared to the duration of the simulation, which leads to long computation times. Unfortunately, this problem becomes worse with higher order methods. A special class of methods has, therefore, been developed for stiff problems (Gear, 1971). The simplest of these is the backward Euler method. In order to make one time step in solving equation (*equation*), is solved from: (*equation*) Since the function is usually a complex non-linear function, this requires iterative solution, making one time step costly in terms of computation time. However, the method is stable for step sizes much larger than the smallest time constant thus resulting in better overall performance for stiff systems. 4.9.4 EQUATIONS OF MOTION FOR MECHANICAL SYSTEMS In section 4.4.5 (***FIRST ED.), a straightforward method for formulation of equa- tions of motion for multi-segment systems was introduced. This method is based on the Newton-Euler equations for each body segment. In two dimensions, these equations are: (*equation*) (PLEASE GIVE THESE EQUATIONS SEPARATE LABELS: ...a and ...b) where is the position vector of the center of mass of the segment in the global reference frame, is the orientation of the body, is the mass, and is the moment of inertia of the segment. The forces and moments acting on the body, are identified using a free body diagram. Note that these equations are the same equations that can be used to solve unknown forces and moments from measured movements. That type of analysis, inverse dynamic analysis, is not considered to be a simulation, since the input (movement and ground reaction force data) is not causative to the output (internal forces and moments). Proper simulation of movement is only possible using forward dynamic analysis, where the inputs are forces that cause the movement. From equation (*equation*) we see that a system with body segments has degrees of freedom, since each body segment has a two- dimensional position vector , and one orientation angle . The next step is to take into account the kinematic connections (joints) between the segments, which reduce the number of degrees of freedom. Identify a suitable set of variables as the degrees of freedom (also referred to as generalized coordinates of the system), and eliminate all other kinematic variables from the equations of motion. On p. 478 (***FIRST ED.) this proce- dure was followed to derive equations of motion for a two-segment model of the arm. Although straightforward, the method is tedious and leads to large algebraic equations even for systems of moderate complexity. Once derived, the equations of motion can be easily simulated using a numerical method, as shown in Listing 2 for the two-segment sys- tem. Note that this system has two degrees of freedom, and that the equations of motion have been formulated so that the two segment orientations can be solved as a function of time. With two degrees of freedom, we obtain four first-order ODEs. Hence, the forward Euler equation appears four times in the program, once for each angle and once for each angular velocity; the state variables of the system. Results of this simulation are shown in Fig. 4.9.4. The two segments had initial angular velocities of zero and 10 rad/s, respec- tively. The simulation shows that the angular momentum is transferred back and forth between the two segments, in a 300 ms cycle. Note that segment 2 moves briefly in the opposite direction during each cycle. The whole system also rotates around the shoulder joint at approximately one revolution (360) every two seconds. Since the Newton-Euler method is cumbersome and difficult to perform systemati- cally, alternative methods are often used. The oldest of these is Lagrange's equation. After identifying a set of generalized coordinates , one uses the kinematic con- nections to formulate the potential energy of the system (which includes gravity and elastic forces) as a function of , and the kinetic energy as a function of and . The Lagrangian is defined as the difference and Lagrange's equations are: (*equation*) from which the generalized accelerations are easily solved to obtain a system of simultaneous second order ODEs, which are identical to what would have been obtained using the Newton-Euler formulation. The generalized forces are obtained from the original forces and moments by applying the principle of virtual work (reference PLEASE CLARIFY). This procedure requires considerable algebraic manipulation, but the proce- dure is more systematic than derivations from the Newton-Euler equations. The algrabraic manipulations can even be automated using symbolic algebra software such as Mathemat- ica or Maple. Even more suited for automated derivation of equations of motion are Kane's equations (Kane and Levinson, 1985). This method is used in the commercial software SD/FAST and AUTOLEV. SD/FAST is also used as the dynamics engine for SIMM and Pro/Mechanica. Both of these packages provide a comprehensive user-interface to help create models, numerical solution, and visualizing results of simulations. Another class of methods is based on the use of Newton-Euler equations with addi- tional constraint equations. Since the constraint equations are algebraic equations which express relationships between positions of two segments, this leads to a mixed set of dif- ferential and algebraic equations (DAE; Haug, 1989). This method for formulating the equations, together with a numerical solution method, has been implemented in the com- mercial products DADS and ADAMS, both of which include a comprehensive user interface. An advantage of DAE methods is that the reaction forces associated with the kinematic constraints (i.e., the joint contact forces in a musculo-skeletal model) are com- puted during the simulation. An advantage of Lagrange's equations is that it is easy to quantify the contributions of each driving force, such as a muscle, to the accelerations induced in each degree of freedom (Zajac and Gordon, 1989). For instance, it can be seen that the hamstrings during stance, function as a knee extensor because the large upper body mass effectively couples knee and hip extension (Andrews, 1985). It seems, there- fore, that DAE methods are more suitable for injury-related questions, while methods using generalized coordinates more readily provide insight into the control aspects of movement. All commercial packages mentioned here allow modelling and simulation in three dimensions as well as in two dimensions. 4.9.5 MUSCLE MODELS Movement simulations can be driven using a variety of methods. In simulations of human movement during free fall or space flight, movement of the body can be produced by well-coordinated limb movements (Pasarello and Huston, 1971). Since the required muscle forces are extremely small it is appropriate to consider joint rotations as the cause of the movement (Yeadon, 1993c). The required muscle forces may be computed, but are not relevant to the problem as long as they do not exceed the force-generating capacity of the muscles. Limb orientation variables are also closer to real-world applications. When gravitational forces and ground contact are both present, joints can no longer be moved arbitrarily and muscle forces are a more suitable input for movement simulation. These can be simply represented by joint torques as a function of time (e.g. Hubbard, 1980). Although simple and effective, this approach has limitations. First, when experimenting with joint torque inputs, one is never certain that the value used is within the capability of human muscles. This is especially problematic when optimizing performance, requiring addition of constraints to prevent unphysiological solutions (see later section on optimiza- tion). Joint torque input also ignores certain mechanisms which are important for control of movement, especially movements of longer duration. When a muscle is stretched, its force output increases due to the mechanical properties. This provides shock absorption, stabilization, and damping of movement. Also, muscles often span two joints (biarticular muscles). This provides a mechanical coupling between joints that may facilitate control and make movement more efficient (Ingen Schenau and Bobbert, 1993). When muscle forces are included in a simulation of movement, it is usually desirable to include several well-known mechanical properties of muscle: the force-length relation- ship, the force-velocity relationship, and the activation dynamics. A muscle model which includes these properties will allow the muscle forces to respond realistically to changes in length and neural stimulation. Muscle models exist in various levels of complexity, but for movement simulation, the three-component Hill model has been used almost exclusively (Winters, 1990; Zajac, 1989). The model consists of a contractile element (CE) and two non-linear elastic elements, the parallel elastic element (PEE) and the series elastic ele- ment (SEE), as illustrated in Fig. NEW1. The CE represents the muscle fibers, the SEE represents tendon and other elastic tissue in series with the fibers, and the PEE represents the passive properties of the fibers and elastic tissue surrounding the muscle fibers. The elastic elements are described by a mathematical force-length which has the property that no force is generated below a certain length (the slack length), and stiffness increases as force increases. A simple, two-parameter model for such an element is: (*equation*) The stiffness parameter can be chosen such that the tendon strain at maximal muscle force is about 3-5%, or about half the level where damage occurs. This is supported by lit- erature data (Winters, 1990), but it should be noted that functional demands on different tendons are different. For instance, finger tendons are very long and accurate control of movement is only possible when strains remain low, compared to for instance the tendons in the lower extremity where elastic energy storage occurs. The contractile element is described by Hill's force-velocity relationship: (*equation*) Note that this is equivalent to the traditional formulation (equation 2.5.1 on p.173 FIRST ED.) since the contraction velocity is equal to minus the time-derivative of con- tractile element length . The length-dependence is usually implemented by making the isometric force dependent on the length . Theoretically, the function has a shape as in Fig. 2.5.13 (on p. 171 ***FIRST ED.), based on the properties of sliding fila- ments. However, in real muscle the function is more smooth and can be parameterized by a simple mathematical function, for instance: (*equation*) According to the sliding filament theory, the width parameter , the maximum amount of shortening or lengthening in the CE, should be equal to 0.56 for human muscle (Walker and Schrodt, 1973). Due to the effects of pennation and series arrange- ment of muscle fibres, tends to be higher for the complete CE of real muscles (Bogert et al., 1998). Equation (*equation*) is substituted in equation (*equation*) to obtain the force-length- velocity relationship, and finally multiplied by a scaling factor , with a value between 0 and 1, which represents the activation, or active state of the muscle. This results in a force- length-velocity-activation relationship for the CE: (*equation*) where the right-hand side is a known mathematical function. This model for the contrac- tile element is simple, but not entirely consistent with certain observations in isolated mus- cle preparations. For instance, equation (*equation*) predicts that the maximal shortening velocity of the CE is , and, therefore, depends on length but not on active state. A more realistic CE model, also including the eccentric (for ) force-velocity rela- tionship, is described by van Soest and Bobbert (1993). For a more detailed review of modeling aspects see Winters (1990). From Fig. NEW1 it can be deduced that the force in the SEE is equal to the sum of forces in the CE and PEE. Also, the length of the SEE is the length of the whole mus- cle-tendon complex minus . Therefore, the total state equation for the muscle is obtained as follows: (*equation*) All functions are known, so the lengthening velocity of the CE can be solved from equation(*equation*): (*equation*) This is a first-order ODE, with as state variable. In experiments on isolated muscle, as described in section 2.5.5 (***FIRST ED.), acti- vation and length are controlled by the experimenter. In that case, is the only unknown and equation (*equation*) alone is sufficient to simulate such experiments, as is shown in the example in Listing 3. Activation is assumed maximal (=1) and the imposed mus- cle length is at first isometric, followed by an isokinetic (constant velocity) ramp shortening, and finally isometric again. Note that the PEE is neglected since the muscle length is assumed to be shorter than the slack length of the PEE. The simulated force (Fig. 4.9.5) shows first the build-up of isometric force during lengthening of the SEE, then at =1 s, the force drops suddenly because the initial shortening occurs only in the SEE. The CE cannot shorten as quickly, due to the force-velocity relationship. At =2 s, the shorten- ing stops and the force stabilizes at the isometric force for the final muscle length. Simula- tions such as these are helpful when designing experiments in muscle mechanics, and to interpret results of such experiments. Muscle properties for human muscles, such as , , , , , and can not be measured directly in human muscles. Animal data may be used, after appropriate scaling to muscle size, or values derived from theoretical considerations (see section 2.5.5 FIRST EDITION). It is important, however, to compare the performance of the model to macroscopic human data such as isometric strength data, and adjust the muscle properties where needed (Pandy and Zajac, 1990; Gerritsen et al., in press). The activation level can be given as input for the muscle model, but it is more real- istic to consider the activation dynamics of muscle. Activation of a muscle requires a polarization wave to reach all parts of the muscle, which eventually results in a concentra- tion of Ca2+ ions which leads to the generation of force. The active state can be identi- fied as the concentration of calcium ions bound to troponin, relative to the maximum. This concentration can not change instantaneously since it is the result of electrical and chemi- cal processes. To simulate this, one can use a model that describes how active state depends on the neural stimulation signal . A simple linear model for activation dynam- ics is: (*equation*) which is a first order ODE with state variable . When is larger than , the active state rises and approaches the value asymptotically with a time constant . The oppo- site happens when drops below . In a slightly more complex version (He et al., 1991), the model takes into account that activation occurs at a higher speed than deactiva- tion: (*equation*) where the activation time constant is usually assumed to be shorter than the deactivation time constant (Winters and Stark, 1987). Other, more complex, activation models have been used, ranging from linear second-order models (e.g., Winters and Stark, 1987) to a non-linear model which includes separate differential equations for motor unit recruitment and stimulation as well as the effect of fibre shortening on the acti- vation process (Hatze, 1981). If a model of activation dynamics is not included, optimiza- tions of a musculoskeletal model (see later section on optimization) may result in It should be kept in mind that this three-component Hill model is based on a limited set of experimental observations, such as isometric, isotonic, and isokinetic contractions, usu- ally performed at maximal activation. Interactions between the effects of length change and activation may therefore be poorly predicted. Also, history-dependent effects such as stretch potentation are not included in the model. The Hill model is half mechanical (PEE and SEE), and half empirical (CE) so it may not be valid when the CE is in a state that did not occur in the experiments on which the model was based. More complex models, based on the cross-bridge mechanism of muscle contraction have been formulated as a set of coupled ODE's (Zahalak, 1986). Such models should be valid for a wider range of condi- tions but require more complex programming, more computation time, and most impor- tantly, also seem to produce certain results that disagree with certain observations (Cole et al., 1995). Cross-bridge models have hardly been applied for movement simulation and the Hill model, in spite of its limitations, is presently considered the most practical muscle model for simulation of human movement. When muscles are connected to a skeleton with many degrees of freedom, the muscle models must be coupled to the equations of motion, resulting in a combined system of equations representing the musculo-skeletal system. This will be discussed in the next sec- tion. 4.9.6 SIMULATION USING MUSCULO-SKELETAL MODELS The output of a muscle model is the magnitude of muscle force, a scalar . The equa- tions of motion (NUMBER OF the equation pair r-dot-dot= and phi-dot-dot=) require a force vector (magnitude and direction) and also a point of application, so that the moment can be computed. Also, the muscle model requires the instantaneous length of the muscle as input. All of these require a model of the musculo-skeletal anatomy. One method for doing this is to model the muscle as a straight line of action between origin and inser- tion (Fig. 4.9.6). If the kinematic state variables are known, the position vectors of origin and insertion can be calculated from the position and orientation of the two body segments. Then, muscle length is simply: (*equation*) and the force vector applied to the insertion point is: (*equation*) The same vector, but in opposite direction, is applied to the origin point . Straight- line muscle models have limitations, and usually require additional wrapping points to model muscles and tendons that follow a curved path over underlying tissue (Delp, 1990). Failure to properly model the path of the muscle may result in underestimation of moment arms or even put the muscle on the opposite side of the joint. For some applications, it may be more practical at this point to include only the effect of muscle forces on movement, i.e., the moments exerted with respect to the joint centres. In the Newton-Euler formulation, this removes muscle forces entirely from (NUMBER OF the first of the equa- tions of motion, so only r-dot-dot=...) so that other variables in that equation (the joint reaction forces) will no longer be correct. For questions related to control of movement, this is acceptable. For injury-related questions, however, accurate joint reaction forces are usually required. In the Lagrange formulation of the equations of motion, the generalized forces only contain the effect of muscles on movement. The contribution of a muscle force to a generalized force can be determined using the principle of virtual work. When the system undergoes a hypothetical small change in the generalized coordinates , the muscle will undergo a virtual length change . The principle of virtual work states that the work done by the muscle should be equal to the work done by the generalized forces, i.e., the dot product of the generalized force vector and the generalized displace- ment . Since muscle does positive work when shortening (negative ), this can be expressed as: (*equation*) Expanding the dot product, and considering that may be a function of all general- ized coordinates , we obtain: (*equation*) Since this must be true for all possible virtual displacements, the contribution of mus- cle force to generalized force is: (*equation*) If a generalized coordinate is a joint angle, the corresponding generalized force will be the moment of the muscle at that joint. If the generalized coordinates are segment orienta- tions, the muscle moment will contribute to the generalized forces for both segments. The partial derivatives , which are the moment arms, may be derived from an anatomical model, for instance using straight lines, or may be measured directly in cadaver specimens or in vivo using three-dimensional imaging techniques. Moment arms can be simply mea- sured as the distance from the muscle's line of action to the joint centre, but this requires the assumption of a fixed joint centre. A more elegant method is to collect data of muscle length change at various joint angles (or other generalized coordinates), and then use regression analysis to find a mathematical relationship describing as a function of the generalized coordinates. The partial derivatives may then be computed analytically. This method was first used by Grieve et al. (1978) to determine moment arms of the ankle plan- tarflexors, and applied later to determine moment arms for hip and ankle muscles (Spoor et al., 1990; Visser et al., 1991). It is important to note that the partial derivative is a moment arm in meters if and are measured in meters and radians, respectively. It is important that such regression models not be extrapolated beyond the range of joint angles in the experimental protocol. Extrapolation may result in unreliable moment arms and in such cases, assuming a simple constant moment arm may result in a better model. Summarizing, the coupling between skeleton and muscles requires three steps of cal- culation at each time step of the simulation: (1) Compute muscle lengths from the joint angles in the skeleton. (2) Compute muscle force from the state equations of the muscles. (3) Apply the muscle forces to the skeleton. For a complete musculo-skeletal model, passive forces should also be considered. Lig- aments, when their function is to guide joint kinematics, need not be included since their function has been incorporated in the joint model. For instance, the cruciate ligaments in the knee can be assumed to be inextensible during normal movements, resulting in a mov- ing centre of rotation at the intersection of the cruciates. These guiding ligaments do not perform mechanical work, and can be omitted from the equations of motion. They will, however, have an effect on the articular contact force. In order to model this, joints can no longer be considered to be kinematic mechanisms and deformation needs to be modelled. This can be done (e.g., Blankevoort and Huiskes, 1996), but this is outside the scope of this chapter. Internal passive force also occurs at the end of the range of motion of a joint. This may be the result of ligaments being stretched, muscles being compressed, or a com- bination of both. The net effect on the movement can be modelled as an extra joint moment, which is a function of one or more joint angles and suitable models can be found in the literature (e.g. Yoon and Manour, 1982). Passive tension in ligaments spanning sev- eral joints, which occurs in some animals (Ingen Schenau and Bobbert, 1993), can also be modelled as a simple passive muscle model with only a SEE and coupled to the skeleton in the same way as muscles. As an example of a combined muscle and skeletal model (a musculo-skeletal model), Listing 4 presents a simulation of a maximal effort ballistic knee extension. The muscle model from Listing 3 is used, the moment with respect to the knee joint is computed assuming a constant moment arm, and the moment is applied to the equation of motion of a single rotating body segment (the shank-foot combination). Segment properties (length, mass, moment of inertia, center of mass) were obtained from Winter (1979). The results (Fig. 4.9.6) show that the force increases at first, but then quickly drops off as the move- ment accelerates due to the effects of the force-length and the force-velocity relationships. Such interactions between skeletal dynamics and muscle dynamics can only be observed in a model in which both are represented and coupled. As an example of an experiment, the dashed line shows what happens with an increased moment arm. Although the initial acceleration is larger, as expected, the joint does not even reach full extension because the force drops off too quickly. This is an example of how simulation can be used to look at the effect of design on performance. 4.9.7 MODELING OF EXTERNAL FORCES External forces are forces acting between the system of interest and the environment. These forces must be included in the right hand side of the equations of motion (4.9.13) for the system. Common external forces are: gravity, contact forces, and aerodynamic forces. Gravity is simply modeled as a constant force vector of magnitude mg and down- ward direction, applied to the center of mass of each rigid body segment in the model, where m is the mass of the body segment and g is the acceleration of gravity (9.81 m/s2). Contact forces can often be included by modeling the contact as a kinematic connection, similar to a joint. This was, for example, appropriate for the foot-ground contact forces in a simulation of vertical jumping (Pandy and Zajac, 1990), provided that the simulation is terminated when tensile forces in the connection begin to occur. A similar model for walk- ing could therefore only be used for approximately 48% of the stride cycle, which is the time that one toe remains connected to the ground (Yamaguchi and Zajac, 1991). A model of cycling on an ergometer can use a similar method to include external forces: hinge joints between pelvis and ground, as well as between crank and ground effectively gener- ate the external forces required to keep the pelvis attached to the seat and to prevent the crankshaft from moving (Neptune and Bogert, 1998). This technique, sometimes referred to as the "hard constraint" method, becomes difficult to implement when a connection is not permanent, for instance the contact between foot and ground during simulation of a complete walking movement. At heel strike, a joint is suddenly added which requires the heel velocity to decrease instantaneously from some value to zero. This discontinuity in velocity is transmitted thorugh the entire system and is accompanied by an infinitely high force with a finite impulse (time integral). These discontinuities and impulses can be cal- culated, as shown by Hatze and Venter (1981) but the equations are complex. It may also be more realistic to model the contact by a force-deformation model, especially when sim- ulation results are compared to experimental force measurements which do not show impulsive forces. Implementation of a force-deformation model, sometimes referred to as a "soft constraint", requires three components: (1) An equation that calculates the deformation of the element from the generalized coordinates of the system. When forces are velocity-dependent, the deformation veloc- ity is also required, which depends in the generalized velocities . (2) A constitutive equation that calculates the force as a function of deformation and deformation velocity . (3) Equations that convert the external force into generalized forces (or forces and moments) that can be applied to the equations of motion. A disadvantage of force-deformation models is that they may make the system of differen- tial equations stiff, especially when the contact stiffness is high and the mass of contacting body segments is small relative to the total mass of the system (Bogert et al., 1989). Steps (1) and (3) are very straightforward to program when the Newton-Euler formulation is used for the equations of motion and more complicated for the Lagrange formulation. The procedures are similar to the implementation of muscle forces as shown in the previous section. The following constitutive equation has been used for foot-ground contact during running (Gerritsen, 1995): Suitable values for the parameters and can be determined by fitting this equation to experimental force-deformation data. Contact forces are often distributed forces, so it may be necessary to use many such contact elements, discributed over the entire contact sur- face. Care must be taken that the stiffness parameters are scaled to be proportional to the contact area represented by one contact element, as was done for a 3-D model of an ath- letic shoe (Wright et al., in press). The resulting model is a discrete element model (DEM) rather than a finite element model (FEM) since there is no interaction between force and deformation in different elements. A DEM is considerably faster to compute than a FEM and usually adequate for incorporating contact in human movement simulations, provided we are not interested in accurate predictions of stresses in the contact material. Aerodynamic forces are important during movements at high speed. A good examply is the simulation of ski jumping (Müller, 1986), which would be completely unrealistic without aerodynamic forces. In fact, the entire purpose of such simulations is to determine how athletes can make the best use of aerodynamic forces. Practically, this is done by collect- ing wind tunnel data and develop regression models that predict the total aerodynamic force on the body as a function of the generalized coordinates . The force is typically represented as a lift and a drag component, each proportional to the square of the velocity by a factor which is a function of the orientation of the body segments relative to the air- flow. Such equations are adequate to model the effect of aerodynamics on the flight trajec- tory, but it should be noted that this does not provide information about the distribution of these forces over the different body segments. This would require finite element models of air flow, which are extremely computationally expensive and has only recently been applied to this problem (Asai et al., 1997). Another application where aerodynamic forces are important is the simulation of javelin flight in order to find optimal release conditions (Hubbard and Rust, 1984). 4.9.8. OPTIMIZATION STUDIES As defined in section 4.9.1, simulation consists of performing experiments on a numerical model. Since it is easy to perform hundreds, sometimes thousands of experiments on a numerical model, a special type of experimentation, optimization, becomes feasible. The first optimization studies using a musculoskeletal model were done by Hatze (1976) to find the optimal combination of muscle activation patterns in five muscles that produced the fastest kicking movement in a two-segment model of the lower extremity. In a more recent simulation of vertical jumping (Pandy and Zajac, 1990) the goal was to find the optimal muscle activation patterns for a maximal-height jump. Solutions to such problems can be found using iterative methods, which require the execution of many simulations while the unknowns are automatically varied by a suitable search algorithm. Traditionally, optimal control methods have been used which solve for arbitrary functions , but more recently such problems have been formulated as parameter optimizations (Pandy et al., 1992; Soest et al., 1993). The unknown functions are described using equations with unknown parameters, for instance switching times (Soest et al., 1993) or function values sampled at certain time intervals (Pandy et al., 1992). Then the problem is simply to find the parameter vector which results in the best performance. Perfor- mance is a scalar function of , which can be calculated in two steps: (1) perform a sim- ulation with as input, and (2) determine the performance (e.g. jumping height) from the results of the simulation. Performance measures are easily defined for maximal-effort tasks such as mentioned above. In submaximal movements, one might consider finding a strat- egy which minimizes energy consumption, muscle forces, or some other optimization cri- terion. This is an elegant method to investigate which criteria govern control of movement. If the result of the optimization is similar to observed human behaviour, support is obtained for the validity of the hypothesized control principle. A good example is the notion that movements are executed with minimal "jerk" or rate of acceleration (Flash and Hogan, 1985). Although muscle coordination is most frequently optimized, other model parameters can be optimized using the same methods. A good example is the problem of designing an optimal bicycle with respect to a certain performance criterion (Bogert, 1994). When opti- mizing design of equipment used by humans, care must be taken to optimize the human's muscle coordination and equipment at the same time. Humans can and will adapt to changes in equipment and this learning effect may be very significant. Such design prob- lems therefore result in optimizations of a large number of simultaneous unknowns. Performance optimization is not always the objective of optimization. In some cases we merely wish to create a model which replicates observations on human subjects. This requires solution of the so-called tracking problem. The inputs are parameterized as before, but the function now quantifies the difference between observed and simulated variables. In humans, we can typically observe only external forces and kinematics. If is the simulated result for variable , and is the corresponding observation, a suitable weighted least-squares formulation for the objective function is (Neptune and Bogert, 1998): where is the variance between multiple observations of variable , averaged over the duration of the movement. Thus, more reliable observations are weighted more heavily and the final result after optimization will be that tracking errors in each variable are pro- portional to the experimental variability of that variable. One reason for solving the track- ing problem for a musculoskeletal model to obtain information about internal muscle forces during a specific movement. This problem can also be solved using inverse dynam- ics and applying a static optimization to solve the distribution problem, as shown earlier in this chapter, but there are two important differences. First, when optimizing a forward dynamics model it is automatically guaranteed that the solutions are consistent with prop- erties of the system. Specifically, this procedure will only produce muscle forces that do not exceed the known capabilities of the muscle, i.e. maximal force as a function of length and velocity and maximal rate of force development. Secondly, if the number of unknown parameters is small enough, a mathematically unique solution is obtained without requiring additional optimization criteria to solve the distribution problem. Effectively, the distribution problem is avoided by assuming simple activation functions of the mus- cles, since these were parameterized using a small number of parameters. This can be jus- tified by assuming minimal complexity of the control system that generates the signals . Solving the tracking problem requires orders of magnitude more computation time than solving the inverse problem and is just becoming feasible for complex musculoskele- tal models (Neptune and Hull, 1998). Simulations developed with optimization of the tracking criterion, once validated, can be used for experiments that are not possible on human subjects (Gerritsen et al., 1996) Many numerical methods for optimization of nonlinear functions have been developed (Press et al., 1992, chapter 10). All of these start from an initial guess of , and then itera- tively improve the solution by repeatedly calculating while varying the inputs . The most efficient algorithms are based on gradient information, i.e. the derivatives of with respect to each of the parameters . Gradient information speeds up the search because it indicates in which direction each parameter needs to be changed in order to maximize the function However, gradient-based algorithms have not performed well for optimization of complex movements, for several reasons. First, the gradients may be unreliable. If a muscle activation parameters include switching times, gradients can be anywhere between zero and infinity due to the division of time into steps of size when performing the sim- ulations (e.g. Soest et al., 1993). The simulation will then remain unchanged if a switching time is changed by a sufficiently small amount. In other words, the function is piecewise constant. This problem may be overcome by using optimization methods that do not require gradients, such as the simplex method (Press et al., 1992; Matlab function fmins). A more serious problem of most optimization methods is the tendency to converge on local optima rather than the global optimum (Fig. 1). To a certain extent, this can be pre- vented by solving the optimization problem many times with different initial guesses. However, some models have many local optima and the global solution may never be found. The simulated annealing (SA) method (Press et al., 1992) is an elegant and power- ful method for global optimization, based on the analogy of minimization of potential energy in a crystal lattice which is a function of the positional coordinates of all nuclei. At high temperature, large random changes in position occur. Lattice configurations with low energy will be encountered from time to time and the system likely to stay longer in such a state. As the temperature decreases, the system will "freeze" in the state of lowest energy. It has been proved analytically that with a sufficiently slow decrease in temperature, the SA method will converge to the global optimum. Unfortunately, this theoretically derived temperature reduction scheme is to slow for practical applications, and "quenching" (a faster temperature reduction) must be used. Nevertheless, the performance of the SA method for human movement optimization has been shown to be faster as well as find a better solution than other methods (Neptune, in press). As an example, we will determine the optimal release conditions for the throw of a javelin, based on a simplification of the model by Hubbard and Rust (1984). The model is a single rigid body, with three degrees of freedom in the plane of motion. Sideways movements will be neglected. We will assume that the javelin is always released at a height of 1.5 meters and with a velocity of 30 m/s. The athlete can choose the three remaining release parameters: orientation of the javelin, angular velocity of the javelin, and direction of the throw (Fig. 2). The Matlab program in listing 1 simulates a single throw. It is written as a Matlab function which takes the three release parameters as input and returns the negative of the distance as output. Type javelin([0.5 0.0 0.1]) to see the result of a throw with initial orientation of 0.5 radians, zero angular velocity and an angle of attack of 0.1 radians. Now type fmins(`javelin',[0.5 0.0 0.1]) which performs an optimization of the distance using the simplex method, with those same release parameters as an initial guess. For this model, the folliwing optimal release conditions were found: = 33.4°, = 0.16°/s, and = -3.7°. 4.9.9. SIMULATION AS A SCIENTIFIC TOOL When simulation is used as a scientific tool, experiments are performed on a numerical model. Direct validation of a numerical model is usually difficult because the same exper- iments can not be done on human subjects. This could be because the experiments are related to severe injury, because human subjects are not sufficiently reproducible, because humans are too fatiguable, or because the outcome is a variable which can not be mea- sured. How then do we ensure the validity of scientific studies using computer models? Only indirect methods are available. For this reason, the term "validation" may be too strong and "evaluation" should be used instead. First of all, a model should be consistent with observations that can be made on humans. When optimization of performance is carried out, and a realistic movement is obtained, the model is generally considered to be valid because no movement data were used to develop the simulation. When solving the tracking problem (EQUATION), it is expected that after optimization all variables are within two standard deviations of the mean. If this is not the case, the model is unable to perform the movement in a sufficiently realistic manner. This should be reason to closely examine the model and make improvements where necessary. Passing this test, however, does not guarantee a valid model. Due to the reduncancy of the locomotor system, the model could have found a different solution than the human in order to achieve the same external movement and force variables. In that case, additional predictions must be elicited from the model which can then be compared to measurements that were not used for development of the model. When the consistency test is not sufficient, as is generally the case when the tracking prob- lem is solved, one must test the response of the model to controlled interventions and com- pare that response to results of the same experiments on humans. Even when the final application is a study on severe injuries, is is often still possible to evaluate the model dynamics using non-desctructive perturbation tests. Care must be taken that these experi- ments test those aspects of the model that are important for the final application. For devel- opment of a model for ankle sprains, small perturbations of surface orientation were applied to model and humans to evaluate the validity of the model (Wright et al., submit- ted). Finally, a model should not be overly sensitive to errors in model parameters. Critical model parameters can be identified by sensitivity analysis: each parameter is adjusted by a small amount and the change in the results of the simulation is examined. In some cases, this will show that certain model parameters are too critical and the results of a simulation study would depend entirely on a random error in such a parameter. In certain cases, opti- mization methods are helpful. In a quasi-static model of the knee joint, it was found that the behaviour was sensitive to the lengths of the ligaments (Blankevoort and Huiskes, 1996). These unknown parameters were then eliminated by solving the tracking problem. Ligament lengths which minimized the difference between simulated and measured move- ments were found and these values were used for subsequent applications of the model. Another powerful safeguard against overly sensitive models is statistical analysis. Experi- ments on humans or animals are never performed on a single individual, because one indi- vidual may not be representative of the population. Statistical analysis is performed to ensure valid generalizations. When using complex musculoskeletal models, the same prin- ciple should apply. These models have many degrees of freedom, many natural frequen- cies, and are often unstable and chaotic. Results from a single model could well be completely irrelevant. In a simulation study on the effect of shoe hardness on impact forces in running, both positive and negative responses were found in a group of models (Wright et al., in press). By examining the model, this could be attributed to two mechanisms which worked in opposite directions. Impact forces tend to increase initially with harder material. This then increases the rate of knee flexion, resulting in a better shock absorption by the body. In subjects with a certain movement style, the latter mechanism resulted in an overcompensation. Statistical analysis showed no significant effect for the population as a whole, confirming earlier results on human subjects. Thus, statistical analysis prevented incorrect generalizations which could have been made when just one model had been used. In this case, the model was sensitive to movement style. Sensitivities to other variations withn the human population, such as anatomical variations, may be eliminated similarly, by creating a population of models with the appropriate range of parameter values. Note, however, that modeling and simulation are sometimes used specifically to determine the influence of inter-subject variations, and in such case the statistical approach is not appro- priate. Summarizing: it should be kept in mind that simulation experiments only tell the truth about the model that was used. Generalization to the human population is always hazard- ous and requires extensive validation and careful examination of the results. 4.0 ADDITIONAL EXAMPLES This section has been divided into basic and advanced examples. In general, the basic questions can be answered by studying the text or by performing some calculations based on information provided by the text. The advanced questions go beyond what can be found in the text. BASIC QUESTIONS: 1. In the simulation of the downhill coasting cyclist (Listing 1), determine which value of the step size is required to get at most a 5% error in the speed of the cyclist at any given time. 2. Modify the simulation program in Listing 1 to add the propulsive force due to pedal- ing. Consider different ways to model this: (a) constant force, and (b) constant power. 3. Modify the program in Listing 1 to compute , the distance traveled by the cyclist, by including the ODE in the model. Then include the effect of change of altitude on the air drag coefficient . 4. Modify the equations of motion for arm movement, and the simulation model (Listing 2) to include constant joint torques at both joints. 5. Modify the equations of motion for arm movement, and the simulation model (Listing 2) to simulate movement in a vertical plane by including gravitational forces. Simulate this two-segment pendulum. What is different, compared to the simulated movement in the horizontal plane? 6. Modify the program in Listing 3 to include the eccentric part of the force-velocity rela- tionship (i.e., ). Develop a suitable mathematical model from the following information: (a) The force-velocity relationship must be continuous at zero velocity, and (b) At high eccentric velocities, the muscle force approaches 1.5 asymptotically. Use the modified program to simulate isokinetic lengthening protocols. 7. How many state variables are there in the model of the ballistic knee extension (List- ing 4)? 8. Imagine an unloaded knee extension task with the object of hitting a target with maxi- mal foot velocity. Assume that the target is always reached at full knee extension, and use a series of simulations to determine the optimal quadriceps moment arm for this task. How does this optimal moment arm compare to actual moment arms reported in the literature (Visser et al., 1991)? 9. Modify the simulation of knee extension to include a model for activation dynamics. Assume an initial condition =0 at =0. How important is activation dynamics for this movement? REFERENCES Andrews, J.G. (1985) A General Method for Determining the Functional Role of a Muscle. J. Biomech. Eng. 107, pp. 348-353. Arnold, V.I. (1978) Mathematical Methods of Classical Mechanics. Springer Verlag, New York, Berlin, Heidelberg. Asai, T., Kaga, M., and Akatsuka T. (1997) Computer Simulation of the V-style Technique in Ski Jumping using CFD. Proc. 6th Int. Symp. Computer Simulation in Biomechanics, Tokyo, Japan, pp. 48-49. Bogert, A.J. van den, Schamhardt H.C., and Crowe, A. (1989) Simulation of Quadrupedal Locomotion Using a Dynamic Rigid Body Model. J. Biomechanics 22, pp. 33-41. Bogert, A.J. van den (1994) Optimization of the Human Engine: Application to Sprint Cycling. Proc. 8th Congress of the Canadian Society for Biomechanics, Calgary, pp. 160-161. Bogert, A.J. van den, Gerritsen, K.G.M., and Cole, G.K. (1998) Human Muscle Modelling from a User's Perspective. J. Electromyogr. Kinesiol. 8, pp. 119-124. Blankevoort L. and Huiskes, R. (1996) Validation of a Three-Dimensional Model of the Knee. J Biomechan- ics 29, pp. 955-961. Cole, G.K., Bogert, A.J. van den. Herzog, W., and Gerritsen, K.G.M. (1996) Modelling of Force Production in Skeletal Muscle Undergoing Stretch. J Biomechanics 29, pp. 1091-1104. Delp, S.J. and Loan, J.P. (1995) A Graphics-Based Software System to Develop and Analyze Models of Musculoskeletal Structures. Comput. Biol. Med. 25, pp. 21-34. Flash T. and Hogan, N. (1985) The Coordination of Arm Movements: An Experimentally Confirmed Mathe- matical Model. J. Neurosci. 5, pp. 1688-1703. Gear, C.W. (1971) Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, Englewood Cliffs, NJ. Gerritsen, K.G.M., Bogert, A.J. van den, and Nigg, B.M. (1995) Direct Dynamics Simulation of the Impact Phase in Heel-Toe Running," J. Biomechanics 28, pp. 661-668. Gerritsen, K.G.M., Nachbauer, W., and Bogert, A.J. van den (1996) Computer Simulation of Landing Move- ment in Downhill Skiing: Anterior Cruciate Ligament Injuries," J Biomechanics 29, pp. 845-854. Gerritsen, K.G.M., Bogert, A.J. van den, and Herzog, W. (in press) Force-Length Properties of Lower Extremity Muscles Derived from Maximum Isometric Strength Tests. Eur. J. Appl. Physiol. Grieve, D.W., Pheasant, S., and Cavanagh P.R. (1978) Prediction of Gastrocnemius Length from Knee and Ankle Joint Posture. In: Asmussen E., Jorgensen K. (eds.) Biomechanics VI-A. University Park Press, Baltimore, pp. 405-412. Groot, G.de, (1995) Air Friction and Rolling Resistance during Cycling. Med. Sci. Sports Exerc. 27, pp. 1090-1095. Hatze, H. (1976) The Complete Optimization of a Human Motion. Math. Biosci. 28, pp. 99-135. Hatze, H. (1981) A Comprehensive Model for Human Motion Simulation and its Application to the Take-off Phase of the Long Jump. J. Biomechanics. 14 (3), pp. 135-142. Hatze, H. and Venter, A. (1981) Practical Activation and Retention of Locomotion Constraints in Neuromus- culoskeletal Control System Models. J. Biomechanics 14, pp. 873-877. Haug, E.J. (1989) Computer-aided Kinematics and Dynamics of Mechanical Systems. Basic Methods. Allyn & Bacon, Boston. He, J., Levine, W.S., and Loeb, G.E. (1991) Feedback Gains for Correcting Small Perturbations to Standing Posture. IEEE Trans. Autom. Cont. AC-3, pp. 322-332. Hubbard, M. (1980) Dynamics of the Pole Vault. J. Biomechanics 13, pp. 965-976. Hubbard M. and Rust H.J. (1984) Simulation of Javelin Flight using Experimental Aerodynamic Data. J Bio- mechanics 17, pp. 769-776. Ingen Schenau, G.J. van and Bobbert, M.F. (1993) The Global Design of the Hindlimb in Quadrupeds. Acta Anat (Basel) 146, pp. 103-108. Kane, T.R. and Levinson, D.A. (1985) Dynamics : Theory and Applications. New York, McGraw-Hill. Müller, W., Platzer, D., and Schmolzer, B. (1996) Dynamics of Human Flight on Skis: Improvements in Safety and Fairness in Ski Jumping. J. Biomechanics 29, pp. 1061-1068. Neptune, R.R. and Bogert, A.J. van den (1998) Evaluation of Mechanical Energy Analyses: Application to Steady-state Cycling. J. Biomechanics 31, pp. 239-245. Neptune, R.R. and Hull, M.L. (1998). Evaluation of Performance Criteria for Simulation of Submaximal Steady-state Cycling using a Forward Dynamic Model. J. Biomech. Eng. 120(3), pp. 334-341. Neptune, R.R. (in press) Optimization Algorithm Performance in Determining Optimal Controls in Human Movement Analyses. J. Biomech. Eng. Pandy, M.G. and Zajac, F.E. (1990) Optimal Muscular Coordination Strategies for Jumping. J Biomechanics 24, pp. 1-10 Pandy, M.G., Anderson, F.C.,and Hull, D.G. (1992) A Parameter Optimization Approach for the Optimal Control of Large-scale Musculoskeletal Systems. J Biomech Eng 114, pp. 450-460. Passerello, C.E. and Huston, R.L. (1971) Human Attitude Control. J. Biomechanics. 4 (2), pp. 95-102. Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P. (1992) Numerical Recipes in FORTRAN, 2nd Ed., Cambridge University Press, Cambridge, UK. Shampine, L.F. and Gordon, M.K. (1975) Computer Solution of Ordinary Differential Equations: The Initial Value Problem. W.H. Freeman, San Francisco, CA. Soest, A.J. van, Schwab, A.L., Bobbert, M.F., and Ingen Schenau, G.J. van (1993) The Influence of the Biar- ticularity of the Gastrocnemius Muscle on Vertical-Jumping Achievement. J Biomechanics 26, pp. 1-8. Soest, A.J.van and Bobbert, M.F. (1993) The Contribution of Muscle Properties in the Control of Explosive Movements. Biol. Cybern. 69, pp. 195-204. Spoor, C.W., Leeuwen, J.L. van, Meskers, C.G.M., Titulaer, A.F., and Huson, A. (1990) Estimation of Instantaneous Moment Arms of Lower-leg Muscles. J. Biomechanics. 23 (12), pp. 1247-1259. Visser, J.J., Hoogkamer, J.E., Bobbert, M.F., and Huijing, P.A. (1990) Length and Moment Arm of Human Leg Muscles as a Function of Knee and Hip-joint Angles. Eur. J. Appl. Physiol. 61 (5-6), pp. 453- 460. Walker, S.M. and Schrodt, G.R. (1973) I-Segment Lengths and Thin Filement Periods in Skeletal Muscle Fibres of the Rhesus Monkey and the Human. Anat. Rec. 178, pp. 63-82. Winter, D.A. (1979) Biomechanics of Human Movement. Appendix A. John Wiley & Sons, New York. Winters, J.M., Stark, L. (1987) Muscle Models: What is Gained and What is Lost by Varying Model Com- plexity. Biol. Cybern. 55, pp. 403-420. Winters, J.M. (1990) Hill-Based Muscle Models: A Systems Engineering Perspective. In: J.M. Winters & S.L.-Y. Woo (eds.) Multiple Muscle Systems: Biomechanics and Movement Organization. Springer-Verlag, New York. Wright, I.C., Neptune, R.R., Bogert, A.J. van den, and Nigg, B.M. (in press) Passive Regulation of Impact Forces in Heel-Toe Running. ". Clin. Biomech. Wright, I.C., Neptune, R.R., Bogert, A.J. van den, and Nigg, B.M. ( submitted) Validation of a 3D Model for the Simulation of Ankle Sprains". J. Biomech. Eng. Yamaguchi, G.T. and Zajac, F.E. (1990) Restoring Unassisted Natural Gait to Paraplegics via Functional Neuromuscular Stimulation: A Computer Simulation Study. IEEE Trans. Biomed. Eng. BME-37, pp. 886-902. Yeadon, M.R. (1993c) The Biomechanics of Twisting Somersaults, Part 3: Aerial Twist. J. Sports Sciences. 11 (3), pp. 209-218. Yoon, Y.S. and Mansour, J.M. (1982) The Passive Elastic Moment at the Hip. J. Biomechanics 15, pp. 905- 910. Zahalak, G.I. (1986) A Comparison of the Mechanical Behavior of the Cat Soleus Muscle with a Distribu- tion-Moment Model. J. Biomech. Eng. 108, pp. 131-140. Zajac, F.E. (1989) Muscle and Tendon: Properties, Models, Scaling, and Application to Biomechanics and Motor Control. Crit. Rev. Biomed. Engng. 17, pp. 359-411. Zajac, F.E. and Gordon, M.E. (1989) Determining Muscle Force and Action in Multi-articular Movement. Exercise and Sport Science Reviews (ed. Pandolf, K.). Williams & Wilkins, Baltimore, MD. 17, pp. 187-230. Figure NEW1: The three-component Hill muscle model. CE: contractile ele- ment, PEE: parallel elastic element, SEE: series elastic element. Fig. 1: Hypothetical performance function f as a function of a single parameter p. Optimization of f from the initial guess p1 may result in finding the local maximum m2, while the more important global maximum is located at m1. Fig. 2: Definition of release conditions of the javelin: orientation , angular veloc- ity , and angle of attack . In this example, the angle of attack is negative, indi- cating that the javelin is oriented clockwise relative to the velocity vector . Figure 4.9.1 Free body diagram of a cyclist coasting downhill. Figure 4.9.2 Simulation of a cyclist coasting downhill. Solid curve: =0.15 Ns/m, dashed curve: =0.135 Ns/m. Figure 4.9.3 Effect of step size on the simulation program in Listing 1. Solid line: =0.01 s, dashed line: =2 s, dotted line: =10 s, dash-dotted line: =20 s. Figure 4.9.4 Arm movement simulated using program in Listing 2. Solid line: segment 1, dashed line: segment 2. Figure 4.9.5 Simulation of a muscle undergoing isokinetic ramp shortening between t=1 and t=2 s. Left: imposed length change, right: force generated by the muscle. See text and listing 3 for details of the model. Figure 4.9.6 Simulation of a maximal effort ballistic knee extension (Listing 4). Solid line: with muscle moment arm =33 mm. Dashed line: =66 mm.