Analytical and numerical analysis of inverse optimization problems: conditions of uniqueness and computational methods

One of the key problems of motor control is the redundancy problem, in particular how the central nervous system (CNS) chooses an action out of infinitely many possible. A promising way to address this question is to assume that the choice is made based on optimization of a certain cost function. A number of cost functions have been proposed in the literature to explain performance in different motor tasks: from force sharing in grasping to path planning in walking. However, the problem of uniqueness of the cost function(s) was not addressed until recently. In this article, we analyze two methods of finding additive cost functions in inverse optimization problems with linear constraints, so-called linear-additive inverse optimization problems. These methods are based on the Uniqueness Theorem for inverse optimization problems that we proved recently (Terekhov et?al., J Math Biol 61(3):423?C453, 2010). Using synthetic data, we show that both methods allow for determining the cost function. We analyze the influence of noise on the both methods. Finally, we show how a violation of the conditions of the Uniqueness Theorem may lead to incorrect solutions of the inverse optimization problem.     

A 'cheap' optimal control approach to estimate muscle forces in musculoskeletal systems

This paper shows a new method to estimate the muscle forces in musculoskeletal systems based on the inverse dynamics of a multi-body system associated optimal control. The redundant actuator problem is solved by minimizing a time-integral cost function, augmented with a torque-tracking error function, and muscle dynamics is considered through differential constraints. The method is compared to a previously implemented human posture control problem, solved using a Forward Dynamics Optimal Control approach and to classical static optimization, with two different objective functions. The new method provides very similar muscle force patterns when compared to the forward dynamics solution, but the computational cost is much smaller and the numerical robustness is increased. The results achieved suggest that this method is more accurate for the muscle force predictions when compared to static optimization, and can be used as a numerically 'cheap' alternative to the forward dynamics and optimal control in some applications.     

A Hybrid Optimization Method for Solving Bayesian Inverse Problems under Uncertainty

In this paper, we investigate the application of a new method, the Finite Difference and Stochastic Gradient (Hybrid method), for history matching in reservoir models. History matching is one of the processes of solving an inverse problem by calibrating reservoir models to dynamic behaviour of the reservoir in which an objective function is formulated based on a Bayesian approach for optimization. The goal of history matching is to identify the minimum value of an objective function that expresses the misfit between the predicted and measured data of a reservoir. To address the optimization problem, we present a novel application using a combination of the stochastic gradient and finite difference methods for solving inverse problems. The optimization is constrained by a linear equation that contains the reservoir parameters. We reformulate the reservoir model’s parameters and dynamic data by operating the objective function, the approximate gradient of which can guarantee convergence. At each iteration step, we obtain the relatively ‘important’ elements of the gradient, which are subsequently substituted by the values from the Finite Difference method through comparing the magnitude of the components of the stochastic gradient, which forms a new gradient, and we subsequently iterate with the new gradient. Through the application of the Hybrid method, we efficiently and accurately optimize the objective function. We present a number numerical simulations in this paper that show that the method is accurate and computationally efficient.     

A dynamic optimization technique for predicting muscle forces in the swing phase of gait

The muscle force sharing problem was solved for the swing phase of gait using a dynamic optimization algorithm. For comparison purposes the problem was also solved using a typical static optimization algorithm. The objective function for the dynamic optimization algorithm was a combination of the tracking error and the metabolic energy consumption. The latter quantity was taken to be the sum of the total work done by the muscles and the enthalpy change during the contraction. The objective function for the static optimization problem was the sum of the cubes of the muscle stresses. To solve the problem using the static approach, the inverse dynamics problem was first solved in order to determine the resultant joint torques required to generate the given hip, knee and ankle trajectories. To this effect the angular velocities and accelerations were obtained by numerical differentiation using a low-pass digital filter. The dynamic optimization problem was solved using the Fletcher-Reeves conjugate gradient algorithm, and the static optimization problem was solved using the Gradient-restoration algorithm. The results show influence of internal muscle dynamics on muscle control histories vis a vis muscle forces. They also illustrate the strong sensitivity of the results to the differentiation procedure used in the static optimization approach.     

Dynamic motion planning of 3D human locomotion using gradient-based optimization

