|Home | About | Journals | Submit | Contact Us | Français|
Finite element (FE) modeling and multibody dynamics have traditionally been applied separately to the domains of tissue mechanics and musculoskeletal movements, respectively. Simultaneous simulation of both domains is needed when interactions between tissue and movement are of interest, but this has remained largely impractical due to high computational cost.
Here we present a method for concurrent simulation of tissue and movement, in which state of the art methods are used in each domain, and communication occurs via a surrogate modeling system based on locally weighted regression. The surrogate model only performs FE simulations when regression from previous results is not within a user-specified tolerance. For proof of concept and to illustrate feasibility, the methods were demonstrated on an optimization of jumping movement using a planar musculoskeletal model coupled to a FE model of the foot. To test the relative accuracy of the surrogate model outputs against those of the FE model, a single forward dynamics simulation was performed with FE calls at every integration step and compared with a corresponding simulation with the surrogate model included. Neural excitations obtained from the jump height optimization were used for this purpose and root mean square (RMS) difference between surrogate and FE model outputs (ankle force and moment, peak contact pressure and peak von Mises stress) were calculated.
Optimization of jump height required 1800 iterations of the movement simulation, each requiring thousands of time steps. The surrogate modeling system only used the FE model in 5% of time steps, i.e. a 95% reduction of computation time. Errors introduced by the surrogate model were less than 1 mm in jump height and RMS errors of less than 2 N in ground reaction force, 0.25 Nm in ankle moment, and 10 kPa in peak tissue stress.
Adaptive surrogate modeling based on local regression allows efficient concurrent simulations of tissue mechanics and musculoskeletal movement.
Computational biomechanics has largely been separated into two distinct modeling domains, finite element analysis (FEA) [e.g. 1, 2] and multi-body dynamics. Due to computational efficiency, muscle driven multibody models have been the primary method used in optimization of movement patterns . While predicting resultant joint loads and muscle forces, musculoskeletal models generally do not provide detailed representation of soft-tissue structures. Therefore, the distribution of muscle forces and joint loads at tissue levels and effects of tissue properties on human movement cannot be studied. Conversely, studies focusing on soft-tissue structures have historically utilized finite element (FE) methods that required significant computational resources and well-defined boundary conditions [4, 5]. From analyzing a specific structure such as a medial collateral ligament (MCL) in the knee to modeling a whole joint or organ such as the foot, FE methods have the ability to yield important soft-tissue information [6–10] not found in musculoskeletal simulations. There is, however, currently no method for incorporating mechanical or sensory effects of soft tissue deformations into predictive musculoskeletal simulations. Creating a framework that spans both domains would allow simulations of this coupled behavior of muscle actuated multi-body dynamics with realistic soft-tissue models.
Combining the benefits of two domains, especially for use in an optimization scheme (usually required for predictive movement simulations), is a methodological and computational challenge. At a similar scale, development of multi-domain analyses incorporating fluid-solid interactions and structural analysis techniques for automotive crash analysis, aerospace applications, and fatigue have illustrated the possibility of multi-domain simulations [11–15]. In musculoskeletal biomechanics, attempts have been made to apply multi-domain techniques but these usually consisted of non-concurrent simulations. Typically, soft-tissue FE models were driven with boundary conditions supplied by a musculoskeletal simulation and effectively served as a post-processing tool . This does not allow prediction of how tissue may affect skeleton movement, either through mechanics (e.g. joint laxity) or through neural pathways (e.g. osteoarthritic pain).
Of notable exception, Koolstra and van Eijden  attempted concurrent simulations of the temporomandibular joint and jaw structure using muscle activations, rigid-body dynamics and soft-tissue deformation. An explicit framework was utilized and the computational cost for each solution was not reported. A major challenge in concurrent domain coupling is that FE simulations are required at each time step of a movement simulation. Typically, a forward dynamic simulation of movement contains thousands of time steps, and an iterative movement optimization may require thousands of such movement simulations until the optimal movement is found. This adds up to millions of FE simulations, which would thus require enormous computational resources in order to solve just one movement optimization problem. In order to obtain solutions, modelers typically focus on one of the modeling domains, while simplifying the other. For instance, surface-surface penetration has been used within multibody dynamics to compute reaction loads in the knee  or between foot and ground . This is a good approximation when tissue deformation is limited to a surface layer but not generally applicable.
Under conditions where the analysis requires iterative simulations of a computationally expensive model, surrogate modeling is often employed. In general, surrogate modeling approaches can be classified as global or local methods. Global methods fit a statistical regression model to a defined set of input/output sets. Accuracy of a global method depends on the number of available data sets and the goodness of fit of the approximation over the whole domain. Examples include response surface techniques and neural networks. Lin et al.  developed a response surface approximation of knee joint contact mechanics and demonstrated its feasibility for potential use in optimization routines. This promising work showed a significant reduction in computational cost associated with the use of the surrogate model but requires an a priori estimate of input data ranges for response surface fitting. In addition, for higher dimensional input spaces, response surface approximations of complicated or highly nonlinear behavior are difficult to capture with a low-order polynomial or other function approximators. User input would also be required to produce a new approximating function whenever the underlying model is changed or updated, such as for patient specific models of joint contact or soft-tissue restraint. Local methods use a set of neighboring points only and include locally weighted regression, spline fitting or radial basis functions. Lazy learning  is one form of locally weighted regression based on linear or polynomial fits to neighboring points. It is particularly attractive because it retains all the original data and can provide error estimates to drive an adaptive sampling scheme for generating additional data. This allows unimportant areas of the domain space to be avoided and the highly nonlinear areas can be densely sampled to accurately describe the response.
The objective of this study was to illustrate that finite element analysis of tissue deformations can be coupled to musculoskeletal movement simulations concurrently and effectively by the use of an adaptive surrogate modeling scheme. To realize this possibility and assess feasibility, an optimal control solution for maximum height jumping was obtained using a musculoskeletal model of the lower extremity, a finite element model of the foot and a corresponding adaptive surrogate model representation of the finite element model. We specifically explored answers to the following questions: (1) Is it possible to perform a forward dynamic movement optimization using this system? (2) How comparable is the movement simulation when using the surrogate foot model against that fully coupled with the FE model? (3) How much computation time is required when using the surrogate foot model?
The musculoskeletal model has been described previously [22, 23]. The model contained seven body segments: trunk, thighs, shanks, and feet. Joints were assumed to be ideal hinges, and there were no kinematic constraints between the feet and ground, resulting in a total of nine kinematic degrees of freedom. Eight muscle groups were included in each lower extremity: Iliopsoas, Glutei, Hamstrings, Rectus Femoris, Vasti, Gastrocnemius, Soleus, and Tibialis Anterior. Each muscle was represented by a 3-element Hill model, as described in McLean et al. , with muscle properties from Gerritsen et al. , and simulated with custom C code. This model has 50 state variables: 9 generalized coordinates, 9 generalized velocities, 16 muscle contractile element lengths, and 16 muscle activations. Equations of motion were generated by SD/Fast (Parametric Technology Corp., Needham, MA):
where M is a mass matrix, C are centrifugal and coriolis effects, G represents gravity, and FMT are the muscle forces, applied via a matrix of moment arms R. The final term represents reaction loads applied to the foot segment by the finite element model of the foot, which will be introduced below. In the absence of friction and viscoelastic effects, these loads are only dependent on kinematic boundary conditions qFEA which are a known function of skeleton pose q.
The model was configured in an initial squat position, where the joint angles were chosen to prevent excessive passive force contribution from extensor muscles. The vertical position and trunk orientation were then calculated in order to satisfy static equilibrium conditions (Figures 1 and and2).2). An optimal steady state, which minimized neural excitation values while maintaining zero accelerations, provided state variables of the muscles (activation and muscle length) that will be used as the initial condition (along with rigid body degrees of freedom) for all forward dynamic simulations.
A plane strain foot model (Figures 1 and and2)2) was implemented in ABAQUS (Simulia, Providence, RI). A sagittal plane cross-section along the second ray of the foot was used to represent the bone and tissue geometry. Out-of-plane thickness was set to an approximate adult foot width 80 mm. Selection of a thickness value allows adequate representation of ground reaction forces from predicted contact pressures. Bones were modeled as rigid and the soft tissue was assumed to be a nonlinearly elastic (hyperelastic) incompressible material. More specifically, coefficients of an Ogden material model, which minimized the differences between model predicted and experimental response of the heel pad under indentation were used . Bones other than the phalanges were combined into one rigid segment, which was controlled by prescribing the vertical position and the orientation of the talus relative to the ground. These were the kinematic boundary conditions needed to run finite element simulations. In effect, the FE model of the foot and the musculoskeletal model were directly coupled by sharing boundary conditions at a point in the talocrural joint and thus, the ankle is modeled as a hinge joint. The phalanges were represented as another rigid segment, which was free to move during simulations. Soft tissue surrounding the metatarsophalangeal joint served to restrain the movements of this segment during passive toe flexion. Elements between the metatarsal head and the proximal phalanx also contributed to passive MTP joint stiffness and were modeled as linearly elastic (E=1e6 Pa, v=0.3). Contact between the plantar aspect of the foot and the ground was modeled as frictionless, eliminating the need to prescribe the horizontal translation of the talus during simulations. Once the vertical translation of the talus and its orientation was passed to the finite element model, the FE simulations were capable of returning the vertical reaction force and moment at the talus to the musculoskeletal model. Stress-strain distribution within the soft tissue and plantar pressures were byproducts of the finite element analysis that can be used to control movement in future studies.
The 2D finite element model was developed to align with the neutral position of the foot in the musculoskeletal model. Ankle joint coordinates, qFEA, were directly coupled between the FE and musculoskeletal models. Coupled time-domain boundary conditions such as acceleration and velocity were not necessary as the FE model did not include mass, inertial effects or time-dependent material properties.
The Lazy Learning Toolbox [26, 21] for Matlab (Mathworks, Inc., Natick, MA) was used as the surrogate modeling tool with two inputs, ankle vertical position and plantar/dorsiflexion rotation, and four outputs, vertical load, plantar/dorsiflexion moment, peak plantar pressure (PPP), and peak von Mises stress. Lazy Learning is a local interpolation method based on the use of nearest neighbor input/output sets present in the database. The linear regression option was utilized and a leave-one-out cross validation error (CVE) was computed, based on a regression model using the nearest N neighbors. The distance from the query point Xq to candidate neighbors Xi was defined as the “Manhattan” distance:
where Xij and Xqj are the jth coordinates of the database point and query point, respectively. In the present application, m=2 is the number of dimensions in the input space and wj are the weighting factors that define its distance metric. In the present application, weights were set to 1.0 and inputs were normalized on the database range which has units of meters (translation) and radians (orientation). The number of nearest neighbors was allowed to range from 7 to 20 and for each of the four model outputs, the number of neighbors was selected based on minimization of CVE.
Cross validation errors of the local regression model were compared to prescribed error tolerances for reaction force and reaction moment. Initially tight tolerances were used to populate the database. Thereafter, the tolerances were set to 200 N and 0.2 Nm, respectively, based on assessment of surrogate model outputs against FE model results using this initial database. When both CVE estimates were below the specified tolerances, the local linear regression model was used to predict output. When either error was above the specified tolerance, an FE simulation was completed and the results were supplied to both the musculoskeletal model and the database (Figure 1). A complete description of the lazy learning algorithms can be found in Atkeson et al. .
To test the efficacy of the multi-domain simulation, an optimization was performed to generate a maximal height jumping movement. Left-right symmetrical neural excitation patterns for the eight muscles were parameterized as 32 parameters along simulation time: the excitation values for 4 time points of 0.09, 0.18, 0.225, and 0.27 seconds. Time values were chosen based on our preliminary jumping simulation studies and included a neural excitation parameter near the expected toe-off (0.225 seconds). It is possible to select a larger number of nodes to identify a finer jumping control scheme but it is not necessary to illustrate the concurrent simulation framework. To start the optimization control variable neural excitations were set to 0.5. Bounds were prescribed on the neural excitations to only allow a range of 0.01 to 1.0. A lower bound of 0.01 was specified to avoid an inherently unstable condition if the muscles were to impart zero force. It should also be noted that once the tolerance values were chosen, the initial database was cleared. This allowed the true contribution of the surrogate model to be assessed over the optimization routine.
Each objective function evaluation consisted of one complete forward simulation using a parameter vector p containing the 32 neural excitation parameters. The forward simulation was terminated at the beginning of the flight phase at which time the objective function was calculated as the center of mass jump height:
where y and vy are the vertical position and velocity of the center of mass, respectively, when the feet leave the ground, and g is the gravitational acceleration. The bounded maximization problem was solved using sequential quadratic programming (Matlab Optimization Toolbox, Mathworks, Inc., Natick, MA). The convergence criteria for the objective function was set to 1/10th of a millimeter.
As a verification of the surrogate model, neural excitation values from the optimized jump were utilized to compare results from a directly coupled musculoskeletal and FE model simulation and a corresponding simulation with the surrogate model included. FE results were utilized at every integration step for the directly coupled simulation. Root mean square (RMS) errors were calculated between the two simulations to compare the objective function (jump height) and the foot model outputs (reaction loads and tissue stress).
The optimization routine successfully achieved a center of mass jump height of 18.2 cm with respect to a standing configuration (Figures 2 and and3).3). During the optimization, approximately 1800 movement simulations were performed, each requiring several thousand evaluations of the foot model. A total of 51 optimization iterations were performed, with each iteration consisting of an initial model evaluation followed by 32 forward simulations for gradient calculation, and a few more simulations to find the minimum for the iteration. Average percentage FEA runs were calculated over all 51 iterations for initial function evaluation and calculation of individual components of the gradient (Figure 4). As an iteration proceeded from function evaluation to gradient calculation, the number of FE simulations decreased from 13% to 2% (Figure 4). This demonstrates that for the relatively small changes to the control variables during the gradient calculations, the surrogate model was effective in learning a specific area of the input space. Over the whole optimization, utilization of the surrogate model resulted in an average reduction of 95% in the number of potential FE simulations, compared to direct coupling between the FE model and musculoskeletal model (Figure 4). This reduction allowed the movement prediction simulations to complete in approximately 4 weeks on a Linux based dual processor Intel® Xeon 3.4 GHz machine with 6GB of memory.
The final database contained over 140,000 FE input/output sets. Each FE simulation required from 4 to 50 seconds to converge, depending on foot position and orientation. Color graded input/output sets, based on a calculated slope using the 10 closest neighbors, highlighted the nonlinear nature of the FE foot model (Figure 5). As this is a maximal effort simulation many data points were required in high load, and thus very stiff, areas of the database. These points tended to be added in the later stages of the movement optimization.
Accuracy of the surrogate model simulation was within 1 mm, obtained from the difference in predicted jump height versus that obtained using FE simulations at every integration step. Optimal neural excitation values were utilized during one forward simulation for the comparison. Surrogate model predicted ankle reaction loads were acceptable when compared to the FE results throughout the maximal height jumping simulation with RMS errors of 1.59 N for the vertical load and 0.244 Nm for the plantar/dorsiflexion moment (Figure 6). Peak contact pressure and von Mises stresses were also predicted throughout the jumping simulations and showed very good agreement between surrogate model predicted values and FE results with RMS errors of 2 and 7 kPa (0.7% and 1.0% of maximum value predicted during jumping), respectively.
The presented modeling methodology successfully optimized jump height, using neural excitation patterns as control variables, with a coupled musculoskeletal lower extremity model and FE model of the foot. The FE model results were stored in a database, from which a surrogate model attempted to predict subsequent FE results using local regression. When the estimated error of the regression model was below the specified tolerance, the surrogate model output was used, otherwise a new FE simulation was performed. As expected, this surrogate modeling scheme was able to gradually eliminate the need for FE simulations, thus removing its computational cost that appears to be a bottleneck in concurrent simulations of musculoskeletal movements and tissue deformations. Furthermore, we have shown that the optimized jumping simulation using the surrogate foot model sufficiently reproduced the same neuromuscular movement simulation directly coupled to the actual FE model. Finally, we demonstrated that the surrogate modeling method provided good estimates of tissue level variables, such as peak stress.
Maximal height jumping has been extensively studied in the literature and was also chosen for this study because of its straightforward objective to be used as a test problem. Utilizing two-dimensional models, Pandy et al.  achieved a jump height of 33 cm and van Soest et al.  a height of 39.2 cm, with jump height defined by the vertical displacement of the center of mass relative to the standing configuration. Comparing to the jump height of 18.2 cm in this study, the difference in performance may be partly due to differences in muscle properties, including generally lower maximum force producing capability for individual muscles of our study. Additional differences include choice of control variables and the assumption of frictionless contact, which has been shown to reduce predicted jump height . Passive toe flexion also likely reduced the maximum achievable jump height but we do not exclude the possibility that the gradient based SQP algorithm did not reach to a global optimum. We should therefore consider this solution to be a local optimum and future work will explore implementation using a global optimization routine. Nevertheless, this movement optimization served as a good vehicle to demonstrate the feasibility of utilizing a multi-domain simulation in a computationally intensive optimization of movement.
Further advances to the presented model will include incorporation of friction along with validated 2D and 3D FE foot models. Interface loads and peak plantar pressures are influenced by friction and the path taken to a given foot-ground orientation will affect the deformation of the plantar tissue. Simulated jumping does not necessitate friction whereas other movement patterns such as gait require the shear force supplied by the foot-ground interaction. When friction is included, path dependent kinematic variable(s) will need to be incorporated into the estimation for accurate approximations using a surrogate model. As one might expect, the dimension of the input space, and the number of data points, will grow substantially as these features are included. The local regression approach with adaptive sampling has the potential to avoid the “curse of dimensionality” by only generating database points where needed. In practical applications, database management can be costly and more sophisticated neighbor searching methods  will be useful as the complexity of the model increases. Even in this study, with the final database size of over 140,000 points (Figure 5), database management contributed significantly to the overall computational time.
Computational benefits of the surrogate model were assessed based on literature reported run times and in-house simulations. Computation times for the optimization were not reported by Pandy et al.  for two-dimensional jumping simulation. When the model was further developed in three dimensions, the optimization routine for maximal jumping required between 1 and 2.5 months (for the single processor machines) . As more complex models and movement patterns are adopted the computational expense further increases with one study citing 10,000 hours of computational time for a gait cycle optimization . Mclean et al.  reported 37 hours of computational time to simulate a cutting maneuver and Neptune et al  required 660 hours to optimize a simulation of running. We are aware that all these simulations considered different number of muscles, nodal parameters for muscle excitations and were conducted using various computational platforms. Nevertheless, it is clear that movement prediction takes considerable computational time even when one does not consider soft tissue deformations through coupling. None of the studies cited above attempted concurrent simulations. As stated earlier, the complete optimization routine for this study, even using a surrogate model of foot deformations, required approximately 4 weeks of computational time (672 hours). Computational expense for the forward simulations during the optimization routine varied dramatically, from ~1 minute up to multiple hours, depending on the percentage of FE simulations and the database size. It would not have been feasible to perform a movement optimization without the surrogate model. One forward simulation using neural excitations from the optimized jump with direct FE analysis at every integration point required 9.3 hours of computational time. If direct coupling had been implemented for the movement prediction problem, FE analyses would be performed at every integration point during the ~1800 function calls (each representing one forward simulation). The associated expense for the complete optimization routine in this study would have required approximately 698 days (16,752 hours) of computational time. Obviously, computational expense is a central consideration when performing optimizations of movement patterns and becomes even more important when soft-tissue deformation is included.
The performance of the locally weighted regression method will depend on algorithm parameters, such as polynomial order, the number of neighbors considered, distance metric in input space, and the tolerance value for adaptive sampling. For this study we selected a set of values that provided reasonable local regression characteristics but an extensive parameter tuning is warranted to decrease requested FE simulations without diminishing regression accuracy. Based on the errors found in the surrogate model after completion of the movement optimization, we suspect that the cross validation error estimates are overly pessimistic and that fewer FE simulations, possibly by relaxing the tolerances, would have been sufficient to achieve good results. While the 0.244 Nm RMS error for the ankle moment output metric appeared to be relatively high, peak errors occurred at the high load areas of the database and still represented less than one percent of the applied moment. Sensitivity of the jump height and accuracy of the surrogate model to changes in the interpolation parameters remain a future direction of this work. Regardless, the linear approximation method proved to be accurate, and through successive refinement in the database, it reduced the potential number of FE simulations during an optimization iteration by 95% (Figure 4). As a result, the disproportionate computational cost associated with the FE model was overcome, while the coupled behavior of the musculoskeletal and tissue models was retained. Implementing a higher order regression technique could lessen the computational FE burden, and thus the number of database points, but may require more time to perform each regression.
Exploration of the ankle-foot complex has clinical applications in the prevention of diabetic foot ulceration. With the coupled simulations, we will be able to explore the closed loop interactions between sensory loss, neuromuscular control, and tissue stress and damage. The proposed methodology, however, is not limited to this specific case. Any coupled, computationally expensive modeling system could potentially benefit from a surrogate modeling approach. The complex behavior of the knee would be a very good application where soft-tissue effects and joint level mechanics could be predicted. Traditionally, computational models of the knee have required substantial resources and boundary conditions that have not included musculoskeletal loading. Models of shoulder, hip and other joints of interest could also be developed. Of particular clinical interest, the defined geometry and material behavior of joint replacements would be well-suited to this method. Other applications in biomechanics could include tissue-fluid interactions and coupling of cellular mechanics to tissue and organ level models.
Direct coupling of finite element analysis to a single forward dynamics of a musculoskeletal model has been shown to be possible . However, predictive simulations that require multiple solutions of the forward dynamics problem can only be possible with a cost-effective approach. To the author’s knowledge, this is the first study to complete a predictive optimization of an active movement with a coupled musculoskeletal and FE model of tissue level mechanics. A surrogate modeling technique was developed to efficiently and adaptively predict the joint reaction loads and important soft-tissue conditions of a corresponding FE model. Far less simplification of joint behavior versus traditional musculoskeletal modeling is an important benefit of this method, and the ability to utilize and predict tissue and joint mechanics adds clinical insight. This optimized muscle-loaded simulation helps to further advance the state of musculoskeletal modeling and is an important step toward development of musculoskeletal simulation strategies more aware of tissue deformations.
This study is supported by the National Institutes of Health grant 1 R01 EB006735-01. The authors would like to thank Scott Sibole for developing the finite element model of the foot and Anna Fernandez and Marko Ackermann for initial work on the jump optimizations.