A joint-space numerical model of metabolic energy expenditure for human multibody dynamic system.

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. (2015); e02721 Published online 24 June 2015 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.2721

A joint-space numerical model of metabolic energy expenditure for human multibody dynamic system Joo H. Kim*,† and Dustyn Roberts Department of Mechanical and Aerospace Engineering, New York University, Brooklyn, NY, USA

SUMMARY Metabolic energy expenditure (MEE) is a critical performance measure of human motion. In this study, a general joint-space numerical model of MEE is derived by integrating the laws of thermodynamics and principles of multibody system dynamics, which can evaluate MEE without the limitations inherent in experimental measurements (phase delays, steady state and task restrictions, and limited range of motion) or muscle-space models (complexities and indeterminacies from excessive DOFs, contacts and wrapping interactions, and reliance on in vitro parameters). Muscle energetic components are mapped to the joint space, in which the MEE model is formulated. A constrained multi-objective optimization algorithm is established to estimate the model parameters from experimental walking data also used for initial validation. The joint-space parameters estimated directly from active subjects provide reliable MEE estimates with a mean absolute error of 3.6 ± 3.6% relative to validation values, which can be used to evaluate MEE for complex non-periodic tasks that may not be experimentally verifiable. This model also enables real-time calculations of instantaneous MEE rate as a function of time for transient evaluations. Although experimental measurements may not be completely replaced by model evaluations, predicted quantities can be used as strong complements to increase reliability of the results and yield unique insights for various applications. Copyright © 2015 John Wiley & Sons, Ltd. Received 18 November 2014; Revised 23 March 2015; Accepted 23 April 2015 KEY WORDS: first law of thermodynamics; generalized coordinates; heat dissipation; joint space;

mechanical work; metabolic energy expenditure; muscle activation