Since humans can walk with an infinite variety of postures and limb movements, there is no unique solution to the modeling problem to predict human gait motions. Accordingly, we test herein the hypothesis that the redundancy of human walking mechanisms makes solving for human joint profiles and force time histories an indeterminate problem best solved by inverse dynamics and optimization methods. A new optimization-based human-modeling framework is thus described for predicting three-dimensional human gait motions on level and inclined planes. The basic unknowns in the framework are the joint motion time histories of a 25-degree-of-freedom human model and its six global degrees of freedom. The joint motion histories are calculated by minimizing an objective function such as deviation of the trunk from upright posture that relates to the human model's performance. A variety of important constraints are imposed on the optimization problem, including (1) satisfaction of dynamic equilibrium equations by requiring the model's zero moment point (ZMP) to lie within the instantaneous geometrical base of support, (2) foot collision avoidance, (3) limits on ground-foot friction, and (4) vanishing yawing moment. Analytical forms of objective and constraint functions are presented and discussed for the proposed human-modeling framework in which the resulting optimization problems are solved using gradient-based mathematical programming techniques. When the framework is applied to the modeling of bipedal locomotion on level and inclined planes, acyclic human walking motions that are smooth and realistic as opposed to less natural robotic motions are obtained. The aspects of the modeling framework requiring further investigation and refinement, as well as potential applications of the framework in biomechanics, are discussed.     

A contact model to simulate human–artifact interaction based on force optimization: implementation and application to the analysis of a training machine

Musculoskeletal multibody models are increasingly used to analyze and optimize physical interactions between humans and technical artifacts. Since interaction is conveyed by contact between the human body and the artifact, a computationally robust modeling approach for frictional contact forces is a crucial aspect. In this contribution, we propose a parametric contact model and formulate an associated force optimization problem to simultaneously estimate unknown muscle and contact forces in an inverse dynamic manner from a prescribed motion trajectory. Unlike existing work, we consider both the static and the kinetic regime of Coulomb’s friction law. The approach is applied to the analysis of a leg extension training machine with the objective to reduce the stress on the tibiofemoral joint. The uncertainty of the simulation results due to a tunable parameter of the contact model is of particular interest.     

The inverse dynamics problem of neuromuscular control

The myoskeletal inverse dynamics problem and the myocybernetic control inverse problem were investigated with respect to their ill-posedness. The first problem consists of finding from observed experimental motion and reaction force data the resultant muscle moments that generated the observed motion, while the second aims at finding the corresponding neural controls. It is shown that both problems belong to the class of incorrectly posed (ill-posed) problems that, by definition, do not possess unique solutions. To illustrate this point, results of a forward dynamics simulation of a comprehensive neuromusculoskeletal model of the human body are presented. These results demonstrate that fairly chaotic neural control perturbations have very little influence on the resulting motion trajectory, at least in the present example. While a regularization procedure may be applied to solve successfully the myoskeletal inverse dynamics problem, the myocybernetic control inverse problem is unsolvable. The latter fact has the important implication that, based on the somatosensory inputs it receives, the pars intermedia in the cerebellum is not able to control individual motor unit stimulation rates and recruitment patterns but only whole muscles by means of a single compound signal. The latter signal is identified as the “common drive.” Presumably at the spinal level, special neural circuits are used to decompose the common drive signal into motor unit recruitment patterns and stimulation rates that are specific for a given mode of contraction and probably obey certain optimality principles. Received: 26 February 1999 / Accepted in revised form: 11 June 1999     

A procedure to estimate normal and friction contact parameters in the stance phase of the human gait

A new approach to estimate normal and tangential contact parameters in the foot-ground contact during human gait was proposed. A correct estimation of the contact parameters would be very important in the resolution of predictive forward dynamic problems. The normal contact forces have been well estimated in the literature. But accurate estimation of tangential forces has not been reached yet. This work proposed a new procedure to accurately estimate friction forces. The approach has been based on the consideration of the modulus of the tangential force instead of its components. This modulus was introduced together with the modulus of the normal contact force and its two associated moments in an optimization algorithm to fit the contact forces provided by the model to the experimental data obtained with a force plate. An inverse dynamics problem was solved as a step previous to the optimization algorithm. The results showed that both the normal and tangential forces and the moments in the horizontal plane were in agreement with the experimental measurements. This work also analyzed the influence on the results of the friction law. The results obtained with the general friction law, which considered dry (static and dynamic) and viscous friction, were compared with results provided by simpler laws. The analysis of the components of the friction forces pointed out the importance of the Stribeck component in the resultant force instead of the viscous friction which played a minimal role. But for modelling the stick-slip transition, the implementation of a general friction law is necessary.     

