KAAP686
Neuromusculoskeletal Modeling
W. Rose
2017-02-22
Notes on Buchanan & Manal (2004), Neuromusculoskeletal Modeling, J. Appl. Biomech.30: 367-395.
EMG-based Forward Dynamics vs. Inverse Dynamics
EMG-based Forward Dynamics
Predict forces & movements from the EMGs. Requires a model of relation between EMG & muscle force, and a model of relation between muscle force & movement. Each of these two stages has many uncertainties. First step (EMG to force) is based on the idea that both EMG and force are results of the neural commands to muscle. It is possible or likely that the transfer function from neural command to EMG varies with time, position, coactivation, etc., and that the transfer function from neural command to force also varies with those factors, but that they don’t vary the exact same way, and therefore the transfer function between EMG and force is state-dependent, i.e. not constant. The second step (force to movement) requires detailed knowledge of muscle lines of action, moment arms, moments of inertia, etc. Furthermore, because position is doubly-integrated version of force (or torque), small error in force can accumulate to make large error in position.See integration_error.xls.
Inverse Dynamics
Measure displacements and external forces, estimate muscle forces. Done by differentiating positions and angles twice to get accelerations and angular accelerations. From those, and from external forces and from limb inertias and moments of inertia, estimate forces and torques across joints. From those forces and torques one might try to estimate muscle forces and neuralcommands.
Problems with this approach: Accuracy of the estimated forces is limited by accuracy of the estimates of limb inertias & moments of inertia. Differentiation and double differentiation amplifies noise. Joint reaction forces and torques are net values. Therefore any antagonist or agonist co-contraction effects cannot be determined or teased apart. Finally, there are not good models to obtain muscle activation from muscle force. Because force tends to be a low-passed version of activation, it is easier to go from activation to force than vice versa.
Muscle Activation Dynamics
Chart above is from muscle_activation_dynamics.xlsx.
How to convert EMG to muscle force? First, convert raw EMG to smoothed EMG by rectify and low-pass filter (“linear envelope” filter). Then normalize by EMGmax, measured during MVIC. This should give e(t) in range (0, +1). This signal is converted to muscle activation a(t) by a nonlinear dynamic process. We do this by using a pure delay cascaded with a linear dynamic system cascaded with a static, or instantaneous (means the same thing) nonlinearity. The delay converts e(t) to e(t-d), where d is magnitude of the delay. The linear dynamic system converts e(t) to u(t). (Like e(t), u(t) will be in the range (0,1).) The linear system has unity gain at DC.[1] Then an “instantaneous” (meaning it has no dynamics, no smoothing, etc) nonlinearity is used to convert u(t) to a(t), where a(t) is also in the range (0, +1). The nonlinearity is introduced because if one plots steady isotonic force vs. EMG for various levels of force, one finds a nonlinear relation, especially near the low end. (See Fig. 4A in paper.) This results in the muscle activation signal a(t).
Converting Muscle Activation to Force
Muscle activation is converted to muscle force by scaling, with some more complicated effects added as well. First the simple part: active muscle force = FA(t) = Fmax*a(t), where Fmax= peak muscle force. We cannot stop here for two reasons: force-velocity relation and force-length relation. In general, a muscle that is shortening produces less force than an isometric muscle. We need to take this into account. Furthermore, a muscle has a length, or range of lengths, at which it produces peak active force, and at shorter or longer lengths it produces less force.
Effects of length
A muscle generates some force due to its passive elastic properties. The active and passive components of the muscle are considered to be in parallel, which means the forces generated by them are additive. Thus muscle force is now
F(t,L) = FA(t) + FP(L)
where FA(t)=Fmax*a(t)*fA(L) and FP(L)=Fmax*fP(L). fA(L) varies from 0 to 1 and and fP(L) varies from 0 to infinity. Both are actually functions of normalized muscle length, i.e. length divided by “optimum length”, where optimum length = L0 = length at which active force is maximal. fA(L) captures the length-dependency of muscle force. This length dependence appears to be instantaneous in the model used here. It can be understood as a consequence of the number of overlapping crossbridges among the sliding actin and myosin filaments, although the real story is more complicated. The optimum length is somewhat activation-dependent (this is an experimental finding): L0 gets somewhat shorter (about 15% or so in some muscles) as activation a(t) increases from 0 to 1. A muscle that is <50% or > 150% of optimum length generates no force, i.e. fA(L/L0) is 0 if L/L0 <0.5 or if L/L0 > 1.5.
Effects of shortening or lengthening velocity
The Hill equation is used to adjust FA(t,L):
FA(t,L) = Fmax*a(t)*fA(L)
is replaced with
FA(t,L) = [(Fmax*b - a*v)/(v+b)]*a(t)*fA(L)
where a and b are constants that describe the shape of the force-velocity relationship.
Modeling the tendon
A smooth and slightly nonlinear relation is used. Force in muscle and force in tendon are considered to be always and instantaneously equal since they are in series.
Pennation angle
Component of muscle force along long axis, in series with the tendon is
Flong(t,L) = F(t,L)*cos(φ(t))where φ=pennation angle. φ=0 means muscle fibers are along the long axis of muscle. Pennation angle changes as muscle length changes. The effect is considered to be instantaneous and is modeled by eq. 33 which simplifies to
L(t)*sin(φ(t)) = L0*sin(φ0), where φ0 = pennation angle at optimum muscle fiber length L0. This can be rewritten
sin(φ(t))/sin(φ0) = L0/L(t). (Reminiscent of Snell’s law for refraction.) It keeps the width of the muscle constant, since the width is L(t)*sin(φ(t)).
Adjusting Parameters
Adjustable parameters that will not change with time include maximum force Fmax, optimum fiber length L0, pennation angle at optimum length φ0, and tendon slack length. These are determined separately for each muscle. Maximal shortening velocity vmax (related to Hill force-velocity eqn) is considered constant, from the literature, same in all muscles modeled.
Running the Model
A second order linear differential equation relates e(t) and u(t) for each muscle. The many other nonlinearities present are also accounted for to get a nonlinear system of differential equations. This can be written as a system of first order (ordinary, not partial) differential equations. The solution can be obtained with a Runge-Kutta algorithm with adaptive stepsize. Runge-Kutta-Fehlberg is one such algorithm, Matlab’s ode45() is a very similar Runge-Kutta routine. Syntax for using ode45 is
[T,Y] = ode45(odefun,[Tstart Tstop],y0)
where ode45() puts values into the vector T and the array Y. T contains the times at which solution is calculated; the columns of Y are the different variables Y1(t), Y2(t), etc. odefun is a handle with the name of a function defining the first derivatives of the Yi’s; the vector y0 has the initial conditions for the Yi’s.
[1] The linear dynamic system used is a second order system, analogous to a mass-spring-dashpot system. It can be modeled like other second order lowpass filters. The particular form chosen by the authors is a second order model of the form
u(t) + b1*u(t-1) + b2*u(t-2) = a*e(t), i.e. u(t)=a*e(t) - b1*u(t-1) - b2*u(t-2).
The requirement of unity gain means that a=1+b1+b2. Therefore there are only 2 adjustable parameters: b1 and b2. We can convert b1, b2 to an equivalent pair of variables g1 and g2 which are equally good at characterizing the dynamic system. The advantage of using g1 and g2 is that the allowed values for g1 and g2 that guarantee filter stability (i.e. non-infinite values) are simple: |g1|<1 and |g2|<1. The relation between b1,b2 and g1,g2 is that 1+b1*x + b2*x2 = (1 + g1*x)(1+g2*x), i.e. b1=g1+g2 and b2=g1*g2.
A graph of the allowed values of g1,g2 and the corresponding allowed values of b1,b2 is in spreadsheet muscle_activation_dynamics.xlsx. The same spreadsheet shows impulse response of the system as g1, g2 are varied. It turns out that if either g1 or g2 or both are >0, the system’s response to an impulse (i.e. the force of a twitch) is oscillatory, which is non-physiologic. If either g1 or g2=0 the response is a simple decaying exponential which is also non-physiologic. If both g1 and g2 are <0 (but >-1 for stability) then physiological twitch forces are generated. See muscle_activation_dynamics.xlsx.