1. INTRODUCTION Metabolic energy expenditure (MEE) is an important performance measure of human motion [1]. Evaluations of MEE are commonly used in the analysis of human movement to provide insight in a wide range of scientific, clinical, and engineering fields, including gait analysis [2–4], sports medicine and exercise prescription [5, 6], analysis and treatment of pathological movement [3], design of lower limb prostheses and orthoses [7, 8], geriatrics [9], and obesity [10]. MEE is also a critical indicator of fatigue [11]. Metabolic energy expenditure rate is usually evaluated through indirect calorimetry estimates of _ 2 ) [5, 6, 12, 13] or through formulas and look-up tables constructed from oxygen uptake rate ( VO experimental data of common tasks [14]. However, the predictive capabilities of such experimentbased methods are limited by specific conditions (e.g., steady state and aerobic metabolic range) and the scope of the data. The measurement equipment also limits the tasks and ranges of motion appropriate for testing and does not provide comprehensive real-time MEE data because of the phase _ 2 reaches steady state [6]. For these delay between the actual onset of energy expended and when the VO *Correspondence to: Joo H. Kim, Department of Mechanical and Aerospace Engineering, New York University, Brooklyn, NY, USA. † E-mail: [email protected] Copyright © 2015 John Wiley & Sons, Ltd.

(1 of 22) e02721

e02721 (2 of 22)

J. H. KIM AND D. ROBERTS

reasons, only the energetics of conventional tasks such as steady-state walking [3, 4, 7, 9, 15] and cycling [16] have been frequently studied, and even in those tasks, the results are not consistent (e.g., disagreements involving the relative costs of stance and swing phases during walking [2]). The limitations of experimental evaluations are also the causes of the conflicted treatment of efficiency in literature, in which the components of the efficiencies and the resulting values differ between studies [17–22]. The constant efficiencies used in these approaches do not represent actual dynamic efficiencies or address different efficiencies between muscles [23]. Furthermore, the efficiencies of uphill and downhill walking [18] are commonly adopted as those of positive and negative mechanical work, respectively [7, 24], although any walking motion includes both positive and negative work phases concurrently [25]. In general, even a simple task is not performed in an identical way by a given subject on different trials. Because different joint trajectories and muscle forces will cause differences in MEE, experimentally derived formulas may produce significant differences from actual measurements in any given task [26]. _ 2 Models of MEE based on individual muscle experiments can be used without measuring VO and encountering the aforementioned limitations on experiments. Additionally, instantaneous MEE rates from models allow integration over periods of interest to estimate the relative costs of different phases of a task [1, 2]. However, multiple musculotendon units that cross a given joint can be activated as synergists and can work in combination to create a net vector resultant of moments of muscle forces about an instantaneous joint center [27]. Because there are approximately 600 functional skeletal muscles in the human body, including redundant agonists and antagonists at most joints, the number of DOFs in the muscle space is too large for deterministic and verifiable evaluation [28]. Furthermore, because of the compliance and flexibility of soft tissues, the muscle space is essentially an infinite dimensional vector space that includes the geometric and mechanical configurations of individual muscles. Several time-varying interactions and properties of muscles, such as normal and frictional contacts, wrapping and sliding with internal environment (e.g., other muscles or bones, where muscles may extend over two to three adjacent joints) [29, 30], pennation angles between muscle fibers and tendons [31], moment arms [29, 32], and restoring forces from passive and active elastic tissues [31, 33, 34], cause complexity and indeterminacy in muscle models [28, 33] by affecting the transmission of muscle forces to the resultant muscle-induced joint actuator torque. These interactions become more complicated when the joint angles deviate far from neutral positions [29], which can make the common modeling assumption of straight-line muscle geometry invalid. In addition to the lack of a scalable, quantitative, molecular-level model of muscle energetics in the current literature [35–37], these complexities indicate that the calculation of MEE directly from muscle forces does not consistently provide valid results and is not always possible. Individual muscle models often use muscle parameters scaled from in vitro studies that may not accurately represent those of in vivo muscles. Therefore, existing energy models based on individual muscles are limited in their validity and applications to certain conditions and tasks, and include only a fraction of the muscles in a given system. Even for normal walking, existing muscle models frequently overestimate MEE by up to 50% [2, 38, 39], and the reason for this discrepancy is not well understood [39, 40]. An alternative approach is to map muscle properties and functions to joint space. The resultant of multiple muscles that contribute to a single anatomical joint movement can be represented by a combination of one or more kinematically equivalent revolute joints as generalized coordinates [29, 41], which constitute the joint space. Relative to muscle space, joint space is finite dimensional with much fewer DOFs yet uniquely determines whole-body segment configurations (assuming a rigid-body model of each segment). Because of its computational efficiency, joint space is commonly used in inverse dynamics for motion analysis [21, 29, 42–45] as well as numerical simulation and optimization [46, 47]. Consequently, an MEE model in joint space could be effectively applied to experimentally measured or simulated movement data [42, 48]. Some studies have confirmed the validity of this approach by demonstrating that joint-space kinematics and kinetics alone can predict MEE during walking [49, 50]. By extending these experimentally based findings with a numerical model, MEE can be estimated with commonly measured kinematic and kinetic parameters. It will also enable the calculation and analysis of MEE at each DOF. An Copyright © 2015 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Biomed. Engng. (2015); e02721 DOI: 10.1002/cnm

JOINT-SPACE METABOLIC ENERGY EXPENDITURE MODEL

e02721 (3 of 22)

additional benefit is that the torque–angle–velocity properties in joint space [51] are generally monotonic and less complicated than the force–length–velocity properties in muscle space [31, 33]. In this study, a general joint-space MEE model is systematically derived by integrating the laws of thermodynamics and the principles of multibody system dynamics, as a top-down (deductive) approach, as opposed to the more common empirically based bottom-up (inductive) approach. The first law of thermodynamics has been adopted in a few inductive muscle energy models [17, 40, 52–55]. However, unlike in existing models, our comprehensive derivation includes an explicit term of MEE, and the relevant work and heat terms are rigorously formulated with respect to system boundaries, reference frames, relative/absolute quantities, and external/internal forces, which are not systematically addressed in the current literature. This study introduces the mapping of energetic components (activation, work, and heat) and properties from muscle space to joint space, in which the MEE model is formulated. A constrained multi-objective optimization algorithm is established and used to illustrate the method of model parameters estimation from experimental walking data that is also used for proof-of-concept validation. 2. MODEL DERIVATIONS AND METHODS The MEE model as a function of joint-space variables is derived using the first law of thermodynamics for rigorous analysis and complete description of energy transformation (through various system energy components) and transfer (through work interaction and heat transfer). The first law of thermodynamics for a closed system [56] is _ ext þ Q_ ext E_ ¼ W

(1)

_ ext is the net work rate carried out by nonwhere E_ is the rate of change in total system energy, W _ ext < 0 if _ ext > 0 and W conservative external forces and moments that cross the system boundary (W ext work is carried out on and by the system, respectively), and Q_ is the net heat transfer rate across ext ext the system boundary (Q_ > 0 and Q_ < 0 if heat transfer is to and from the system, respectively). Although this statement of the first law is trivial, identifying and formulating the terms are not, because the first law can be applied in various ways by selectively including components that are relevant to a given problem. Therefore, it is critical in model derivations to identify all relevant work, heat, and energy components in the first law and formulate them with a consistent approach. The principles of multibody system dynamics and the formulation of activation mapping are used to refine each term. The unknown parameters in the model are estimated from experimental walking data using a constrained multi-objective optimization algorithm.

2.1. Identification of energy components and MEE model derivation for the human multibody dynamic system In the proposed model, the entire human body (not only muscles) is identified as the system of interest in order to explicitly include the inherent metabolic energy level as one of the energy components. This is distinct from previous approaches [17, 52–54], in which the metabolic energy term does not appear explicitly in the system energy derivations, which is partially due to the lack of a validated, comprehensive molecular-level model of energy transfer between muscle, and biochemical energy in a useful form [35, 37]. (The principles of chemical energy transformation into directed molecular movement proposed in the existing studies are mostly hypothetical and are not experimentally validated [37].) In this perspective, the proposed deductive approach from the other end can be regarded as a viable alternative. In this study, the human body as the system of interest is constrained to be a closed system [19], so there is no time rate of change of system energy because of mass transfer (e.g., food intake and excrement) through the system boundary. When applied to the segmental motion of the actuated human multibody dynamic system, the k p relevant energy components include the kinetic ( E_ ), potential ( E_ ; due to external conservative p t forces), internal potential ( U_ ; strain due to internal conservative forces), internal thermal ( U_ ) energy rates, and the metabolic energy (as a form of internal biochemical energy) rate. Here, the Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (4 of 22)


rates of work carried out by external and internal conservative forces are included in the system met energy as potential energy terms. Denoting the MEE rate as E_ , the rate of change of the metabolic met energy level within the system is E_ , resulting in the following motion-relevant energy components: k p p t met _ ext þ Q_ ext ¼W E_ þ E_ þ U_ þ U_ E_

(2)

In this equation, the effects of heat on segment volumes (thus the pressure–volume work) and contact properties (friction and elasticity) of the system elements are assumed to be negligible. k When combined with the work-energy principle [57] of a multibody dynamic system, where E_ þ p p _ int þ W _ ext and W _ int is the rate of work carried out by non-conservative internal forces E_ þ U_ ¼ W and moments, the MEE rate is met _ int þ U_ t Q_ ext ¼W E_

(3)

ext t represents the heat dissipation in the form of increased internal The combined term U_ Q_ thermal energy, heat transfer from the human body, or both. Despite their inclusion of similar terms, this formulation is clearly distinct from the first law. In particular, (a) The left-hand side in this formulation is confined to the MEE, while in the first law, various system energy components (including the metabolic energy rate) should be included explicitly. (b) Unlike the first law that is written in terms of external work, the MEE in this formulation is written explicitly in terms of internal work. (c) The comprehensive energy components in the derivation allow inherent inclusion, without the p need for explicit term as in the first law, of the contribution of soft tissue work (U_ ) in the internal work term, which is known to be non-negligible in some human tasks (e.g., walking [58]). _ int and U_ t Q_ ext result from two In this formulation, apart from dissipative components, both W sources—skeletal muscle activation (mus) and basal metabolism (basal) that include the activation of smooth and cardiac muscles to energize vital organs that happens at rest without visible body movement: _ mus þ W _ basal ; U_ t Q_ ext ¼ Q_ mus þ Q_ basal _ int ¼ W W

(4)

Then the heat and work rate terms due to basal metabolism can be combined to constitute the basal metabolic rate (BMR): _ BMR ¼ W

basal

þ Q_

basal

(5)

Finally, the total MEE rate can be written as met _ mus þ Q_ mus þ BMR ¼W E_

(6)

_ mus is the product of skeletal muscle force and shortening/lengthening It should be noted that W velocity as relative quantities, which should be distinguished from the rates of center-of-mass mus is due to work or total work that are commonly used in the literature (e.g., [58, 59]). The Q_ skeletal muscle actuation and is a function of muscle force and shortening/lengthening velocity [49, 60, 61]. Therefore, unlike the first law, which is stated in terms of external work and heat, this formulation of the MEE rate depends only on the relative internal quantities and is not dependent on a specific reference frame. This allows the consistent use of generalized coordinates (joint angles) and the associated quantities (velocities and torques) in each term and provides all the benefits of the joint-space dynamics as discussed previously. Copyright © 2015 John Wiley & Sons, Ltd.



e02721 (5 of 22)

2.2. Activation mapping from muscle space to joint space The Hill-type musculotendon model [31, 62] (Figure 1) includes the force-producing contractile element (CE), the series elastic (SE) and parallel elastic (PE) elements, and the tendon. Generally, the magnitude of the force FM(t) generated by jth muscle (j = 1, …, m where m is the total number of muscles) at time t is related to its maximum isometric force Fmax, muscle activation level aM(t) (∈ [0, 1] for ∀ t), and its functional configurations, including muscle architecture, length, contraction velocity, and passive tension [31, 33, 63]: M M (7) F M ðt Þ ¼ F max F LV le ; ve M aM ðt Þ þ F PE le M M is the where F LV le ; ve M is the normalized force–length–velocity curve of the active CE, F PE le force–length curve of passive elastic elements, lM and vM are the muscle length and velocity, respectively, and the upper tildes indicate their normalization by peak-force length and maximum velocity, respectively [33]. The joint angles representing the relative orientations between adjacent body segments serve as the generalized coordinates in joint space [29, 41–43, 46–49, 51]. If the whole-body configuration can be represented using n revolute joints (n DOFs), the joint space is composed of vectors of n generalized coordinates (qi, i = 1, …, n) as:

½q1 ; …; qn T ∈ Rn

(8)

where minor translational motions that exist in actual anatomical joints are considered negligible in this study. Although individual muscle soft tissue is flexible and infinite dimensional, the geometric constraints imposed on a musculotendon system through anatomical joint structure result in its configuration that is dependent on the muscle and musculotendon lengths. In addition, the kinematic constraints imposed among the associated musculotendon systems result in one effective DOF for each revolute joint (Figure 1). Therefore, each joint angle is a function of the associated muscle and musculotendon lengths: M M MT MT MT ði ¼ 1; …; nÞ qi ¼ qi lM 1 ; l2 ; …; lm ; l1 ; l2 ; …; lm

(9)

where lMT is the musculotendon length. The joint velocity is then calculated as the total time derivative of this joint angle function. Using the chain rule, the joint velocity at each joint can be written as a function of the associated muscle and musculotendon velocities:

Figure 1. Mapping of state variables and control inputs from muscle space to joint space, where u is muscle excitation. A diagram of a Hill-type musculotendon model is included, where θ is pennation angle. Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (6 of 22)


M M MT MT MT M M M MT MT MT q_ i ¼ q_ i lM ði ¼ 1; …; nÞ (10) 1 ; l2 ; …; lm ; l1 ; l2 ; …; lm ; v1 ; v2 ; …; vm ; v1 ; v2 ; …; vm

where dots indicate time derivatives and vMT is the musculotendon velocity. These functions formulate the mapping (coordinate transformation) of obtaining the state variables in joint space from muscle space (Figure 1). Because n is much smaller than m, the joint space is uniquely determined from the muscle space, but the reverse is not true. The muscle-induced actuator torque τ(t) at each joint is calculated as the vector sum of the moments of the associated muscle forces about the joint axis of rotation. As seen earlier, the anatomical joint structure and the muscles’ contact and wrapping with the surrounding tissues, which also cause curved muscle geometry unlike the common straight-line model assumption, impose geometric constraints on the mapping of joint angles from the associated muscle and musculotendon configurations. Consequently, the orientation and the position of the point of application (including the varying moment arm length) of each muscle force relative to the rotational axis vector are all functions of muscle and musculotendon lengths [29]. Therefore, each joint actuator torque is a function of the associated muscle forces and lengths and musculotendon lengths: M M M M M MT MT MT ði ¼ 1; …; nÞ (11) τi ¼ τi FM 1 ; F 2 ; …; F m ; l1 ; l2 ; …; lm ; l1 ; l2 ; …; lm

This function represents the mapping of obtaining the control inputs in joint space (joint torque) from muscle space (muscle forces) (Figure 1). It should be noted from the variables in this equation that through the proposed mapping, the activation (u and aM(t)), contraction (FM, aM(t), lM, and vM; Equation (7)), and musculotendon (including SE, PE, and tendon compliance and CE contraction) dynamics (Figure 1), which govern the control of individual muscles, are incorporated into the actuator torque at each joint as the resultant of multiple associated muscles. In this approach, the force–length–velocity properties in muscle space [31, 33] are also inherently incorporated into joint space as torque–angle–velocity properties, which are generally monotonic and smooth [51]. To avoid the complexities and indeterminacies in obtaining joint actuator torque from multiple musculotendon dynamics, the joint actuator torque serves as a generalized torque that is directly associated with the corresponding generalized coordinate (i.e., joint angle) (Figure 2). The jointspace dynamics includes the resultant effects of actuator torque, inertia, internal reaction, external loads, and passive and active rotational stiffness (due to SE, PE, tendons, and other surrounding tissues). The remainders of the muscle forces that correspond with (zero) moment equilibrium contribute to the active stiffness as muscle cocontraction (i.e., coactivation). Cocontraction can be quantified as an index or ratio of the muscle forces from electromyography (EMG) signals of major agonist and antagonist muscle groups or their relative joint moments [64, 65]. Although it is inefficient from an energy standpoint, certain levels of cocontraction are necessary to stabilize a joint during some activities [9, 17, 64, 66]. Unlike in muscle space, the contributions of agonistic and antagonistic muscle activations are not distinguishable in joint space. In this study, the cocontraction at a joint

Figure 2. Leg segments above and below the knee illustrating mapping from muscle-space model to dynamically equivalent joint-space model. Copyright © 2015 John Wiley & Sons, Ltd.



e02721 (7 of 22)

is divided into two components: torque dependent and torque independent. Torque-dependent cocontraction provides structural and motor control stability [66, 67] through active joint stiffness [67–70]. Because the structure and motor control for increased joint actuator torque should be stabilized through increased active joint stiffness [69, 70], the activation level associated with torque-dependent cocontraction is approximately proportional to the joint actuator torque. Torque-independent cocontraction is due to muscle activation that is not essential for the given activity, and its activation level is arbitrary and independent of the resultant joint actuator torque. In line with the joint-space concepts of generalized coordinates and torques, the generalized activation is defined in this study as a measure of the effective total activations of the associated muscles at each joint. The aforementioned calculation of actuator torque from the associated muscle forces at a joint and the contraction dynamics (Equation (7)) for each muscle indicate that, apart from cocontraction, the generalized activation has a monotonically increasing relationship with the joint actuator torque. This enables the approximation of the generalized activation as a linear polynomial function of the joint actuator torque, which is also in analogy with the contraction dynamics (Equation (7)) in muscle space. Because the calculation of actuator torque also depends on the muscle and musculotendon lengths as discussed earlier, the coefficients in the linear function of generalized activation depend on the joint is, the angle. The kinematic parameter functions, that M LV eM PE eM of the active CE and the force–length curve F l ; ve l force–length–velocity curve F of passive elastic elements, in the contraction dynamics, together with the aforementioned kinematic mappings, indicate that the coefficients of the generalized activation linear function depend on the joint angle and velocity. In addition, because the maximum force of a given muscle depends on the muscle length and contraction velocity (as well as physiological cross-sectional area) [17, 31], the maximum (across angles and velocities) actuator torque τ max at a given joint should be adjusted depending on the joint angle and velocity [51]. Therefore, the linear-approximated contribution of the joint actuator torque to the generalized activation a(t) depends on the joint angle, joint velocity, and τ max: _ τ max Þfjτ ðt Þj þεðq; q; _ τ max Þg þ ca2 ðtÞ aðt Þ ¼ ca1 ðq; q;

(12)

where the 0th-order coefficient term in the linear function is divided according to the dependency on the joint-space state variables and the muscle force–length–velocity curve characteristics are taken into account in the lumped coefficient functions ca1 and ca2 such that 0 ≤ a(t) ≤ 1 for ∀ t. The function ca1 also includes the effects of torque-dependent cocontraction. The independent coefficient function ca2 includes both the activation component for active stiffness during torque-independent cocontraction _ τ max Þ corresponds to the minimum and the passive torque from elastic elements. The function εðq; q; torque-dependent cocontraction that exists even when τ(t) is zero.

2.3. Generalized heat dissipation rates A common assumption is that depending on the energetic characteristics, the heat that active skeletal muscles dissipate can be divided into activation–maintenance and shortening–lengthening heat [33, 61, 62, 71]. 2.3.1. Activation–maintenance heat rate. The activation–maintenance heat dissipation rate for a given muscle is proportional to muscle activation [60, 72, 73]. The generalized activation– am maintenance heat dissipation rate Q_ i of each joint can be derived using the generalized activation function: am _ τ max jτ i ðt Þj þhamε q; q; _ τ max þ ham1 ðt Þ ði ¼ 1; …; nÞ (13) Q_ ðt Þ ¼ ham q; q; i

i

i

i

i

i

where ham is the coefficient function associated with the joint actuator torque and torque-dependent i is the coefficient associated with minimum torque-dependent cocontraction at cocontraction, hamε i is the time-varying coefficient associated with the active stiffness due to zero torque, and ham1 i torque-independent cocontraction. Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (8 of 22)


2.3.2. Shortening–lengthening heat rate. The shortening–lengthening heat dissipation rate can be approximated as proportional to the activation level and the shortening–lengthening speed of the muscle CE [60, 61, 73, 74]. Because of the redundancy in the combined CE speeds of agonist and antagonist as well as the respective speeds of the tendon, SE, and PE, the CE speed cannot be uniquely mapped as a joint velocity. Instead, the muscle CE speed should be mapped to joint-velocitydependent and joint-velocity-independent terms. Using a similar approach as in the muscle-to-joint activation mapping, the joint-velocity-dependent term is approximately proportional to the joint velocity, and the joint-velocity-independent term is arbitrary and associated with cocontraction. Therefore, muscle CE speed can be approximated as a linear function (using similar process as in the generalized activation) of joint velocity. The generalized shortening–lengthening heat dissipation rate as the product of the activation function, muscle CE speed, and proportionality coefficients is sl sl2 _ τ max _ τ max _ i ðtÞj þ hsl3 Q_ i ðt Þ ¼ hsli q; q; þ hsl1 jτ i ðt Þq_ i ðt Þj þ hslε i i q; q; i i ðt Þjτ i ðt Þj þ hi ðt Þjq i ðt Þ ði ¼ 1; …; nÞ (14) sl2 sl3 where hsl1 i , hi , and hi are the time-varying coefficients associated with the active stiffness due to torque-independent cocontraction. The coefficient hslε represents the minimum torque-dependent i cocontraction. The coefficient hsli is associated with joint mechanical power, and separate functions are used to address the different energetic effects of positive and negative work [17, 40, 75]:

_ τ max _ τ max _ τ max _ τ max hsli q; q; q; q; ðif τ i ðt Þq_ i ðtÞ ≥ 0Þ; hsli q; q; q; q; ðif τ i ðt Þq_ i ðtÞ < 0Þ (15) ¼ hslþ ¼ hsl i i i i i i

_ am and the last three terms in 2.3.3. Cocontraction heat rate. The coefficients of the last term in Q i sl Q_ i , which are all associated with torque-independent cocontraction, can be regarded as voluntarily arbitrary variables. This is in contrast to the rest of the coefficients, which are system parameters specific to a given subject. Because of the nature of the mapping from muscle space to joint space, it is reasonable to re-group the arbitrary terms as a single independent variable that represents the cc generalized torque-independent cocontraction heat dissipation rate Q_ i ðt Þ: cc sl2 _ i ðtÞj þhsl3 ðtÞ þ hsl1 Q_ i ðt Þ ¼ ham1 i i ðt Þjτ i ðt Þj þhi ðt Þjq i ðt Þ ði ¼ 1; …; nÞ

(16)

Likewise, the coefficients that represent the minimum torque-dependent cocontraction can be regrouped as a single function: ε _ τ max _ τ max _ τ max ¼ hamε þ hslε ði ¼ 1; …; nÞ (17) q; q; Q_ i q; q; i i i i q; q; i In these forms, the redundant number of variables can be reduced significantly, which is another benefit of the joint-space approach. 2.4. Total MEE model in joint space _ mus at time t can be calculated as the summation Using the generalized coordinates, we obtain that W _ mus of the rate of work carried out by the joint actuator torque (or moment) for each DOF W i ðt Þ ¼ τ i ðt Þq_ i ðt Þ (i = 1, …, n), which is also called the rate of joint work [58, 59]. The total heat dissipation rate is the sum of the earlier components at each joint and can be incorporated into the total MEE rate (W) (Figure 3): met _ mus ðt Þ þ Q_ mus ðt Þ þ BMR E_ ðt Þ ¼ W n n n X X X _ τ max _ τ max ¼ τ i ðt Þq_ i ðt Þ þ ham q; q; hsli q; q; jτ i ðt Þj þ jτ i ðt Þq_ i ðt Þj i i i i¼1

þ

n X

i¼1

ε _ τ max þ Q_ i q; q; i

i¼1

Copyright © 2015 John Wiley & Sons, Ltd.

n X

i¼1

(18)

cc Q_ i ðt Þ þ BMR

i¼1



e02721 (9 of 22)

Figure 3. Overall structure of energy transformation and transfer for the human body. The gray boxes indicate the components that are included in the current model. Quantities in the white boxes, although not currently included in the model, are included here for complete system definition. MEE, metabolic energy expenditure; BMR, basal metabolic rate.

From dynamic systems perspectives, this equation represents the numerical model of the ε cc MEE rate in terms of the state variables (qi and q_ i), control inputs (τ i, Q_ i , and Q_ i ), and system sl parameters (ham i , hi , and BMR). The time-varying state variables and control inputs depend on the specific motion kinematics and dynamics and thus can incorporate any transient or steady-state motion. According to the second law of thermodynamics, any heat transfer to mus the body cannot be stored or recycled as MEE. Therefore, the components in Q_ ðt Þ are always non-negative, where the absolute values reflect the contributions of both positive and negative variables. It should be noted from the earlier derivations that the joint variables, velocities and torques in this model, and experimentally measured inverse-dynamics data (see succeeding text) are the results of the activation–contraction–musculotendon dynamics of the associated muscles.

2.5. Experiments and data processing 2.5.1 Parameter estimation. While the numerical MEE model is derived analytically using the first principles and muscle heat characteristics, its model parameter values should be estimated from experimental data. Twenty-one healthy subjects (Table I) were studied from a publicly available experimental walking data set [76]. Each subject was instrumented with passive optical markers and was captured during several trials at a preferred walking speed, while force plates captured ground reaction forces. The motion capture system (Motion Analysis Corporation, Santa Rosa, CA, USA) and force plate (Bertec Corporation, Columbus, OH, USA) data of each of the four to five trials from each subject were processed (total number of trials N = 102) with the MotionMonitor software [76], and the exported joint velocities and torques were used as inputs to a custom MATLAB (MathWorks, Natick, MA, USA) algorithm (see succeeding text). The available joint torques were calculated from the built-in morphometric data and inversedynamics modules. All walking experiments provided three-dimensional lower-body kinematics and kinetics with n = 10 DOFs consisting of one ankle, one knee, and three hip DOFs at each leg (Figure 4). As for metabolic data, the MEE values were estimated using a common walking energy expenditure formula from the American College of Sports Medicine (ACSM) [5, 6, 12], and the BMR was approximated as a function of basic subject-specific parameters (mass, height, age, and gender) using the Mifflin–St Jeor equation [77]. Because the focus of this initial study is rigorous derivation of numerical models, the purpose of using these approximated metabolic data was to illustrate the proposed method of parameter estimation (details in the succeeding text) for proof of concept. Estimating the parameter values more accurately from indirect calorimetry with a large number of subjects and trials will be pursued in a future study as described in the following. However, it should be noted that all the kinematic and kinetic data were directly measured from experiments. Copyright © 2015 John Wiley & Sons, Ltd.



11 subjects 10 subjects S23 S12 IC_M IC_F

MEE, metabolic energy expenditure.

Validation

Estimation

Subjects M F M F M F

Sex 1.79 ± 0.09 1.69 ± 0.04 1.79 1.63 1.75 1.56

Height (m) 82.1 ± 9.7 61.1 ± 4.5 68.5 52.6 75 52.2

Mass (kg) 23 ± 4 21 ± 1 21 20 32 23

Age 1.33 ± 0.10 1.34 ± 0.14 1.33 1.18 0.94 0.76

Trial 1

1.34 ± 0.15 1.35 ± 0.13 1.33 1.21 1.10 1.03

Trial 2

1.33 ± 0.12 1.36 ± 0.13 1.27 1.36 1.31 1.26

Trial 3

Speed per trial (m/s)

1.34 ± 0.11 1.37 ± 0.15 1.32 1.28 1.46 1.42

Trial 4

Table I. Descriptive statistics for all subjects used in parameter estimation and validation of the MEE model (mean ± standard deviation).

1.31 ± 0.12 1.35 ± 0.11 1.34 1.30 1.65 1.48

Trial 5

e02721 (10 of 22) J. H. KIM AND D. ROBERTS



e02721 (11 of 22)

Figure 4. Marker set used for three-dimensional motion analysis (source: MotionMonitor [76]) and snapshot of walking with global and local coordinate frames.

2.5.2. Validation. Two separate experiment groups of subjects (Table I) were used for initial validation and illustration of the model. In the first group, the ACSM MEE values were evaluated for two subjects (S12 and S23) that were excluded from the set used for parameter estimation. In the second group, MEE was measured through indirect calorimetry for two subjects (IC_M and IC_F) during walking. Informed consent was obtained, and the study was approved by a local Institutional Review Board. The BMR was obtained from indirect calorimetry during quiet sitting as the average of the last three minutes over a six minute interval. To provide variations between subjects’ trials, the protocol consisted of walking (on a treadmill for indirect calorimetry then over ground for kinematic and kinetic data collection with a pair of photogates to ensure consistent speeds) with rest periods of 2 min between five different speed levels: 70%, 85%, 100%, 115%, and 130% of the preferred walking speed. The short duration of each trial and the sufficient resting time between the trials of a given subject ensure that the effects of fatigue, sweating, and the _ 2 (L/min) was measured for each minute variations of the BMR are minimized. Average VO during the 5-min period at each speed, and the average of the last three of these measurements was converted [12] to MEE at each speed. Whole-body kinematic data were collected using a 12-camera motion capture system (Motion Analysis Corporation), and ground reaction forces were measured concurrently using four built-in force plates (AMTI, Watertown, MA, USA). Data processing (including inverse dynamics and morphometric approximations) was completed in Visual3D software (C-Motion, Germantown, MD, USA), and all joint variables and torques were exported per stride for subsequent processing in MATLAB.

2.6. Model approximations To take into account the different energetic properties across different subjects and DOFs in comparable and uniform scales, the heat coefficients and cocontraction heat rates are normalized as the following dimensionless quantities (indicated with upper tilde): max _ cc Þ ham ecc _ q; q; τ am i i max max sl max e sl q; q; e _ ðtÞ ¼ Qi ðtmax _ τi _ _ ¼ ¼ h ; ; h τ q; q; τ (19) Q hi q; q; i i i i i q_ max fτ i q_ i g i where q_ max is the maximum voluntary angular velocity and fτ i q_ i gmax is the maximum voluntary i mechanical power at joint i for a given subject; in this study, the respective maximum values in the available walking data are used. In other words, the heat coefficients are non-dimensionalized for each DOF and for each subject. As a simplifying assumption, each of these three dimensionless quantities is approximated as a constant function with respect to the DOF (including right–left symmetry) and subject. In addition, the dimensionless quantities for activation–maintenance and shortening– lengthening heat are assumed to be time invariant, which is a valid approximation for use in normal walking in which the variations due to the activation dynamics and the force–length–velocity properties of muscle are negligible [63]. Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (12 of 22)


X ε The Q_ i is assumed to be negligible as compared with other terms. The τ max at each DOF is assumed to be consistent across each group of male and female subjects, and is obtained as the mean value of those found in previous studies [51, 78]. To approximate the whole-body model MEE by using the given lower-body kinematic and kinetic data, a scaling factor of 1.05, in accordance with the reduced MEE (approximately 6%) during walking experiments with external stabilization without arm swing [79], is multiplied to the lower-body work and heat terms in the MEE model. Also, because of the short duration of each trial and sufficient resting time between the trials of a given subject, the effects of fatigue, sweat loss, and the variations of the BMR during each trial and across trials are assumed to be negligible. 2.7. Constrained multi-objective optimization for parameter estimation An optimization algorithm is used to estimate the heat coefficient functions that are unique system cc parameters of a given subject and Q_ i ðt Þ that are arbitrarily controlled state variables. Because the cc arbitrary distribution of Q_ i ðt Þ to each DOF cannotX be determined uniquely without additional cc experimental data from EMG, only the summation Q_ i ðt Þ as a lumped term at each time is determined from the algorithm. The cost function has two components. First, the differences between the model MEE values and the ACSM MEE values from the earlier experiments are minimized as a form of least squares problem. Because only the time-average rate or the total amount of MEE (Ēmet) can be obtained from experimental measurement or existing estimation methods [6, 14], the proposed MEE rate is X integrated over the time duration T and implemented into cc the average residual. Second, the total Q_ i ðt Þ is minimized with the condition of similar time-profile patterns across different subjects, which is a valid assumption for a given natural task such as normal walking [17, 64]. The multi-objective cost function is h i2 !2 met _ met ðt Þdt n E E ∫ X X T cc am slþ sl 1 1 cc e k ∫T þ J he i ; he i ; he i ; Q_ i ðt Þ ¼ Q_ i ðt Þdt ðk ¼ 1; …; N Þ (20) 2 2 N i¼1

where the scale of one-half and the quadratic function in the second term are used for dimensional consistency with the least squares term. The cost function is minimized in a constrained multi-objective optimization algorithm, from which the dimensionless heat quantities and the non-dimensionalized torque-independent cocontraction heat rates are estimated. Based on the second law of thermodynamics, several inequality constraints representing the irreversibility of the MEE and heat dissipation are imposed. The non-negativity of heat dissipation implies cc _ τ max _ τ max _ τ max ≥ 0; hslþ ≥ 0; hsl ≥ 0; and Q_ i ðt Þ ≥ 0 ði ¼ 1; …; n for ∀t Þ (21) ham q; q; q; q; q; q; i i i i i i

In addition, because the expended metabolic energy cannot be recharged by negative muscle work, the net MEE (excluding the BMR [18]) at each DOF should be non-negative at each time instant: _ τ max _ τ max q; q; (22) τ i ðt Þq_ i ðt Þ þ ham jτ i ðt Þj þ hsli q; q; jτ i ðt Þq_ i ðt Þj ≥ 0 ði ¼ 1; …; n for ∀t Þ i i i The algorithm is implemented into the fmincon subroutine in MATLAB for general and reliable constrained nonlinear optimization on a 1.8-GHz processor with 4.0-GB RAM. The resulting optimization variables are used to complete the MEE model terms, which are validated against separate experiments. Note that under the aforementioned model approximations, the heat coefficients are time-invariant system parameters that are independent of motion, and thus, their estimated values from steady-state normal walking experiments can be used in the evaluation of general tasks (along with the associated time-varying state variables and control inputs) that are transient or non-periodic. Copyright © 2015 John Wiley & Sons, Ltd.



e02721 (13 of 22)

3. RESULTS The speed (1.34 ± 0.13 m/s), stride length (1.42 ± 0.10 m), and gait cycle time (T = 1.06 ± 0.04 s) determined from the experimental data confirm that these gait parameters are within normal ranges for healthy adults. The walking data show characteristic patterns [17] for the two main inputs to the model: joint velocity and torque, as well as the product of these terms (mechanical power) (Figure 5). The dimensionless heat quantities and the total torque-independent cocontraction heat are estimated from the constrained multi-objective optimization algorithm using the experimental data as inputs into the model (Table II). Independent experiments with two subjects (S12 and S23) using ACSM formula and two subjects (IC_M and IC_F) using indirect calorimetry were used to validate the model. The mean of the time-averaged model MEE rates for all five trials in all four validation subjects are within 3.6 ± 3.6% absolute error of the validation values. These results are also cross-compared with each other for inter-subject correlation analysis (Figure 6) where a significant correlation is noted (R2 = 0.97, p < 0.0005). The model also enables the prediction of instantaneous MEE rate as a function of time (Figure 7). Here, the gait cycle is defined as right-heel contact to the next right-heel contact and includes double-support (DS) and single-support phases when referencing both legs, and stance and swing phase when referencing a single leg. The two DS phases account for just 28% of the gait cycle time but are responsible for 49% of the model MEE per gait cycle on average across all validation subjects’ trials. The results of the model can also be used to illustrate the relative contributions of the three main _ mus , Q_ mus , and BMR) and the relative contributions of the terms components of MEE rate ( W X am X sl X cc mus component included in the Q_ Q_ , Q_ , and Q_ . The first trial of subject S12 is i

i

i

used as an illustrative example (Figure 8).

Figure 5. Right leg joint torque (top), angular velocity (middle), and mechanical power (bottom) profiles for main DOFs for all subjects and trials included in the model and validation subjects S12 and S23. The thick line represents the mean, and the dotted lines enclosing the shaded gray area represent ±1 standard. Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (14 of 22)


Table II. Results of optimization implemented to estimate dimensionless heat quantities and torque-independent cocontraction heat rate. Parameter (i = 1, …, n) am _ τ max he i q; q; i slþ e _ τ max hi q; q; i sl e _ τ max hi q; q; i e_ cc Time-average ∑Q i

Value 0.054 0.283 1.423 0.004

Figure 6. Plot showing inter-subject correlation between time-averaged model metabolic energy expenditure (MEE) rates and the validation MEE rates from the American College of Sports Medicine walking equation (S12 and S23) and indirect calorimetry (IC_M and IC_F).

Figure 7. Whole-body instantaneous model metabolic energy expenditure (MEE) rates of one stride at preferred walking speeds of validation subjects versus % gait cycle. Horizontal lines indicate the timeaveraged model MEE rates. Copyright © 2015 John Wiley & Sons, Ltd.



e02721 (15 of 22)

Figure 8. Components of model metabolic energy expenditure (MEE) rate for subject S12, trial 1.

The instantaneous efficiency as a function of time is calculated for the validation subjects as the ratio of mechanical power to the gross (BMR included) model MEE rate at each time step (Figure 9). In addition, for comparative purposes, different forms of mean efficiencies [18–20] are calculated using time-averages from the model and compared with their respective counterparts in the literature (Figure 9). The mean efficiencies as the ratio of total, positive, and negative joint work to the time-averaged gross model MEE from the validation subjects are 0.05, 0.29, and 0.34, respectively. The mean positive work efficiencies as the ratio of positive joint work to the sum of net or gross MEE and negative work magnitude [19, 20] calculated from the model for the validation subjects are (mean ± standard deviation) 0.22 ± 0.02 and 0.30 ± 0.02 for gross and net, respectively. Two existing mechanical measures [80] are also calculated to compare and validate the MEE model results from additional perspectives: total body segmental work (W/kg) and sum of the normalized absolute moment impulses per second acting on the joints of the lower extremities. These mechanical measures and their variations are commonly used to indicate the energy expenditure of walking, with or without direct relation with MEE [1, 11, 47], because they are relatively straightforward to calculate from experimental data. The model MEE rate and the two mechanical measures for the validation experiments are compared along with experimental data from Burdett et al. [80] (Figure 10). The correlation between the total body segmental work and

Figure 9. Gross model efficiency versus % gait cycle and analogous mean model total, positive, and negative joint work efficiencies for subject S12, trial 1. Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (16 of 22)


Figure 10. Total body segmental work comparison (left) and sum of normalized moment impulse comparison (right) between data from Burdett et al. [80] (circles are shown only for comparable walking speeds, but line of best fit and R2 value are for full data set) and the current model (triangles) from all validation subjects.

the model MEE rate is significant in both the data from Burdett et al. [80] (R2 = 0.63, p < 0.0005) and the model (R2 = 0.89, p < 0.0005). Similarly, the correlation between the sum of normalized moment impulse and the model MEE rate is significant in both the data from Burdett et al. [80] (R2 = 0.72, p < 0.0005) and the model (R2 = 0.20, p = 0.048). 4. DISCUSSION This study introduces a general numerical model of MEE in joint space and the estimation of the model parameters from a constrained multi-objective optimization algorithm using normal walking data. Both the discrete and continuous (Figure 5) descriptive results of the experimental kinematic/kinetic data are consistent with comparable literature that shows similar results for healthy adults [2, 17, 80]. While the reduction of the total number of joint-space parameters to am slþ sl hi , and e hi (Table II) is based on the aforementioned the three dimensionless quantities e hi , e assumptions (i.e., constant–function approximation with respect to DOF and subject, and time on subject and DOF are incorporated through invariance), the dependencies of the resulting ham i the maximum voluntary angular velocity at each joint for a given subject. 4.1. Model accuracy and validity The model MEE rates evaluated from independent tests with separate validation subjects are within 3.6 ± 3.6% absolute error relative to the validation values determined through the ACSM walking formula (S12 and S23) and indirect calorimetry experiments (IC_M and IC_F). The inter-subject correlation for the model and validation MEE values of all four validation subjects are significant (R2 = 0.97 and p < 0.0005; Figure 6). It is worth noting that the data points for subjects IC_M and IC_F demonstrate correlation over a larger range of MEE because they represent MEE at five different walking speed levels. This demonstrates the proof-of-concept ability of the proposed model to consistently predict reliable MEE values across different ranges of subjects and walking speeds. The possible sources of errors include not only the approximations used in the proposed model derivation but also the inaccuracies in the respective validation MEE values. Some of the errors in subjects S12 and S23 are attributable to the approximations used in the ACSM formula [5, 6, 12] and using estimated BMR as part of the model, because the Mifflin–St Jeor equation [77] itself has been shown to only be accurate to within ±10% of the measured value [81]. The errors in the indirect calorimetry subjects (IC_M and IC_F) may be due to the assumptions inherent in the indirect calorimetry measurement method, inaccuracies in the measurement equipment, and Copyright © 2015 John Wiley & Sons, Ltd.



e02721 (17 of 22)

_ 2 to MEE conversion [12]. These errors are also consistent with the approximations used in the VO _ the trial-to-trial VO2 difference for a subject’s given level of walking speed. Nevertheless, in comparison with existing muscle-space models that consistently overestimate MEE by up to 50% [2, 38, 39], the errors from this joint-space model are relatively low and not consistently positive. The improved accuracy in the proposed model is likely due to the estimations of the joint-space parameters directly from active human motion data, as opposed to the in vitro parameters that most muscle-space models are based on. 4.2. Instantaneous model MEE rates The instantaneous model MEE rates as functions of time (Figure 7) can only be validated indirectly, because experimental time-profile data of whole-body MEE rates for direct comparison does not exist in the current literature. For this purpose, various existing measures with the associated inherent assumptions can be used as discussed in the following. For instance, the mean model MEE rates of the validation subjects match well (Figure 6) with the respective validation values, as discussed earlier, and the profiles are similar between the subjects, with periodic attributes as would be expected. The similar patterns of the MEE rate curves across the subjects are consistent with the similar patterns of the MEE rate function’s independent variables, that is, joint velocity and torque profiles (Figure 5), across the subjects. It is also clear that the peaks in the MEE rate profile in all cases correspond with the DS phases of walking. Although the DS phases account for just 28% of the gait cycle time, they are responsible for 49% of the model MEE. This is consistent with previous studies demonstrating that the step-to-step transition that occurs in parallel with the DS phase is a major determinant of the metabolic cost of walking [2, 15]. Because mechanical power peaks during this DS phase when more muscles are engaged to support the body, it follows that metabolic cost is higher at this point because of the increase in muscular activity. Another validation aspect is the passive dynamic characteristics of normal human walking, in which energy expenditure and balance stability have incompatible relationship with each other. In this perspective, the current results are in line with the relatively high and low instability indices during single-support and DS phases, respectively [44, 82]. The model also provides the breakdown of the whole-body MEE rate into the various components (Figure 8). The similarity between the mechanical power and whole-body MEE rate profiles is consistent with the strong correlations (p < 0.001) between the experimental mean metabolic power and the anterior/posterior center-of-mass positive and negative mechanical power, respectively, during DS [25]. The difference between the two curves is explained by the BMR and the various components of the heat dissipation rates. The small positive values of the total torqueindependent cocontraction heat for the duration of X the gait cycle (Table II) results from the cc Q_ i ðtÞdt is minimized in one of the cost constrained multi-objective optimization, in which ∫T functions. This is consistent with the previously mentioned assumption that torque-independent cocontraction during normal walking is minimal and shows similar patterns of time profiles across different subjects [17, 64]. Again, it should be noted that this torque-independent cocontraction is distinct from the torque-dependent cocontraction, which is included in the models of activation– maintenance and shortening–lengthening heat dissipation rates. 4.3. Instantaneous metabolic efficiency The instantaneous metabolic efficiency during a dynamic task can be derived from the proposed model as a function of time, while the comparison with various mean efficiencies provides indirect validation of the model from additional perspectives (Figure 9). The mean model efficiency for total joint work (0.05) is near zero, which is consistent with analogous results obtained from total mechanical work for level walking in the literature [18]. The mean efficiency of positive joint work as determined by the model is 0.29, which is within the typical range reported in the literature (0.08–0.38 for whole body [21, 83] and 0.14–0.40 for isolated skeletal muscles and fibers [75, 84, 85]), and close to the uphill walking efficiency of 0.25 based on total mechanical work [7, 18, 24]. Although a comparable range of mean efficiency of negative joint work, other than that Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (18 of 22)


for downhill walking (1.2, [18]), does not exist in the literature, the model value of 0.34 is similar to the efficiency of eccentric contraction of an isolated skeletal muscle (0.35, [75]), and its higher absolute value than that of positive work is consistent with both whole body [17, 18] and individual muscle [75, 84, 85] studies. More rigorous comparison is provided by the mean positive work efficiencies calculated as the ratio of positive joint work to the sum of MEE and negative joint work magnitude [19, 20]. The mean of the time-averaged model positive work efficiencies with net MEE (0.30) and gross MEE (0.22) is close to the respective values (0.38 and 0.28) from the experiments and shows a similar relative difference [20]. These validations are more indirect in nature, because actual values of efficiency only exist in the literature for steady-state walking, but not within a single walking cycle or as a function of time. 4.4. Validation against mechanical measures The model MEE rate results can be further validated by comparison with other mechanical measures investigated in the literature. The total body segmental work and the model MEE rate from the validation subjects relative to the subjects in Burdett et al. [80] show significant correlations and a similar linear trend (Figure 10). A similarly significant correlation is also noted (Figure 10) between the sum of the normalized moment impulse and the normalized MEE rate for the model results and the experimental data [80]. These correlations between the two mechanical measures and the model MEE rate may be partially due to the structures of the respective formulations. The total body segmental work in Burdett et al. [80] resembles the mechanical power component in the current MEE model. The sum of the normalized moment impulse is also similar to the activation–maintenance heat dissipation component in the model. Therefore, the current MEE model inherently includes the naturally weighted contributions of these two mechanical measures. It should be noted that the inclusion of these similar terms in the proposed MEE model is a natural consequence of the derivation using the first law of thermodynamics and activation mapping to joint space. On the other hand, the two mechanical measures are time-averaged for a walking cycle, and the significant correlations are only shown to be valid under their inherent assumptions of normal walking, during which the magnitudes of joint torques and velocities consistently represent active and periodic movements. For general tasks, particularly those including durations of zero mechanical work (e.g., load carriage without movement), the total body segmental work may not be a reliable indicator of the MEE. In addition, the sum of the normalized moment impulse alone does not take into account the MEE associated with the time rate of change of configuration. The proposed MEE model is a function of time and can be applied to any general tasks, including those that are not periodic or have durations of zero mechanical work. 4.5. Limitations of the current model and future work This study focused on rigorous mathematical modeling of MEE and the method of its computational parameter estimation and evaluation. Although the results in this proof-of-concept study demonstrate the initial validity of the proposed formulations and methods, they were obtained from a model that was established under several aforementioned assumptions, which will be addressed in future work including indirect calorimetry measurements for a larger number of subjects and trials. First, although steady-state normal walking was chosen for validation because of the wealth of comparable literature, future work will extend experimental validation methods of the model to non-steady-state walking and other general tasks. Second, the dimensionless quantities in the generalized heat coefficients were approximated using constant terms for all DOFs and subjects, despite a possibility that these quantities differ depending on DOF and subject and contain higher order, time-dependent terms. More comprehensive estimation of subject-specific and DOF-specific heat coefficients and their sensitivity to the model will be investigated in future work. This will also include subject-specific modeling of cocontraction energetics at each DOF and their quantitative validations through EMG measurements. Third, the model has been validated against aerobic activities only, under the closed-system assumption of ideal contraction conditions such as ample amounts of resources that contraction depends on (e.g., calcium and Adenosine Tri-Phosphate (ATP)). The model does not currently account for the MEE because of anaerobic respiration that Copyright © 2015 John Wiley & Sons, Ltd.



e02721 (19 of 22)

happens during fatigue in extended activity and during high-power activities like running. Finally, future work should investigate an inverse-dynamics method that partitions the resultant joint torque into muscle-induced actuator torque and passive reaction moment (including elastic and frictional elements [86, 87]), and includes reliable contact models for more accurate calculation of the MEE components. 5. CONCLUSIONS A general joint-space numerical model of MEE was derived from the laws of thermodynamics applied to the human multibody dynamic system. The joint-space parameters estimated directly from active human subjects in a constrained multi-objective optimization algorithm provided the model with relatively high accuracy. Although the model was first validated in this study using various existing measures during normal walking, there are no inherent limitations on its ability to evaluate MEE during general tasks. This model enables real-time calculations of instantaneous MEE rate and efficiency as functions of time for transient as well as steady-state evaluations and focused analysis of different phases of general non-periodic tasks. The proposed MEE model can be applied with or without experimental data. When the experimental kinematic/kinetic data are provided as inputs (such as in typical gait analysis), the model will predict the MEE as output _ 2 measurements for complex non-periodic functional tasks that may not be without VO experimentally verifiable by traditional means. When experimental data are not provided, the model can be used for simulations, in which the kinematic/kinetic variables as well as the MEE are not specified and are computationally evaluated for various what-if scenarios. Although experimental measurements may not be completely replaced by model evaluations, predicted quantities can be used as strong complements and will substantially increase the reliability of the results and yield unique insight. While this study focused on rigorous mathematical modeling of MEE and the method of its computational parameter estimation and illustration, future work will address the simplifying assumptions in the current model with more extensive indirect calorimetry measurements and include more comprehensive validation methods. ACKNOWLEDGEMENTS

The authors would like to acknowledge the work of Amanda Setiawan and Zhenxiang Chen in processing the data and the research staff at the Leon Root, M.D. Motion Analysis Laboratory (Director: Dr. Howard Hillstrom) at the Hospital for Special Surgery for assistance with indirect calorimetry experiments. This work was supported by a US National Science Foundation (NSF) Graduate Research Fellowship (grant no. DGE-1104522), an NSF Career–Life Balance Supplement, an International Society of Biomechanics (ISB) Student Matching Dissertation Grant to D. R., and an NSF grant (no. CMMI-1436636) to J.H.K. REFERENCES 1. Srinivasan M. Fifteen observations on the structure of energy-minimizing gaits in many simple biped models. Journal of the Royal Society Interface 2011; 8(54):74–98. 2. Umberger BR. Stance and swing phase costs in human walking. Journal of the Royal Society Interface 2010; 7(50):1329–1340. 3. Waters RL, Mulroy S. The energy expenditure of normal and pathologic gait. Gait & Posture 1999; 9(3):207–231. 4. Long LL, Srinivasan M. Walking, running, and resting under time, distance, and average speed constraints: optimality of walk–run–rest mixtures. Journal of the Royal Society Interface 2013; 10(81):1–10. 5. Thompson WR, Gordon NF, Pescatello LS (Eds). ACSM’s Guidelines for Exercise Testing and Prescription. Lippincott Williams & Wilkins: Philadelphia, PA, USA, 2009. 6. Glass S, Dwyer GB (Eds). ACSM’s Metabolic Calculations Handbook. Lippincott Williams & Wilkins: Philadelphia, PA, USA, 2006. 7. Donelan JM, Li Q, Naing V, Hoffer JA, Weber DJ, Kuo AD. Biomechanical energy harvesting: generating electricity during walking with minimal user effort. Science 2008; 319(5864):807–810. 8. Au SK, Weber J, Herr H. Powered ankle-foot prosthesis improves walking metabolic economy. IEEE Transactions on Robotics 2009; 25(1):51–66. 9. Mian OS, Thom JM, Ardigò LP, Narici MV, Minetti AE. Metabolic cost, mechanical work, and efficiency during walking in young and older men. Acta Physiologica 2006; 186(2):127–139. Copyright © 2015 John Wiley & Sons, Ltd.


e02721 (20 of 22)


10. Brychta R, Wohlers E, Moon J, Chen K. Energy expenditure: measurement of human metabolism. IEEE Engineering in Medicine and Biology Magazine 2010; 29(1):42–47. 11. Ackermann M, van den Bogert AJ. Optimality principles for model-based prediction of human gait. Journal of Biomechanics 2010; 43(6):1055–1060. 12. de V Weir JB. New methods for calculating metabolic rate with special reference to protein metabolism. Journal of Physiology 1949; 109(1–2):1–9. 13. Seale JL, Rumpler WV, Conway JM, Miles CW. Comparison of doubly labeled water, intake-balance, and directand indirect-calorimetry methods for measuring energy expenditure in adult men. American Journal of Clinical Nutrition 1990; 52(1):66–71. 14. Garg A, Chaffin DB, Herrin GD. Prediction of metabolic rates for manual materials handling jobs. American Industrial Hygiene Association Journal 1978; 39(8):661–674. 15. Donelan JM, Kram R, Kuo AD. Mechanical work for step-to-step transitions is a major determinant of the metabolic cost of human walking. Journal of Experimental Biology 2002; 205(23):3717–3727. 16. Neptune RR, van den Bogert AJ. Standard mechanical energy analyses do not correlate with muscle work in cycling. Journal of Biomechanics 1997; 31(3):239–245. 17. Winter DA. Biomechanics and Motor Control of Human Movement. Wiley: Hoboken, NJ, USA, 2009. 18. Margaria R. Positive and negative work performances and their efficiencies in human locomotion. Internationale Zeitschrift für Angewandte Physiologie 1968; 25:339–351. 19. Prilutsky BI. Work, energy expenditure, and efficiency of the stretch-shortening cycle. Journal of Applied Biomechanics 1997; 13(4):466–470. 20. Umberger BR, Martin PE. Mechanical power and efficiency of level walking with different stride rates. Journal of Experimental Biology 2007; 210(18):3255–3265. 21. Farris DJ, Sawicki GS. The mechanics and energetics of human walking and running: a joint level perspective. Journal of the Royal Society Interface 2012; 9(66):110–118. 22. Sousa ASP, Silva A, Tavares JMRS. Biomechanical and neurophysiological mechanisms related to postural control and efficiency of movement: a review. Somatosensory and Motor Research 2012; 29(4):131–143. 23. Woledge RC, Curtin NA, Homsher E. Energetic Aspects of Muscle Contraction. Academic Press: Waltham, MA, USA, 1985. 24. Ruina A, Bertram JEA, Srinivasan M. A collisional model of the energetic cost of support work qualitatively explains leg sequencing in walking and galloping, pseudo-elastic leg behavior in running and the walk-to-run transition. Journal of Theoretical Biology 2005; 237(2):170–192. 25. Sousa AS, Silva A, Tavares JMR. Interlimb relation during the double support phase of gait: an electromyographic, mechanical and energy-based analysis. Proceedings of the Institution of Mechanical Engineers. Part H, Journal of Engineering in Medicine 2013; 227(3):327–333. 26. Hall C, Figueroa A, Fernhall B, Kanaley JA. Energy expenditure of walking and running: comparison with prediction equations. Medicine and Science in Sports and Exercise 2004; 36(12):2128–2134. 27. Van Bolhuis BM, Gielen CCAM. A comparison of models explaining muscle activation patterns for isometric contractions. Biological Cybernetics 1999; 81(3):249–261. 28. Ting LH, Chvatal SA, Safavynia SA, Lucas McKay J. Review and perspective: neuromechanical considerations for predicting muscle activation patterns for movement. International Journal for Numerical Methods in Engineering 2012; 28(10):1003–1014. 29. Sherman, MA, Seth, A, and Delp, SL, 2013, What is a moment arm? Calculating muscle effectiveness in biomechanical models using generalized coordinates. Proc. ASME Int. Design Engineering Technical Conf. and Computers and Information in Engineering Conf. (IDETC/CIE), Portland, OR, USA. 30. Modenese L, Phillips ATM, Bull AMJ. An open source lower limb model: hip joint validation. Journal of Biomechanics 2011; 44(12):2185–2193. 31. Zajac F. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering 1989; 17(4):359–411. 32. Karduna AR, Williams GR, Iannotti JP, Williams JL. Kinematics of the glenohumeral joint: influences of muscle forces, ligamentous constraints, and articular geometry. Journal of Orthopaedic Research 1996; 14(6):986–993. 33. Millard M, Uchida T, Seth A, Delp SL. Flexing computational muscle: modeling and simulation of musculotendon dynamics. Journal of Biomechanical Engineering-Transactions of the ASME 2013; 135(2):021004–1–021004–11. 34. Geraldes DM, Phillips ATM. A comparative study of orthotropic and isotropic bone adaptation in the femur. International Journal for Numerical Methods in Engineering 2014; 30(9):873–889. 35. Walcott S. A differential equation model for tropomyosin-induced myosin cooperativity describes myosin–myosin interactions at low calcium. Cellular and Molecular Bioengineering 2013; 6(1):13–25. 36. Ma S, Zahalak GI. A distribution-moment model of energetics in skeletal muscle. Journal of Biomechanics 1991; 24(1):21–35. 37. Vologodskii A. Energy transformation in biological molecular motors. Physics of Life Reviews 2006; 3(2):119–132. 38. Anderson FC, Pandy MG. Dynamic optimization of human walking. Journal of Biomechanical EngineeringTransactions of the ASME 2001; 123(5):381–390. 39. Miller, RH, Ackermann, M, and van den Bogert, AJ, 2013, Revisiting the prediction of walking mechanics and energetics in three dimensions by minimum metabolic cost. Proc. Am. Society of Biomechanics Conf., Omaha, NE, USA.




e02721 (21 of 22)

40. Miller RH. A comparison of muscle energy models for simulating human walking in three dimensions. Journal of Biomechanics 2014; 47(6):1373–1381. 41. Yang J, Feng X, Xiang Y, Kim JH, Rajulu S. Determining the three-dimensional relation between the skeletal elements of the human shoulder complex. Journal of Biomechanics 2009; 42(11):1762–1767. 42. Hawkins D, Molé P. Modeling energy expenditure associated with isometric, concentric, and eccentric muscle action at the knee. Annals of Biomedical Engineering 1997; 25(5):822–830. 43. Robertson DGE, Winter DA. Mechanical energy generation, absorption and transfer amongst segments during walking. Journal of Biomechanics 1980; 13(10):845–854. 44. Mummolo C, Mangialardi L, Kim JH. Quantifying dynamic characteristics of human walking for comprehensive gait cycle. Journal of Biomechanical Engineering-Transactions of the ASME 2013; 135(9):091006–1–091006–10. 45. Lin Y-C, Kim HJ, Pandy MG. A computationally efficient method for assessing muscle function during human locomotion. International Journal for Numerical Methods in Engineering 2011; 27(3):436–449. 46. Xiang Y, Chung H-J, Kim JH, Bhatt R, Rahmatalla S, Yang J, Marler T, Arora JS, Abdel-Malek K. Predictive dynamics: an optimization-based novel approach for human motion simulation. Structural and Multidisciplinary Optimization 2010; 41(3):465–479. 47. Kim JH, Xiang Y, Yang J, Arora JS, Abdel-Malek K. Dynamic motion planning of overarm throw for a biped human multibody system. Multibody System Dynamics 2010; 24(1):1–24. 48. Kim JH, Yang J, Abdel-Malek K. Planning load-effective dynamic motions of highly articulated human model for generic tasks. Robotica 2009; 27(5):739–747. 49. Silder A, Besier T, Delp SL. Predicting the metabolic cost of incline walking from muscle activity and walking mechanics. Journal of Biomechanics 2012; 45(10):1842–1849. 50. Minetti AE, Alexander RM. A theory of metabolic costs for bipedal gaits. Journal of Theoretical Biology 1997; 186(4):467–476. 51. Frey-Law L, Laake A, Avin KG, Heitsman J, Marler T, Abdel-Malek K. Knee and elbow 3D strength surfaces: peak torque-angle-velocity relationships. Journal of Applied Biomechanics 2012; 28(6):726–737. 52. Ward-Smith AJ. A mathematical theory of running, based on the first law of thermodynamics, and its application to the performance of world-class athletes. Journal of Biomechanics 1985; 18(5):337–349. 53. Morton RH. The critical power and related whole-body bioenergetic models. European Journal of Applied Physiology 2006; 96(4):339–354. 54. Peronnet F, Thibault G. Mathematical analysis of running performance and world running records. Journal of Applied Physiology 1989; 67(1):453–465. 55. Carr, CE, 2005, The bioenergetics of walking and running in space suits. Ph.D. Dissertation, Massachusetts Institute of Technology. 56. Black WZ, Hartley JG. Thermodynamics. Harpercollins College Division: New York, NY, USA, 1991. 57. Langhaar HL. Energy Methods in Applied Mechanics. R.E. Krieger Publishing Co.: Malabar, Fla, 1989. 58. Zelik KE, Kuo AD. Human walking isn’t all hard work: evidence of soft tissue contributions to energy dissipation and return. Journal of Experimental Biology 2010; 213(24):4257–4264. 59. Sasaki K, Neptune RR, Kautz SA. The relationships between muscle, external, internal and joint mechanical work during normal walking. Journal of Experimental Biology 2009; 212(5):738–744. 60. Davy DT, Audu ML. A dynamic optimization technique for predicting muscle forces in the swing phase of gait. Journal of Biomechanics 1987; 20(2):187–201. 61. Bhargava LJ, Pandy MG, Anderson FC. A phenomenological model for estimating metabolic energy consumption in muscle contraction. Journal of Biomechanics 2004; 37(1):81–88. 62. Hill AV. The heat of shortening and the dynamic constants of muscle. Proceedings of the Royal Society B-Biological Sciences 1938; 126(843):136–195. 63. Anderson FC, Pandy MG. Static and dynamic optimization solutions for gait are practically equivalent. Journal of Biomechanics 2001; 34(2):153–161. 64. Falconer K, Winter DA. Quantitative assessment of co-contraction at the ankle joint in walking. Electromyography and Clinical Neurophysiology 1985; 25(2–3):135–149. 65. Knarr BA, Zeni JA, Higginson JS. Comparison of electromyography and joint moment as indicators of co-contraction. Journal of Electromyography and Kinesiology 2012; 22:607–611. 66. Da Fonseca ST, Vaz DV, de Aquino CF, Brício RS. Muscular co-contraction during walking and landing from a jump: comparison between genders and influence of activity level. Journal of Electromyography and Kinesiology 2006; 16(3):273–280. 67. Osu R, Franklin DW, Kato H, Gomi H, Domen K, Yoshioka T, Kawato M. Short- and long-term changes in joint co-contraction associated with motor learning as revealed from surface EMG. Journal of Neurophysiology 2002; 88(2):991–1004. 68. Milner T. Adaptation to destabilizing dynamics by means of muscle cocontraction. Experimental Brain Research 2002; 143(4):406–416. 69. Milner TE, Cloutier C, Leger AB, Franklin DW. Inability to activate muscles maximally during cocontraction and the effect on joint stiffness. Experimental Brain Research 1995; 107(2):293–305. 70. Hunter IW, Kearney RE. Dynamics of human ankle stiffness: variation with mean ankle torque. Journal of Biomechanics 1982; 15(10):747–752. 71. Hill AV. Production and absorption of work by muscle. Science 1960; 131(3404):897–903.



e02721 (22 of 22)


72. Hase K, Yamazaki N. Development of three-dimensional whole-body musculoskeletal model for various motion analyses. JSME International Journal Series C Dynamics, control, robotics, design and manufacturing 1997; 40(1):25–32. 73. Hatze H, Buys JD. Energy-optimal controls in the mammalian neuromuscular system. Biological Cybernetics 1977; 27(1):9–20. 74. Umberger BR, Gerritsen KGM, Martin PE. A model of human muscle energy expenditure. Computer Methods in Biomechanics and Biomedical Engineering 2003; 6(2):99–111. 75. Ryschon TW, Fowler MD, Wysong RE, Anthony A-R, Balaban RS. Efficiency of human skeletal muscle in vivo: comparison of isometric, concentric, and eccentric muscle action. Journal of Applied Physiology 1997; 83(3):867–874. 76. Wilson, JD, and Kernozek, T, “3D gait data by the MotionMonitor®” [online]. (Available from: http://www. innsport.com/related-products/data-sets/uw-l-gait-data-set.aspx). [Accessed: 10-Sep-2012]. 77. Mifflin M, St Jeor S, Hill L, Scott B, Daugherty S, Koh Y. A new predictive equation for resting energy expenditure in healthy individuals. American Journal of Clinical Nutrition 1990; 51(2):241–247. 78. Anderson DE, Madigan ML, Nussbaum MA. Maximum voluntary joint torque as a function of joint angle and angular velocity: model development and application to the lower limb. Journal of Biomechanics 2007; 40(14):3105–3113. 79. Ortega JD, Fehlman LA, Farley CT. Effects of aging and arm swing on the metabolic cost of stability in human walking. Journal of Biomechanics 2008; 41(16):3303–3308. 80. Burdett RG, Skrinar GS, Simon SR. Comparison of mechanical work and metabolic energy consumption during normal gait. Journal of Orthopaedic Research 1983; 1(1):63–72. 81. Frankenfield D, Roth-Yousey L, Compher C. Comparison of predictive equations for resting metabolic rate in healthy nonobese and obese adults: a systematic review. Journal of the American Dietetic Association 2005; 105(5):775–789. 82. Kim JH, Joo CB. Numerical construction of balanced state manifold for single-support legged mechanism in sagittal plane. Multibody System Dynamics 2014; 31(3):257–281. 83. Gaesser GA, Brooks GA. Muscular efficiency during steady-rate exercise: effects of speed and work rate. Journal of Applied Physiology 1975; 38(6):1132–1139. 84. Smith NP, Barclay CJ, Loiselle DS. The efficiency of muscle contraction. Progress in Biophysics and Molecular Biology 2005; 88(1):1–58. 85. He Z-H, Bottinelli R, Pellegrino MA, Ferenczi MA, Reggiani C. ATP consumption and efficiency of human single muscle fibers with different myosin isoform composition. Biophysical Journal 2000; 79(2):945–961. 86. Damm P, Dymke J, Ackermann R, Bender A, Graichen F, Halder A, Beier A, Bergmann G. Friction in total hip joint prosthesis measured in vivo during walking. PLoS ONE 2013; 8(11):e78373. 87. Bergmann G, Graichen F, Bender A, Kääb M, Rohlmann A, Westerhoff P. In vivo glenohumeral contact forces— measurements in the first patient 7 months postoperatively. Journal of Biomechanics 2007; 40(10):2139–2149.



Multibody model of the human upper extremity for fracture simulation.

Metabolic energy expenditure and terrain coefficients for walking on snow.

Resting metabolic rate of esophageal carcinoma patients: a model for energy expenditure measurement in a homogenous cancer population.

A dynamic model for human population growth.

Validation of a portable indirect calorimetry system for measurement of energy expenditure in sick preterm infants.

The adaptive metabolic response to exercise-induced weight loss influences both energy expenditure and energy intake.

A mathematical model of the human metabolic system and metabolic flexibility.

Use of a Wireless Network of Accelerometers for Improved Measurement of Human Energy Expenditure.

A Single Oral Administration of Theaflavins Increases Energy Expenditure and the Expression of Metabolic Genes.

Effect of protein overfeeding on energy expenditure measured in a metabolic chamber.

Resting energy expenditure (REE) in an old-old population: implications for metabolic stress.

A dynamic multibody model of the physiological knee to predict internal loads during movement in gravitational field.

Energy expenditure in caving.

Estimating energy expenditure.

Aging and energy expenditure.

Energy expenditure during bicycling.

Human resting energy expenditure in relation to dietary potassium.

Human resting energy expenditure in relation to dietary potassium.

Measurement Methods for Physical Activity and Energy Expenditure: a Review.

The measurement of energy expenditure.

1H-NMR-based metabolic analysis of human serum reveals novel markers of myocardial energy expenditure in heart failure patients.

[The evaluation of serum creatine phosphokinase for estimation of energy expenditure in human subjects].

organ energy expenditure to metabolic fluxes in mouse and human: experimental data integrated with mathematical modeling.

Numerical solution of a microbial growth model applied to dynamic environments.