Muscular synergism--I. On criteria for load sharing between synergistic muscles

The use of optimization techniques to predict individual muscle forces in redundant biomechanical systems implies the formulation of a criterion for load sharing between the muscles. In part I of this paper, the characteristics and performance of several linear and non-linear criteria reported in the literature have been compared for static-isometric knee flexion. The results show that linear criteria inherently predict discrete muscle action (orderly recruitment of muscles) whereas non-linear criteria can predict synergistic action. All criteria predict that relatively more force is allocated to muscles with large moment arms. When muscle stresses (or ratios of muscle force to maximum muscle force) are used as the decision variables in the objective function, then relatively more force is allocated to muscles with large maximum possible force as well. Future formulations of the optimization should consider the differences in fiber type composition among the muscles. Such an approach is presented in part II of the paper.     

A method for solution of the multi-objective inverse problems under uncertainty

We describe a method to solve multi-objective inverse problems under uncertainty. The method was tested on non-linear models of dynamic series and population dynamics, as well as on the spatiotemporal model of gene expression in terms of non-linear differential equations. We consider how to identify model parameters when experimental data contain additive noise and measurements are performed in discrete time points. We formulate the multi-objective problem of optimization under uncertainty. In addition to a criterion of least squares difference we applied a criterion which is based on the integral of trajectories of the system spatiotemporal dynamics, as well as a heuristic criterion CHAOS based on the decision tree method. The optimization problem is formulated using a fuzzy statement and is constrained by penalty functions based on the normalized membership functions of a fuzzy set of model solutions. This allows us to reconstruct the expression pattern of hairy gene in Drosophila even-skipped mutants that is in good agreement with experimental data. The reproducibility of obtained results is confirmed by solution of inverse problems using different global optimization methods with heuristic strategies.     

Exploring novel objective functions for simulating muscle coactivation in the neck

Musculoskeletal modeling allows for analysis of individual muscles in various situations. However, current techniques to realistically simulate muscle response when significant amounts of intentional coactivation is required are inadequate. This would include stiffening the neck or spine through muscle coactivation in preparation for perturbations or impacts. Muscle coactivation has been modeled previously in the neck and spine using optimization techniques that seek to maximize the joint stiffness by maximizing total muscle activation or muscle force. These approaches have not sought to replicate human response, but rather to explore the possible effects of active muscle. Coactivation remains a challenging feature to include in musculoskeletal models, and may be improved by extracting optimization objective functions from experimental data. However, the components of such an objective function must be known before fitting to experimental data. This study explores the effect of components in several objective functions, in order to recommend components to be used for fitting to experimental data. Four novel approaches to modeling coactivation through optimization techniques are presented, two of which produce greater levels of stiffness than previous techniques. Simulations were performed using OpenSim and MATLAB cooperatively. Results show that maximizing the moment generated by a particular muscle appears analogous to maximizing joint stiffness. The approach of optimizing for maximum moment generated by individual muscles may be a good candidate for developing objective functions that accurately simulate muscle coactivation in complex joints. This new approach will be the focus of future studies with human subjects.     

Sensitivity of muscle force estimations to changes in muscle input parameters using nonlinear optimization approaches.

The purpose of this study was to analyze the sensitivity of muscle force calculations to changes in muscle input parameters. Force sharing between two synergistic muscles was derived analytically for a one-degree-of-freedom system using three nonlinear optimization approaches. Changes in input parameters that are within normal anatomical variations often caused changes in muscular forces exceeding 100 percent. These results indicate that errors in muscle force calculations may depend as much on inadequate muscle input parameters as they may on the choice of the objective and constraint functions of the optimization approach.     

Non-linear global optimization via parameterization and inverse function approximation: an artificial neural networks approach

In this article, a novel technique for non-linear global optimization is presented. The main goal is to find the optimal global solution of non-linear problems avoiding sub-optimal local solutions or inflection points. The proposed technique is based on a two steps concept: properly keep decreasing the value of the objective function, and calculating the corresponding independent variables by approximating its inverse function. The decreasing process can continue even after reaching local minima and, in general, the algorithm stops when converging to solutions near the global minimum. The implementation of the proposed technique by conventional numerical methods may require a considerable computational effort on the approximation of the inverse function. Thus, here a novel Artificial Neural Network (ANN) approach is implemented to reduce the computational requirements of the proposed optimization technique. This approach is successfully tested on some highly non-linear functions possessing several local minima. The results obtained demonstrate that the proposed approach compares favorably over some current conventional numerical (Matlab functions) methods, and other non-conventional (Evolutionary Algorithms, Simulated Annealing) optimization methods.     

Particle swarm optimizer for weighting factor selection in intensity-modulated radiation therapy optimization algorithms

In inverse treatment planning of intensity-modulated radiation therapy (IMRT), the objective function is typically the sum of the weighted sub-scores, where the weights indicate the importance of the sub-scores. To obtain a high-quality treatment plan, the planner manually adjusts the objective weights using a trial-and-error procedure until an acceptable plan is reached. In this work, a new particle swarm optimization (PSO) method which can adjust the weighting factors automatically was investigated to overcome the requirement of manual adjustment, thereby reducing the workload of the human planner and contributing to the development of a fully automated planning process. The proposed optimization method consists of three steps. (i) First, a swarm of weighting factors (i.e., particles) is initialized randomly in the search space, where each particle corresponds to a global objective function. (ii) Then, a plan optimization solver is employed to obtain the optimal solution for each particle, and the values of the evaluation functions used to determine the particle’s location and the population global location for the PSO are calculated based on these results. (iii) Next, the weighting factors are updated based on the particle’s location and the population global location. Step (ii) is performed alternately with step (iii) until the termination condition is reached. In this method, the evaluation function is a combination of several key points on the dose volume histograms. Furthermore, a perturbation strategy – the crossover and mutation operator hybrid approach – is employed to enhance the population diversity, and two arguments are applied to the evaluation function to improve the flexibility of the algorithm. In this study, the proposed method was used to develop IMRT treatment plans involving five unequally spaced 6 MV photon beams for 10 prostate cancer cases. The proposed optimization algorithm yielded high-quality plans for all of the cases, without human planner intervention. A comparison of the results with the optimized solution obtained using a similar optimization model but with human planner intervention revealed that the proposed algorithm produced optimized plans superior to that developed using the manual plan. The proposed algorithm can generate admissible solutions within reasonable computational times and can be used to develop fully automated IMRT treatment planning methods, thus reducing human planners’ workloads during iterative processes.     

Dynamics of granzyme B-induced apoptosis: mathematical modeling

Apoptosis is mediated by an intracellular biochemical system that mainly includes proteins (procaspases, caspases, inhibitors, Bcl-2 protein family as well as substances released from mitochondrial intermembrane space). The dynamics of caspase activation and target cleavage in apoptosis induced by granzyme B in a single K562 cell was studied using a mathematical model of the dynamics of granzyme B-induced apoptosis developed in this work. Also the first application of optimization approach to determination of unknown kinetic constants of biochemical apoptotic reactions was presented. The optimization approach involves solving of two problems: direct and inverse. Solving the direct optimization problem, we obtain the initial (baseline) concentrations of procaspases for known kinetic constants through conditional minimization of a cost function based on the principle of minimum protein consumption by the apoptosis system. The inverse optimization problem is aimed at determination of unknown kinetic constants of apoptotic biochemical reactions proceeding from the condition that the optimal concentrations of procaspases resulting from the solution of the direct optimization problem coincide with the observed ones, that is, those determined by biochemical methods. The Multidimensional Index Method was used to perform numerical solution of the inverse optimization problem.     

Velocity-dependent cost function for the prediction of force sharing among synergistic muscles in a one degree of freedom model

Prediction of accurate and meaningful force sharing among synergistic muscles is a major problem in biomechanics research. Given a resultant joint moment, a unique set of muscle forces can be obtained from this mathematically redundant system using nonlinear optimization. The classical cost functions for optimization involve a normalization of the muscle forces to the absolute force capacity of the target muscles, usually by the cross-sectional area or the maximal isometric force. In a one degree of freedom model this leads to a functional relationship between moment arms and the predicted muscle forces, such that for constant moment arms, or constant ratios of moment arms, agonistic muscle forces increase or decrease in unison. Experimental studies have shown however that the relationship between muscle forces is highly task-dependent often causing forces to increase in one muscle while decreasing in a functional agonist, likely because of the contractile conditions and contractile properties of the involved muscles. We therefore, suggest a modified cost function that accounts for the instantaneous contraction velocity of the muscles and its effect on the instantaneous maximal force. With this novel objective function, a task-dependent prediction of muscle force distribution is obtained that allows, even in a one degree of freedom system, the prediction of force sharing loops, and simultaneously increasing and decreasing forces for agonist pairs of muscles.     

Are natural coordinates a useful tool in modeling planar biomechanical linkages?

The paper analyzes the use of natural coordinates in modeling and simulation of musculoskeletal models of the human body. A biomechanical model of the lower extremity of the human body was constructed for the analysis. It consists of three anatomical segments described by eight natural coordinates located at the joints. The developed model was applied to solve three classic dynamics problems of human motion: inverse dynamics, direct dynamics, and static optimization. The analysis covers the raising of a leg together with; the time characteristics of the resultant net torques at the basic joints of the leg (inverse dynamics), the time histories of natural coordinates and their velocities (direct dynamics) as well as the time-varying muscle force patterns (static optimization). In order to check the numerical efficiency of modeling in the natural coordinates' environment the results were compared with the ones received through generalized coordinates. Some conclusions drawn from this comparison and final remarks referring to the biomechanical identification of the analyzed motor task were included in the paper.     

A new motor model representing the stretch-induced force enhancement and shortening-induced force depression in skeletal muscle

A motor model that consists of two Maxwell elements with a force generator and one Voigt element is proposed in this paper. The motor model can achieve a hyperbolic force-velocity relation when we alter weight functions applied to the Maxwell elements and the force generator. Rate coefficients are introduced to determine the weight function and to improve the motor performance and the time course of the motor force. The weight functions are used as a controller of the motor. We assume that the mechanical impulse applied to the motor affects the rate coefficients and found that the amount of the mechanical impulse is related to the amount of force depression following motor shortening and to the amount of force enhancement following motor stretching. The time courses of the motor force following shortening and stretching quantitatively resemble those in other muscle experiments. The maximum energy efficiency of the motor that we obtained was 50% with an ATP hydrolysis type and 25% with an AC-DC motor type.     

A neuromusculoskeletal tracking method for estimating individual muscle forces in human movement

A neuromusculoskeletal tracking (NMT) method was developed to estimate muscle forces from observed motion data. The NMT method combines skeletal motion tracking and optimal neuromuscular tracking to produce forward simulations of human movement quickly and accurately. The skeletal motion tracker calculates the joint torques needed to actuate a skeletal model and track observed segment angles and ground forces in a forward simulation of the motor task. The optimal neuromuscular tracker resolves the muscle redundancy problem dynamically and finds the muscle excitations (and muscle forces) needed to produce the joint torques calculated by the skeletal motion tracker. To evaluate the accuracy of the NMT method, kinematics and ground forces obtained from an optimal control (parameter optimization) solution for maximum-height jumping were contaminated with both random and systematic noise. These data served as input observations to the NMT method as well as an inverse dynamics analysis. The NMT solution was compared to the input observations, the original optimal solution, and a simulation driven by the inverse dynamics torques. The results show that, in contrast to inverse dynamics, the NMT method is able to produce an accurate forward simulation consistent with the optimal control solution. The NMT method also requires 3 orders-of-magnitude less CPU time than parameter optimization. The speed and accuracy of the NMT method make it a promising new tool for estimating muscle forces using experimentally obtained kinematics and ground force data.     

Impact of Stoichiometry Representation on Simulation of Genotype-Phenotype Relationships in Metabolic Networks

Genome-scale metabolic networks provide a comprehensive structural framework for modeling genotype-phenotype relationships through flux simulations. The solution space for the metabolic flux state of the cell is typically very large and optimization-based approaches are often necessary for predicting the active metabolic state under specific environmental conditions. The objective function to be used in such optimization algorithms is directly linked with the biological hypothesis underlying the model and therefore it is one of the most relevant parameters for successful modeling. Although linear combination of selected fluxes is widely used for formulating metabolic objective functions, we show that the resulting optimization problem is sensitive towards stoichiometry representation of the metabolic network. This undesirable sensitivity leads to different simulation results when using numerically different but biochemically equivalent stoichiometry representations and thereby makes biological interpretation intrinsically subjective and ambiguous. We hereby propose a new method, Minimization of Metabolites Balance (MiMBl), which decouples the artifacts of stoichiometry representation from the formulation of the desired objective functions, by casting objective functions using metabolite turnovers rather than fluxes. By simulating perturbed metabolic networks, we demonstrate that the use of stoichiometry representation independent algorithms is fundamental for unambiguously linking modeling results with biological interpretation. For example, MiMBl allowed us to expand the scope of metabolic modeling in elucidating the mechanistic basis of several genetic interactions in Saccharomyces cerevisiae.