The Neuromusculoskeletal Modeling Pipeline: MATLAB-based model personalization and treatment optimization functionality for OpenSim

Abstract

Neuromusculoskeletal injuries including osteoarthritis, stroke, spinal cord injury, and traumatic brain injury affect roughly 19% of the U.S. adult population. Standardized interventions have produced suboptimal functional outcomes due to the unique treatment needs of each patient. Strides have been made to utilize computational models to develop personalized treatments, but researchers and clinicians have yet to cross the “valley of death” between fundamental research and clinical usefulness. This article introduces the Neuromusculoskeletal Modeling (NMSM) Pipeline, two MATLAB-based toolsets that add Model Personalization and Treatment Optimization functionality to OpenSim. The two toolsets facilitate computational design of individualized treatments for neuromusculoskeletal impairments through the use of personalized neuromusculoskeletal models and predictive simulation. The Model Personalization toolset contains four tools for personalizing 1) joint structure models, 2) muscle–tendon models, 3) neural control models, and 4) foot–ground contact models. The Treatment Optimization toolset contains three tools for predicting and optimizing a patient’s functional outcome for different treatment options using a patient’s personalized neuromusculoskeletal model and direct collocation optimal control methods. Support for user-defined cost, constraint, and model modification functions facilitate simulation of a vast number of possible treatments. An NMSM Pipeline use case is presented for an individual post-stroke with impaired walking function, where the goal was to predict how the subject’s neural control could be changed to improve walking speed without increasing metabolic cost. First the Model Personalization toolset was used to develop a personalized neuromusculoskeletal model of the subject starting from a generic OpenSim full-body model and experimental walking data (video motion capture, ground reaction, and electromyography) collected from the subject at his self-selected speed. Next the Treatment Optimization toolset was used with the personalized model to predict how the subject could recruit existing muscle synergies more effectively to reduce muscle activation disparities between the paretic and non-paretic legs. The software predicted that the subject could increase his walking speed by 60% without increasing his metabolic cost per unit time by modifying existing muscle synergy recruitment. This hypothetical treatment demonstrates how NMSM Pipeline tools could allow researchers working collaboratively with clinicians to develop personalized neuromusculoskeletal models of individual patients and to perform predictive simulations for designing personalized treatments that maximize a patient’s post-treatment functional outcome.

Introduction

Approximately 19% of the U.S. adult population is afflicted by a movement impairment caused by stroke, osteoarthritis, traumatic brain injury, limb amputation, cerebral palsy, Parkinson’s disease, or spinal cord injury [1,2,3]. In addition to the physical challenges caused by a movement impairment, these individuals are more likely to experience increased health care costs, lower work productivity, a higher risk of developing heart disease or diabetes, a reduction or loss of independence, and a decreased quality of life [1, 4]. At present, treatment of these life-altering impairments typically involves selection of a standardized treatment based on subjective clinical assessment rather than design of a customized treatment based on objective evidence-based methods. Consequently, existing treatments for movement impairments do not address the unique clinical situation and needs of the individual patient to the extent that they could, and should.

Not surprisingly, the standardized treatment paradigm is often unable to achieve the level of post-treatment functional recovery that patients desire and expect. The percentage of patients who are dissatisfied with their functional outcome following a standardized treatment supports this viewpoint. For those who experience stroke, satisfaction with functional outcome is closely tied to ability to return to a pre-stroke level of movement function [5]. However, stroke patient satisfaction with mobility improvements generally declines after hospital discharge [5], with only 65% of individuals post-stroke recovering the ability to walk [6]. Even then, they tend to walk at a slow self-selected speed with an asymmetric gait pattern that is less efficient metabolically than for healthy individuals [6,7,8]. Furthermore, 20% of stroke survivors remain dependent on others to perform activities of daily living (ADLs) [9]. For those with physician-diagnosed arthritis, roughly 44% have functional limitations that affect their ability to perform ADLs [10]. For those who receive a total knee replacement (TKR), the best predictor of post-surgery patient satisfaction is functional outcome [11]. However, 22% of TKR patients are dissatisfied with their functional outcome one year after surgery [12], while 16 to 30% are dissatisfied with their ability to perform specific ADLs [13]. For those who receive a total hip replacement, 7% are dissatisfied with their functional outcome [14], with expectations for improved walking and stair climbing ability not being met in over 30% of patients [15]. For those who receive opening wedge high tibial osteotomy surgery, 22% are dissatisfied with their outcome, primarily due to unmet expectations for their post-surgery ability to perform functional activities [16]. For those with cerebral palsy (CP), 33% do not improve walking speed as a result of multi-level orthopedic surgery [17], while 21% of parents of children with CP do not feel that their child’s orthopedic surgery was worth the time, effort, pain, and cost [18]. More than 50% of children with CP who undergo multi-level surgeries do not see improvements in their gait patterns [19]. For those with lower limb amputation, between 12 and 32% cannot perform desired ADLs with their prosthesis, while only 30% to 52% find the same ADLs very easy to perform [20]. Those with upper limb amputation use their prosthesis primarily for cosmetic rather than functional reasons, with the percentage of individuals who rate their upper extremity prosthesis as beneficial for performing various ADLs ranging from 5% to only 25% [21]. These statistics across a myriad of clinical conditions that impair movement highlight the need for a better treatment design paradigm to help a larger percentage of patients recover their desired level of post-treatment function.

One possible way to improve treatment design and thus functional outcomes for individuals with movement impairments is to employ engineering design optimization methods that utilize physics-based computational models. Physics-based computational models have an advantage over “black box” machine learning models since they extrapolate better to situations outside the boundaries of the training data, which is where optimal solutions often lie. Engineering design optimization methods have dramatically improved function and reliability of commercial products in numerous industries including automotive, aerospace, farming, construction, and structural engineering, among others [22,23,24,25,26,27]. Using these virtual prototyping methods, engineers can evaluate countless designs–including highly novel ones–in a time- and cost-efficient manner to determine an optimal design for a given problem, thereby reducing or eliminating trial-and-error design iterations performed on physical prototypes. Researchers have already demonstrated the benefits of combining engineering design optimization with personalized models to develop effective personalized interventions for a small number of clinical problems. [28,29,30,31,32,33]

While a similar computational approach holds promise for designing personalized treatments for movement impairments, at least three significant challenges must first be overcome for such an approach to become a reality. First, researchers must be able to personalize the neuromusculoskeletal computer models used in the computational treatment design process [34]. Published musculoskeletal modeling studies have typically used scaled generic models rather than models whose relevant characteristics have been personalized to a subject’s movement data [34]. This decision is understandable given the challenges of collecting the experimental movement (and possibly imaging) data needed for the model personalization process coupled with the challenges of performing the personalization process itself for any aspect of a model (e.g., joint functional axes). However, differences between patients with the same clinical condition are precisely what make the clinical treatment design process so challenging. Second, researchers must be able to perform the treatment design optimization process using a patient’s personalized neuromusculoskeletal model. Such a process requires performing repeated predictive simulations of patient function following different implementations of a planned intervention. To date, the majority of such simulations have been performed using “home brew” computational software [35,36,37,38,39,40,41,42,43,44,45,46,47], making it difficult for the broader research community to engage in the process, to repeat studies performed by others, and to transfer knowledge about the process to others. Third, researchers must have access to computational software that makes the model personalization and treatment optimization processes easy, repeatable, and transferable. While recent software advances such as OpenSim Moco [48] have made it possible to perform predictive simulations more easily, neither OpenSim [49, 50] nor Moco [48] provide means to personalize relevant aspects of a patient’s neuromusculoskeletal computer model or to use a personalized model to perform the treatment design optimization process.

Several existing musculoskeletal modeling software tools have sought to address these challenges to various degrees [51,52,53,54,55,56,57,58]. These tools primarily focus on personalizing the musculoskeletal model creation process, including personalizing bone geometry from medical imaging data, personalizing muscle attachment points on bones, and generating scaled generic musculoskeletal models using automated methods. The SimCP simulation platform, designed to predict post-surgery walking function for individuals with cerebral palsy, goes beyond other tools in several important ways [51]. It supports personalization of Hill-type muscle–tendon model properties, personalization of neural control properties using a muscle synergy structure, and prediction of the joint moment “capability gap” that would remain following implementation of a proposed surgical plan. Nonetheless, none of these tools can personalize multiple aspects of a patient’s neuromusculoskeletal computer model and then use the personalized model within a predictive simulation framework to design clinical interventions.

This article introduces the Neuromusculoskeletal Modeling (NMSM) Pipeline, a MATLAB-based software package that provides researchers and clinicians working together with state-of-the-art computational tools for designing effective personalized treatments for movement disorders. The software package is built on the foundation of the OpenSim musculoskeletal modeling software and adds two toolsets to OpenSim [49, 50]. The Model Personalization toolset uses patient movement data and gradient-based optimization to personalize four subcomponents of a patient’s neuromusculoskeletal computer model: 1) joint structure models, 2) muscle–tendon models, 3) neural control models, and 4) foot–ground contact models. The Treatment Optimization toolset uses a patient’s personalized model with direct collocation optimal control methods [35, 41, 48, 59,60,61,62] to predict how the patient’s neural control and/or anatomy should be changed, or how an external device (e.g., an exoskeleton) or implant should be designed and/or controlled, to maximize improvement in the patient’s post-treatment movement function. By emphasizing functionality, ease of use, flexibility, and computational speed, the NMSM Pipeline will hopefully accelerate involvement of the neuromusculoskeletal modeling research community in the computational design of clinically implementable personalized neurorehabilitation and surgical interventions.

The remainder of this article is outlined as follows. First, a design overview of the NMSM Pipeline is presented, including an overview of the philosophy and goals of the software and an introduction to the two included toolsets. Second, implementation details, design decisions, and important information for each tool are presented. Third, an example treatment design use case is provided including specifics about the use of each tool and final results. Fourth, a discussion section is included to provide information about future work, limitations, and access to the software.

Design overview

Philosophy and goals

The goal of the NMSM Pipeline software is to help neuromusculoskeletal modeling researchers engage in clinical research—and cross the valley of death [63], thereby moving neuromusculoskeletal modeling a significant step closer to becoming a clinically useful tool. To achieve this goal, the NMSM Pipeline uses a personalized medicine approach, where physics- and physiology-based personalized neuromusculoskeletal computer models constructed from patient movement data are used to design personalized treatments for movement impairments [64]. Such impairments arise from clinical conditions like osteoarthritis, stroke, cerebral palsy, spinal cord injury, Parkinson’s disease, limb amputation, and even cancer. Because the open-source OpenSim [49, 50] and commercial Anybody [65] software commonly used for musculoskeletal modeling research lack the ability to personalize neuromusculoskeletal computer models to patient movement data, these software tools cannot currently be used to optimize treatments for individual patients. To overcome this critical limitation, the open-source NMSM Pipeline software uses MATLAB (The MathWorks, Natick MA) to build on the OpenSim foundation and provides two new toolsets—a Model Personalization toolset that facilitates the creation of personalized neuromusculoskeletal computer models, and a Treatment Optimization toolset that uses the personalized models to design personalized treatments that maximize a patient’s post-treatment movement function.

To develop these two new MATLAB-based toolsets, we defined a high-level software design philosophy. First, the NMSM Pipeline was designed to be both accessible for basic users and extensible for advanced users so that neuromusculoskeletal modeling researchers of all levels can benefit from it. Second, modeling and simulation methods were based on laws of physics and principles of physiology and neuroscience rather than on machine learning, since no training data exist for post-treatment movement conditions, and machine learning methods do not, in general, extrapolate well outside the boundaries of the training data. Third, existing OpenSim functionality and model entities were utilized to the fullest extent possible to avoid duplication. Fourth, NMSM Pipeline functionality was designed to have the “look and feel” of OpenSim so that those with experience using OpenSim can quickly become proficient using the Pipeline. Fifth, the NMSM Pipeline was designed to work with any OpenSim model, providing complete flexibility in the musculoskeletal anatomy and movement situations that can be modeled and simulated.

To achieve the goal of accessibility for basic users and extensibility for advanced users, we designed NMSM Pipeline tools to work using Extensible Markup Language (XML) settings files, similar to OpenSim. This approach allows users to run each NMSM Pipeline tool with only a single line of code. For basic users, we developed two OpenSim GUI plugins — one for the Model Personalization toolset with its four tools, and one for the Treatment Optimization toolset with its three tools, where each toolset appears under the OpenSim “Tools” menu. The GUI plugin window for each tool allows users to set up tool runs by configuring the most common settings for each tool (Fig. 1), with default values being used automatically for advanced tool settings. Tool settings selected in the GUI plugin, along with default advanced tool settings, can be exported by the user to an XML settings file that can be read back into the GUI and modified as desired, just as with OpenSim tools. For advanced users, advanced tool settings are exposed in each tool’s XML settings file, again similar to OpenSim. Advanced users can create an initial XML settings file using the GUI plugin, or they can start from a pre-existing XML settings file and edit the values of basic and advanced tool settings using any text editor. For both basic and advanced users, each tool is run from the MATLAB command line using a tool-specific function that takes a tool-specific XML settings file as an input. Thus, no “Run” button is provided in any of the GUI plugin tool windows. Tool outputs can be plotted in MATLAB using included utility functions.

Fig. 1

Screenshot of the JMP GUI tool sub-window for adding a JMP task. The GUI tools allow for modifying common settings and creating XML settings files. Advanced settings can be modified directly in the XML file after creation

Software engineering principles were used for NMSM Pipeline development to facilitate extensibility. Namely, atomic functions—functions that complete only a single operation—were used extensively. All code was written using clear, English-language variable and function names and organized into logical groupings. GitHub’s source control, automated testing framework, and continuous integration tools were used to maintain this extensive project and its documentation throughout the development process. These choices make NMSM Pipeline code easy to understand and easy to modify.

Because NMSM Pipeline tools are based on laws of physics and principles of physiology and neuroscience, we selected computational methods that best facilitate the intended personalized medicine approach. Gradient-based optimization was used for all tools within the NMSM Pipeline primarily for computational speed reasons but also for the ability of gradient-based optimizers to find better solutions than the initial guess, even if the final solution is not the global optimum. For the Model Personalization toolset, built-in MATLAB optimization algorithms (primarily lsqnonlin and fmincon) were used exclusively, while for the Treatment Optimization toolset, the open-source IPOPT optimizer [66] was used within the commercial GPOPS-II direct collocation optimal control software [60]. This approach eliminates GPOPS-II and direct collocation optimal control from the Model Personalization process, making model personalization a more reliable and less complex process than it would be if direct collocation optimal control were required [41, 67, 68]. The Model Personalization toolset automatically uses MATLAB’s built-in parallel processing capabilities for calculating finite difference derivatives needed by MATLAB optimization algorithms, while the Treatment Optimization toolset automatically uses parallel processing only for repeated OpenSim analyses performed through MATLAB via C++ Mex files.

Through the utilization of OpenSim’s native models and entities, we included a number of convenient facilities in the NMSM Pipeline. OpenSim native file types were used for all data interactions including OpenSim model files (.osim), marker data files (.trc), kinematic motion and ground reaction data files (.mot), time series data files (.sto), and XML settings files (.xml). Only one new file type was developed– the NMSM Pipeline model file type (.osimx) to store personalized model properties that are not included in native OpenSim models. Similarly, easy-to-use OpenSim GUI tools such as the Scale tool, Inverse Kinematics (IK) tool, Inverse Dynamics (ID) tool, and Muscle Analysis tool were leveraged as part of the typical pipeline run-through using experimental movement data. For some NMSM Pipeline tools, groups of muscles must be identified, and the muscle groups functionality available in OpenSim models was used for this purpose to avoid the need to store this group information elsewhere.

NMSM pipeline toolsets

The NMSM Pipeline is composed of two toolsets – a Model Personalization toolset and a Treatment Optimization toolset. The Model Personalization toolset personalizes four aspects of a scaled generic OpenSim model using experimental movement data and gradient-based optimization (Fig. 2). Those four aspects are 1) joint models, 2) muscle–tendon models, 3) neural control models, and 4) ground contact models. The personalization process is important since personalized rather than generic models are needed to support model-based treatment design for clinical conditions where patients exhibit significant heterogeneity (e.g., stroke) [69, 70]. Furthermore, without a personalization process, neuromusculoskeletal models do not reliably predict internal muscle and joint contact forces, body motion, or metabolic cost [41, 71,72,73,74,75,76,77,78]. Each stage in the NMSM Pipeline Model Personalization process has been shown to play an important role in predicting these quantities reliably. [41, 71,72,73, 78, 79]

Fig. 2

Overview of the Model Personalization toolset. This toolset creates neuromusculoskeletal models with personalized joint, muscle–tendon, neural control, and ground contact parameters through the use of gradient-based optimization

The Treatment Optimization toolset uses a subject’s personalized neuromusculoskeletal model to predict their post-treatment movement function using direct collocation optimal control (Fig. 3). Over the past 10 years, musculoskeletal modeling researchers have converged on direct collocation optimal control as the best method for generating predictive simulations of human movement [35,36,37, 37, 40, 41, 48, 59, 67, 80,81,82,83,84,85,86,87,88]. The Treatment Optimization Toolset uses the MATLAB-based GPOPS-II direct collocation optimal control software [60] to perform predictive simulations of movement for novel conditions where no experimental data exist. Similar to MATLAB, GPOPS-II is inexpensive for academic use. The optimal control solver can adjust neural and/or torque control signals along with neural control, anatomical, internal implant, and/or external device parameters in a patient’s personalized neuromusculoskeletal model to achieve a specified treatment goal (e.g., achieve a desired reduction in metabolic cost). For any patient, this approach can potentially provide the best treatment prescription, the best dosage, and an estimate of the patient’s capacity for improvement.

Fig. 3

Overview of the Treatment Optimization toolset. This toolset generates a tracking simulation and two subsequent predictive simulations through the use of direct collocation optimal control

Optimization cost function terms for both toolsets were designed to minimize squared errors and/or cost term values and utilize a common formulation with physically meaningful cost function parameters. Specifically, cost function terms utilize a maximum allowable error and, in some cases, a user-provided error center (Fig. 4). For example, if a cost function term seeks to minimize errors between model and experimental joint moments, and matching to within 2 Nm is acceptable, then the user would set their cost term error center to 0 and the associated maximum allowable error to 2. Each cost function term has equal weight with all other cost function terms included in the optimal control problem formulation, where continuous cost function terms involving time series data (e.g., joint angles, joint moments) are normalized by final time so that they do not contribute to the total cost more heavily than do terminal cost function terms involving only a single quantity (e.g., metabolic cost). This method for formulating all cost function terms allows for contextualized decision-making when determining weights for different cost function terms.

Fig. 4

The generic cost function term formulation used throughout the Model Personalization and Treatment Optimization toolsets. The formulation utilizes a maximum allowable error and error center to make the calculated cost physically meaningful. If the difference between the current value and the error center is greater than the maximum allowable error, then the quantity in parentheses will be greater than 1, and squaring this quantity will amplify the cost (e.g., 1.1^2 = 1.21). In contrast, if the difference between the current value and the error center is less than the maximum allowable error, then the quantity in parentheses will be less than 1, and squaring this quantity will attenuate the cost (e.g., 0.9^2 = 0.81). For cost function terms that track time-varying experimental data, the time-varying experimental value is used implicitly as the error center. For cost function terms that minimize errors relative to a static reference value, the specified error center value is used explicitly as the error center

Neural control models for both toolsets are described using muscle synergies. Muscle synergies provide a physiologically-relevant low dimensional representation of a large number muscle activations and can be personalized to represent a specific patient’s neural control strategy [89,90,91,92,93,94,95,96,97]. Each muscle synergy is composed of a time-varying synergy activation and an associated time-invariant synergy vector. Each synergy vector contains a weight for each muscle in the model that defines how the associated synergy activation contributes to the total activation of the muscle. Muscle synergies are commonly calculated using non-negative matrix factorization [89, 90], and typically between three and six muscle synergies are required to account for over 95% of the signal variability for common tasks such as walking, running, and reaching. [98]

Implementation details

Model personalization

The Model Personalization toolset calibrates parameter values in an OpenSim neuromusculoskeletal model so that the model closely reproduces a wide range of experimental movement data. The toolset consists of four tools and was designed so that the tools can be used largely independently from one another, though a preferred tool workflow exists. Personalization is likely to be necessary if the musculoskeletal model is to be used for clinical treatment design purposes [34]. It also improves models for use in predictive simulations by maximizing the self-consistency of collected experimental data and easing the task of finding a realistic, dynamically consistent motion.

The four tools of the Model Personalization toolset personalize different anatomical, physiological, and neurological aspects of a patient’s OpenSim model. Joint Model Personalization (JMP) calibrates the properties of OpenSim constraint-based joint models, Muscle–tendon Model Personalization (MTP) calibrates the properties of Hill-type muscle–tendon models, Neural Control Model Personalization (NCP) calibrates the properties of muscle synergy control models, and Ground Contact Model Personalization (GCP) calibrates the properties of elastic foundation foot–ground contact models. OpenSim models personalized using the Model Personalization toolset can be used in both the NMSM Pipeline’s Treatment Optimization Toolset and, in the near future, OpenSim Moco. For use with OpenSim Moco, an NMSM Pipeline.osimx model file must be converted to an OpenSim.osim model file using an provided MATLAB model conversion utility program.

The Model Personalization toolset requires input data generated using standard OpenSim tools. First, a generic musculoskeletal model is scaled using the OpenSim Scale Model tool. To facilitate modeling and simulation of full-body walking motions, we recommend starting with the RCNL2024 OpenSim model (distributed with the NMSM Pipeline software), which is a modified version of the Rajagopal full-body OpenSim model [99] where the ankle joint kinematic model has been changed to match van den Bogert et al. (1994) [100] and the knee joint kinematic model changed to eliminate all rotational offsets while producing the same motion. A scaled generic model is one of the required inputs for the JMP tool. After the joints of the model are personalized to the individual’s functional joint axes using the JMP tool, the user must calculate IK motions and ID loads using the OpenSim Inverse Kinematics and Inverse Dynamics tools, respectively. The user must then run the Muscle Analysis tool to calculate muscle–tendon lengths and muscle moment arms if the NMSM Pipeline MTP and NCP tools are to be used. For the JMP tool, marker files can be manually concatenated to personalize joint parameters, values that characterize the transformation of the parent and child bodies of a joint, using multiple motion trials (e.g., marker data from isolated ankle, knee, and hip motion trials can be concatenated with marker data from a walking motion trial to calibrate the functional axes and joint centers of all lower body joints simultaneously). The MTP and NCP tools accept data from one or more motion trials, and the GCP tool accepts data from only a single motion trial, ideally the same trial to be used for Treatment Optimization. If multiple motion trials are used with either the MTP or NCP tools, IK motions, ID loads, muscle–tendon lengths and velocities, and muscle moment arms need to be calculated for each motion trial.

Joint model personalization

The Joint Model Personalization tool optimizes joint parameters, body scaling, and marker placement to minimize IK marker distance errors. [100,101,102,103,104,105,106]. Reducing inverse kinematics marker distance errors reduces downstream errors in calculated inverse dynamic joint moments [72], muscle–tendon lengths and velocities, muscle moment arms, and ultimately muscle activations and forces. These quantities are used by subsequent Model Personalization tools.

JMP tool inputs are a scaled generic OpenSim model and marker data from one or more motion trials, where the tool settings file allows the user to specify one or more tasks (Fig. 5). The design variables available in a JMP settings file include joint locations and orientations in parent and child bodies, body scale factors applied uniformly, and marker location offsets for each direction independently. The only cost function term is the sum of the squared marker distance errors normalized by the motion trial duration and number of markers. Only one constraint is included in this optimization problem to prevent extreme normalized muscle fiber length values. The output of JMP is a new OpenSim model with optimized joint parameters, body scale factors, and/or marker placements.

Fig. 5

Flowchart showing the inputs to and outputs from the JMP tool. The diagram key also applies to Figs. 6 through 9

Several important considerations should be kept in mind when using the JMP tool. A JMP run can be split into multiple tasks, each with its own set of design variables and an associated marker file. Careful JMP setup is required to avoid obtaining non-anatomical results. Such a result could be found by personalizing a joint parameter for a coordinate that does not move through sufficient range of motion (> ~ 25 degrees) during the motion trial [107]. This issue can be rectified by personalizing the associated joint parameters with separate tasks using isolated joint motion trials, or by concatenating data from an isolated joint trial with data from the desired functional trial (e.g., walking). Another cause of non-anatomic results is selection of joint parameters that produce relative bone motions that are inconsistent with joint function. For example, the parent frame of the hip joint is often defined to be at the center of the acetabulum, while the child frame is often defined to be at the center of the femoral head. Allowing joint parameters to modify the location of the hip joint in either the parent or the child frame could move the center of the femoral head out of the center of the acetabulum, breaking that anatomical relationship. In addition, muscle paths and/or relative body sizes may also become non-anatomical due to inclusion of inappropriate design variables. JMP tool results should be visualized in the OpenSim GUI and problem formulations modified to achieve an anatomically realistic post-JMP model.

Muscle–tendon model personalization

The Muscle–tendon Model Personalization tool finds an optimal set of subject-specific muscle–tendon properties and muscle activations from EMG, joint kinematic, and joint moment data by balancing optimization cost function terms related to muscle properties, similarity of properties among grouped muscles, and matching of EMG-driven and experimental inverse dynamics joint moments [75, 76, 108, 109]. Muscle activation and force predictions are sensitive to optimal muscle fiber length and tendon slack length [71, 110,111,112]. Therefore, reliable personalization of these parameters is essential for generating reliable predictions of muscle activations and forces during predictive simulations of movement.

The inputs to the MTP tool are a post-JMP OpenSim model as well as IK motion, ID load, muscle–tendon length and velocity, and muscle moment arm data from one or more motion trials of interest (Fig. 6). The design variables that can be optimized include activation dynamics and muscle–tendon length parameters for each muscle [41, 75, 113]. Cost function terms used by the MTP tool minimize muscle–tendon model parameter deviations from initial values, ID joint moment tracking errors, passive muscle forces, inconsistencies among grouped normalized muscle fiber lengths, inconsistencies among grouped muscle excitation scale factors, and inconsistencies among grouped electromechanical delays. No constraints are applied in this optimization problem. The outputs of the MTP tool are a .osimx model file containing optimized muscle–tendon properties and estimated muscle activations, both for use by other tools in the NMSM Pipeline.

Fig. 6

Flowchart showing the inputs to and outputs from the MTP tool. The diagram key is included in Fig. 5

Several important considerations should be kept in mind when using the MTP tool. Muscle-tendon Length Initialization is included as an option within the MTP tool and is typically used to find initial values for muscle–tendon length parameters that act as reasonable reference values for cost function deviation terms involving these parameters. Muscle-tendon Length Initialization also estimates the maximum isometric force of each muscle based on published regression relationships [114]. Since scaling maximum isometric force and muscle excitation at the same time produces a non-unique muscle force solution, and since muscle excitations are scaled using EMG scale factor design variables, maximum isometric force is held constant during an MTP tool run [71]. Synergy Extrapolation (SynX) [115, 116] is also included in the MTP tool to estimate the activations of muscles without experimental EMG data. SynX uses muscle synergy concepts to estimate missing activations and is more accurate than traditionally used static optimization at predicting missing EMG signals [117]. A typical use case for SynX is estimating the activations of deep muscles from which EMG data cannot be easily collected. SynX can use either principal component analysis or non-negative matrix factorization to calculate the muscle synergies needed to predict muscle activations associated with missing EMG data. The MTP tool uses a Hill-type muscle–tendon model with a rigid tendon and assumes that muscles function primarily on the ascending region of their normalized active force–length curves [118]. Users can modify cost terms to specify the region of the normalized force–length curve over which muscles are assumed to operate. Some MTP cost function terms, such as inconsistencies among grouped normalized muscle fiber lengths, calculate the difference between the grouped values relative to their mean. These cost terms are intended to promote similarity in muscle properties within a group under the assumption that similar muscles, such as the three uniarticular vastus muscles, should have similar properties. These muscle groups must be defined within the OpenSim model. Similar cost function terms are also present in NCP.

Neural control model personalization

The Neural Control Model Personalization tool finds muscle synergies that are as consistent as possible with ID joint moments and, when available, MTP-estimated muscle activations. The NCP tool fits muscle synergies at the level of muscle activations (i.e., after electromechanical delay and activation dynamics) for regions of the body where either all muscle activations are available from the MTP tool (e.g., the right lower extremity) or no muscle activations are available [119, 120]. The tool can also be used to find muscle synergies in multiple regions of the body simultaneously — some with and some without previously calculated muscle activations. During the NCP optimization process, grouped muscle properties are maintained and activation-related cost function terms are minimized. Muscle synergies are typically calculated from muscle activations alone using non-negative matrix factorization, and thus the muscle activations reconstructed from the resulting synergies will not closely reproduce ID joint moments when input into a subject’s musculoskeletal model. NCP rectifies this issue by optimizing synergy activations and synergy vectors to be consistent with both ID joint moments and available muscle activations simultaneously.

The inputs to the NCP tool are a post-JMP OpenSim model as well as data for IK motions, ID loads, muscle–tendon lengths and velocities, muscle moment arms, and, optionally, MTP-calculated muscle activations from one or more motion trials of interest (Fig. 7). The OpenSim model must contain muscle groups that define lists of muscles for which synergy sets are to be fitted (e.g., all right leg muscles). The design variables for the NCP tool are a set of time-varying synergy activations and an associated set of time-invariant synergy vectors. The cost function terms for the NCP tool minimize ID joint moment matching errors, muscle activation matching errors for muscles used in MTP, muscle activations for muscles not used in MTP, inconsistencies among grouped normalized muscle fiber lengths, and inconsistencies among grouped muscle activations. The NCP tool includes a constraint that the sum of the muscle weights in each synergy vector must equal a constant times the number of muscles described by that synergy. Upon output, the synergy vectors and synergy activations are normalized such that the largest weight in each synergy vector is one. The NCP tool produces a set of synergy groups inside a new or pre-existing NMSM Pipeline model file (.osimx) as well as data files (.sto) containing synergy activations and their associated synergy vectors.

Fig. 7

Flowchart showing the inputs to and outputs from the NCP tool. The diagram key is included in Fig. 5

Several important considerations should be kept in mind when using the NCP tool. Synergy-driven Tracking Optimization (TO), part of the Treatment Optimization toolset, requires that NCP tool results be used as the initial guess. NCP results that distribute errors equally between MTP muscle activation matching and ID joint moment matching will facilitate the process of finding good synergy-driven TO results. In addition, to facilitate obtaining good NCP results, users should perform their own synergy analysis on the muscle activations produced by MTP to determine the optimal number of synergies to use. The number of synergies should be selected to achieve the desired variability accounted for (VAF), which is often set at 90 or 95% [41, 93, 121, 122]. The NCP tool can be run without having MTP tool results as inputs through the use of the activation minimization cost function term. Since running the NCP tool without MTP tool results produces a non-unique solution, the maximum allowable error of the activation minimization cost function term can be set to achieve a unique solution that contains activations of expected magnitude. Shared synergy vectors between synergy groups can also be enforced by the NCP tool. This feature allows for two muscle synergy groups, such as those for right and left leg muscles, to have nearly identical synergy vectors.

Ground contact model personalization

The Ground Contact Model Personalization tool finds physical properties for elastic foundation foot–ground contact models that closely reproduce experimental ground reaction forces and moments while allowing slight changes in foot kinematics. The elastic foundation is composed of a uniform grid of linear springs with nonlinear damping and friction placed across the bottom of the foot. The tool can personalize ground contact models for each foot separately, where physical properties will differ between the two feet, or both feet together, where physical properties will be the same for the two feet. Currently, the tool works with only a two-segment foot model composed of a hindfoot segment (including calcaneus, cuboid, navicular, cuneiform, and metatarsal bones) and a toes segment (including phalange bones). While generic foot–ground contact models can reproduce experimental ground contact forces accurately, they reproduce experimental ground reaction moments less accurately, introducing errors into the lower body joint moments generated by predictive simulations of walking [79]. Personalized foot–ground contact models can overcome this issue by reproducing experimental ground reaction moments as well.

The inputs to the GCP tool are a post-JMP lower-body or full-body OpenSim model along with IK motion and associated ground reaction data for a walking trial, where the IK and ground reaction data must possess the same time increments (Fig. 8). The ground reaction data must contain three forces, three moments, and a point about which the moments are calculated, which is typically either a time-invariant point fixed in the laboratory or a time-varying center of pressure. The tool automatically strips off both feet from the lower-body or full-body OpenSim model and creates a 6 degrees of freedom (DOF) custom joint connecting each specified hindfoot body to ground. The design variables for the GCP tool are stiffness coefficients unique to each spring with a nonlinear viscous damping coefficient, a dynamic friction coefficient, a viscous friction coefficient, and a spring resting length common to all springs. Additional design variables are B-spline nodal points defining small deviations from experimental foot kinematics with respect to ground. The cost function terms for the GCP tool include minimization of foot marker position errors, foot marker velocity errors, hindfoot and toes segment rotational coordinate errors, hindfoot segment translational coordinate errors, vertical ground reaction force errors, horizontal ground reaction force errors, ground reaction moment errors, and spring coefficient deviations from neighboring spring coefficient values. No constraints are present in this optimization problem. The output of the GCP tool is a new NMSM Pipeline model file (.osimx) containing foot–ground contact model parameter values.

Fig. 8

Flowchart showing the inputs to and outputs from the GCP tool. The diagram key is included in Fig. 5

Several important considerations should be kept in mind when using the GCP tool. The tool can reproduce all six components of ground reactions with good fidelity. To facilitate reproducing ground reaction moments, we use the time-varying projection of a midfoot superior marker onto the ground as the point about which ground reaction moments are calculated. This automatic recalculation of the ground reaction moments helps keep the moments well-scaled during a GCP run, since the moment arm of the ground reaction force vector never becomes excessively large. The center of pressure is not used for joint moment matching since it becomes inaccurate at the transitions into and out of contact due to division by small vertical ground reaction force values. At this time, the GCP tool requires the use of a two-segment foot model with the distal toes segment connected to the hindfoot segment via an oblique pin joint. Before personalization, the GCP tool automatically extracts two-segment OpenSim foot models from the post-JMP lower-body or full-body OpenSim model and automatically calculates hindfoot kinematics with respect to ground that are consistent with the input lower-body or full-body kinematics. Since the foot–ground contact model is sensitive to foot kinematic changes below the measurement accuracy of a markered motion capture system, the GCP tool allows for small changes to the foot motion trajectory to facilitate ground reaction matching. Each foot model possesses seven coordinates: six ground-to-hindfoot coordinates via a 6 DOF custom joint and one toes coordinate via a pin joint. Coordinate deviations are parameterized with B-spline nodes to reduce the number of design variables. The GCP tool automatically estimates the number of B-spline nodes needed based on a specified lowpass filter cutoff frequency. Another important feature of the GCP tool is its ability to use a different time window for calibrating each foot. This feature facilitates the use of the GCP tool with ground reaction data obtained from labs possessing only two in-ground force plates. The GCP tool also allows for symmetric physical properties between contact surfaces on the two feet, and users can modify the number of springs along each dimension of the grid.

Treatment optimization

The Treatment Optimization toolset generates predictive simulations of a patient’s post-treatment movement function by optimizing specified treatment design parameters. The toolset consists of three tools designed using a “theme and variation” approach, where the tools are intended to be used in a specific order, with each tool serving a distinct purpose. Each tool uses the GPOPS-II direct collocation optimal control software for MATLAB [60] and maintains a consistent structure for data inputs, problem design, cost function terms (see Table S1 for an exhaustive list), constraint terms (see Table S2 for an exhaustive list), and outputs, with variations (Fig. 9). The Tracking Optimization tool (TO) generates “starting point” predictive simulations of patient movement that are dynamically consistent and closely reproduce all relevant experimental movement data available from the patient (e.g., joint motion, joint moment, ground reaction, and EMG data). The Verification Optimization tool (VO) starts from the Tracking Optimization solution and verifies that the controls found produce the correct motion without tracking the motion directly, which ensures that the subsequent treatment design problem is well formulated before simulated treatments are explored. Finally, the Design Optimization tool (DO) starts from the Verification Optimization solution, applies and (if desired) optimizes simulated treatments, and predicts the patient’s resulting post-treatment movement along with associated controls and loads. All three tools work with torque controlled, muscle synergy controlled, and combined torque and muscle synergy controlled problem formulations. For problems utilizing muscle synergy control, all three tools require surrogate musculoskeletal geometry created from the final OpenSim model produced by the Model Personalization process. The surrogate geometry creation process is semi-automated and involves using a provided MATLAB utility function to generate sample model poses spread around a specified nominal motion, performing a single Muscle Analysis on the sampled model poses using the OpenSim GUI, and fitting surrogate muscle-tendon length, velocity, and moment arm models from the Muscle Analysis outputs automatically at the start of a TO, VO, or DO tool run.

Fig. 9

Flowchart showing the inputs to and outputs from the tools included in the Treatment Optimization toolset. Since these tools follow a theme and variation approach, various data combinations can be included or excluded based on the goals of the user. The diagram key is included in Fig. 5

Several important considerations should be kept in mind when working with the Treatment Optimization toolset. Direct collocation optimal control problems are solved by modifying states and controls simultaneously to minimize an objective function and fulfill constraints. For Treatment Optimization, the states consist of joint positions and velocities for selected OpenSim model coordinates and the controls consist of corresponding joint accelerations (due to the use of an implicit dynamics formulation [59]), torque controls, synergy controls, or both. All Treatment Optimization tools accept both a “tracked quantities” and an “initial guess” directory. The tracked quantities directory contains reference data (e.g., experimental joint positions, joint moments, muscle activations, and ground reactions) to be tracked by specified cost function terms while the initial guess directory contains data specifying the initial guess for the problem’s states and controls. Treatment Optimization XML settings files allow for joint position, velocity, and acceleration range scale factors as well as muscle synergy and torque controller range scale factors. These range scale factors describe the search space of the optimal control problem by confining the states and controls to be within a certain range relative to their initial range (see Fig. 10). After a completed run, all Treatment Optimization tools output the same files to a user-specified results directory. Specifically, the solution’s states and controls, along with IK motions, ID loads, and ground reactions calculated from the solution’s states and controls, are output to the specified directory. Saved states and controls can then be passed to the next Treatment Optimization tool as an initial guess to pick up where the previous tool left off.

Fig. 10

Illustration showing how range scale factors are used to set bounds on the states and controls within GPOPS-II. The range scale factor value is used to set the maximum and minimum searchable spaces for the states and controls as offsets beyond the initial range. In this figure, a range scale factor of 0.5 was used

Because Treatment Optimization tools use implicit dynamics [41, 59], several constraint terms must be included for successful use of the toolset. One constraint is that the joint acceleration controls must be the first time derivatives of the joint velocity states, and the joint velocity states must be the first time derivatives of the joint position states. These constraints are handled internally by GPOPS-II. Another constraint is that the net joint loads produced by the torque and synergy controls must be consistent with the model’s inverse dynamic loads. This constraint must be handled by the user by including kinetic consistency constraint terms for all controlled coordinates included in the problem formulation. For models with non-physical residual forces and/or torques acting on the root segment (i.e., the segment connecting the model to ground), if the user wants their Treatment Optimization solutions to be dynamically consistent, then root segment residual load constraints must be added by the user to ensure that the final residual loads are close to zero. In addition, if the motion of interest is expected to be periodic, such as for walking, running, cycling, or stair-climbing, then joint position and velocity along with ground reaction force and moment periodicity constraints can be used to ensure that a near-periodic motion is found.

Tracking optimization

The Tracking Optimization tool uses a personalized model to produce a dynamically consistent movement simulation that closely reproduces all available experimental motion data, including joint motions, joint moments, ground reaction forces and moments, and muscle activations. To achieve a dynamically consistent motion, the tool spreads out matching errors between the different experimental quantities based on user-specified maximum allowable errors. The tool accepts a post-JMP OpenSim model (.osim file) and personalized NMSM Pipeline model (.osimx file) along with experimental IK motions, ID loads, ground reactions, muscle–tendon lengths and velocities, muscle moment arms, and, if using synergy controls, NCP results for the trial of interest. The following cost function terms are commonly included to achieve a successful Tracking Optimization run: generalized coordinate tracking, generalized speed tracking, inverse dynamic moment tracking, muscle activation tracking, external force tracking, external moment tracking, body orientation relative to ground tracking, and center of pressure tracking. Common constraint terms include kinetic consistency, state position periodicity, state velocity periodicity, ground reaction force periodicity, ground reaction moment periodicity, and root segment residual load bounding.

Several important considerations should be kept in mind when using the Tracking Optimization tool. Finding a good result with the Tracking Optimization tool is the biggest challenge in the NMSM Pipeline. The goal is to trade off errors between tracked experimental data sources by modifying the maximum allowable error in each cost function term based on the importance of and confidence in the respective experimental data. Higher quality Model Personalization results and highly consistent tracked experimental quantities (e.g., tracked muscle activations should closely reproduce tracked ID moments) make finding a good Tracking Optimization solution easier.

Verification optimization

The Verification Optimization tool is used to perform a “dry run” Design Optimization without including any “design elements”, the goal being to ensure a good initial guess and to verify the appropriateness of the optimal control problem formulation. The tool accepts the same inputs as the TO tool, but in general, the results directory of a TO run is used for both the initial guess and tracked quantities. The cost function terms for the VO tool include controller tracking and, for uncontrolled joints, generalized coordinate tracking. Similar to the TO tool, the VO tool includes constraint terms for kinetic consistency, state position periodicity, state velocity periodicity, and root segment residual load bounding. If a VO problem is well formulated, it should converge quickly since the initial guess is already at the optimal solution.

Several important considerations should be kept in mind when using the Verification Optimization tool. Users must decide whether each coordinate should be controlled or tracked. If controlled, the coordinate should be assigned the same controller (synergy, torque, or both) as in the TO tool, and the controls should be tracked by the VO tool. Conversely, if the coordinate does not have a controller or has a non-essential controller, then the coordinate should be tracked by the VO tool. An example of a non-essential controller is a torque controller on an upper body joint for a full body gait model with only lower body muscles. The results of a VO run should closely match the outputs of the TO run used to define the initial guess and tracked quantities. If the results do not match closely, then a problem formulation issue is likely present, and users should review their TO and VO settings. Examples of potential problem formulation issues are having a generalized coordinate that is neither controlled nor tracked, having a generalized coordinate that is both controlled and tracked, or having a generalized coordinate that is tracked with a different maximum allowable error than was used by TO. Each of these situations can produce a VO solution that differs from the TO initial guess. After an adequate VO solution is found, users can move on to the DO tool by adding a “treatment” to their VO problem formulation.

Design optimization

The Design Optimization tool predicts or optimizes how a planned treatment will affect a patient’s post-treatment movement function. The tool accepts the same inputs as the other Treatment Optimization tools, but typically VO results are used as the initial guess. Similar to VO, DO allows for controller tracking and generalized coordinate tracking alongside other DO-specific cost function terms. Similar to TO and VO, DO can use kinetic consistency, state position periodicity, state velocity periodicity, and root segment residual load constraint terms as well as other constraints to facilitate the design of a personalized intervention.

Several important considerations should be kept in mind when using the Design Optimization tool. Unlike the other Treatment Optimization tools, the DO tool allows for both fixed and free final time problem formulations. For free final time problems, the time vector can be constrained to a user defined range and tracked control quantities are stretched or compressed in time to match the current time vector. In addition to all cost function and constraint terms available in the TO and VO tools, the DO tool includes a number of additional built-in cost function terms. These terms include minimizing the following types of errors with respect to a target value: braking and propulsive impulse, metabolic cost per unit of time or distance [123], whole-body angular momentum, and joint mechanical energy generation or absorption. Three key features of the DO tool are support for model modification functions, user-defined cost function terms, and user-defined constraint terms. By employing model modification functions alongside GPOPS-II’s static parameters feature, users can change the OpenSim musculoskeletal model, modify the neural control model, or adjust parameter values in an assistive device. Possible model modifications include changing bone geometry (e.g., osteotomy parameters), changing muscle–tendon properties (e.g., strengthening muscles or adjusting optimal muscle fiber lengths), changing neural control properties (e.g., altering synergy vector weights), modifying ground contact model parameters, optimizing muscle attachment locations or muscle wrapping surface parameters, redefining joint parameters, and modifying external device design and/or control parameters. User-defined cost function and constraint terms allow users to write custom MATLAB functions, with specific function signatures, that are included as additional cost function and constraint terms given the current iteration’s states, controls, and static parameters.

Example application

Treatment design rationale

To demonstrate the capabilities of the NMSM Pipeline, we used it to investigate a hypothetical treatment to improve walking speed and function for an individual post-stroke. Published studies have hypothesized that time-varying synergy activations originate in the brain while time-invariant synergy vectors reside in the spinal cord [124, 125]. If this hypothesis is true, then one might expect a stroke to affect a patient’s synergy activations but not synergy vectors [125]. This expectation in turn suggests that an individual post-stroke may have similar synergy vectors between the two legs but lack the ability to recruit one or more of them effectively on the paretic side. Furthermore, published studies have reported that a stroke results in altered neural coordination even on the non-paretic side [126,127,128]. Taken together, these observations lead to the hypothesis that a stroke may impair recruitment of one or more synergy activations in the paretic leg, resulting in altered recruitment of one or more synergy activations in the non-paretic leg, and that restoration of these impaired and altered synergy activations could result in improved self-selected walking speed and metabolic cost. [128,129,130]

To investigate this hypothesis, we used the NMSM Pipeline to design a hypothetical treatment that would modify paretic and non-paretic leg synergy activations for an individual post-stroke with the goal of increasing his self-selected walking speed without increasing his metabolic cost. Previously published experimental walking data collected from an individual post-stroke with a self-selected walking speed of 0.5 m/s were used for the investigation [41]. Research has shown that both healthy individuals and individuals post-stroke tend to choose a self-selected walking speed that achieves a desired metabolic cost per unit time [131]. Furthermore, to achieve successful community ambulation (as defined by the walking speed required to cross a typical street before the light changes [132]), individuals should be able to walk at a self-selected speed of 0.8 m/s [132, 133]. Thus, the specific goal of the hypothetical treatment was to increase the subject’s self-selected walking speed to 0.8 m/s while maintaining the same metabolic cost per unit time estimated for the subject’s original self-selected walking speed of 0.5 m/s.

Experimental data collection

As described in Meyer et al. (2016) [41], the experimental data used for the present study were collected from a high-functioning hemiparetic subject with chronic stroke-related walking dysfunction (sex male, age 79 years, LE Fugl-Meyer Motor Assessment 32/34 pts, right-sided hemiparesis, height 1.7 m, mass 80.5 kg). Experimental data collection was approved by the University of Florida Health Science Center Institutional Review Board (IRB-01) and the Malcom Randall VA Medical Center Research and Development Committee, and the subject gave written informed consent prior to participation. Motion capture (Vicon Corp., Oxford, UK), ground reaction (Bertec Corp., Columbus, OH, USA), and electromyography (EMG) data (Motion Lab Systems, Baton Rouge, LA, USA) were collected simultaneously while the subject walked on a split-belt instrumented treadmill with belts tied. The subject walked at his self-selected speed of 0.5 m/s and fastest comfortable speed of 0.8 m/s. To facilitate ground contact model personalization, the subject wore Adidas Samba Classic sneakers, which possess a flat sole with neutral midsole. Motion capture data utilized a full-body marker set with one marker at each toe tip, three markers on each hindfoot, shank, and thigh, three markers on the pelvis, and one marker on each acromion, elbow, and wrist (Reinbolt et al., 2005; Fregly et al., 2007). Ground reaction data from each treadmill force plate included three force and moment components along with the location of the fixed point used for moment calculations. EMG data included a combination of surface and fine-wire electrodes collected from 16 muscles on each leg, including deep muscles adductor longus, iliopsoas, tibialis posterior, extensor digitorum longus, and flexor digitorum longus.

Processing of marker motion, ground reaction, and EMG data followed a previously published protocol [41]. Marker motion and ground reaction data were filtered at a variable cut-off frequency of 7/tf Hz, where tf is the period of the gait cycle being analyzed, using a fourth-order zero phase lag Butterworth filter [134]. EMG data were high-pass filtered at 40 Hz, [75] demeaned, rectified, and low-pass filtered at a variable cut-off frequency 3.5/tf Hz using a fourth-order zero phase lag Butterworth filter. Each processed EMG signal was normalized by the maximum value over all walking cycles and then offset so that the minimum value was zero. For both walking speeds, the single most periodic walking cycle (heel strike to heel strike for the paretic right leg) was selected for analysis, and all data for this one cycle were resampled to 101 time points. After resampling, EMG data were padded with 200 ms (18 time points) before the start of the cycle to account for electromechanical delay.

Model personalization process

We utilized all four tools in the Model Personalization toolset to personalize a full-body neuromusculoskeletal model of the subject being studied. The starting point for the model personalization process was a scaled RCNL2024 full-body OpenSim model, where scaling was performed using the OpenSim Model Scaling tool with experimental marker data collected from the subject during a static standing trial. All subsequent NMSM Pipeline tool runs were performed on a Dell PC workstation possessing 48 cores.

First, the JMP tool was used to personalize lower body joint parameters and lower body scaling by minimizing IK marker distance errors. The selected 0.8 m/s gait cycle marker data were used for all JMP tool tasks. To avoid non-anatomical muscle paths around the hip complex, the left and right hip joint parameters were not personalized. The left and right knee and ankle joints were personalized, but the subtalar joint was not personalized because the joint did not move through a sufficient range of motion to personalize the joint’s parameter values accurately. The JMP tool reduced marker distance errors for the thigh, shank, and foot markers by 48.4% for the gait cycle of interest. Additional JMP results are included in Table 1 and in the supplementary material. Convergence was achieved after approximately 90 min of CPU time.

Table 1 Marker distance errors before and after Joint Model Personalization

Second, the MTP tool was used to personalize the properties of all lower body muscles by minimizing errors in matching inverse dynamic joint moments and changes in initial muscle–tendon parameter values while also predicting missing EMG signals for deep or “small” muscles (e.g., gemellus). Before the MTP tool was run, the OpenSim Inverse Kinematics, Inverse Dynamics, and Muscle Analysis tools were run from the OpenSim GUI using the selected 0.8 m/s walking cycle data. Using the outputs of these GUI tools, an NMSM Pipeline utility script processed EMG data, calculated muscle–tendon velocity from muscle–tendon length data, and automatically organized the trial data into the directory structure required by the MTP tool. Muscle groups were added to the post-JMP OpenSim model to group muscles based on normalized fiber length similarity and activation similarity. Muscle–tendon Length Initialization was used to find initial values for optimal muscle fiber length and tendon slack length. At that point, muscle–tendon model and activation dynamics parameter values were personalized. This process was repeated for each leg separately via separate runs of the MTP tool. After the MTP tool was run for both legs, the left leg was found to match lower body joint moments with an average root mean square (RMS) error of 3.5 Nm and the right leg with an average RMS error of 2.1 Nm. Joint-specific moment matching results are included in Table 2 with additional plots available in the supplementary material. A custom script was used to combine the left and right leg results for use by the NCP tool. Convergence was achieved after approximately 30 min of CPU time.

Table 2 Lower body joint moment matching errors following Muscle–tendon Model Personalization

Third, the NCP tool was used to calculate lower body muscle synergies (i.e., synergy vectors and activations) by minimizing errors in matching inverse dynamic joint moments and muscle activations produced by the MTP tool. Given our treatment design rationale above, we assumed that the subject’s synergy vectors were the same while his synergy activations were different for the two legs. As a result, the NCP tool was used with the bilateral synergy vector symmetry option enabled. Synergy sets were defined as muscle groups denoting the right and left legs. A preliminary synergy analysis performed on the subject’s processed EMG data indicated that six synergies were sufficient to account for signal variability even with the added synergy vector symmetry constraint. Muscle-tendon Length Initialization was not performed as part of this NCP tool run because all muscles had MTP-derived muscle activations and thus the activation minimization term was not used. The muscle synergies found by the NCP tool reproduced all lower body joint moments with an average RMSE of 2.2 Nm and all MTP-derived muscle activations with an average RMSE of 0.034. Additional NCP results are included in Table 3 and in the supplementary material. Convergence was achieved after approximately 3 h of CPU time.

Table 3 Lower body joint moment matching errors following Neural Control Model Personalization

Fourth, the GCP tool was used to personalize foot–ground contact model parameter values by minimizing errors in matching ground reaction forces, ground reaction moments, and foot kinematics. A contact element grid size of 5 by 11 was selected to match previous successful investigations [79]. The cost functions for the three tasks performed by the tool were dominated by terms that emphasized matching the vertical GRF, vertical and horizontal GRF, and vertical and horizontal GRF as well as ground reaction moments, respectively, with each task starting from the solution of the previous task to track additional quantities dependent on the previous quantities. To facilitate smooth ground reaction forces during Treatment Optimization, we used the viscous friction coefficient as a design variable and turned off dynamic friction. The viscous friction model penalizes excessive foot velocity with respect to the ground during ground contact, which could lead to excessive horizontal GRFs in predictive simulations. The GCP tool produced a personalized non-uniform distribution of spring stiffness coefficients and matched ground reaction forces, ground reaction moments, foot positions, and foot orientations with RMSEs of 4.6 N, 2.1 Nm, 6.4 mm, and 1.7 degrees, respectively. [135, 136] Additional GCP results are included in Table 4 and in the supplementary material. Convergence was achieved after approximately 4 h of CPU time.

Table 4 Ground reaction force/moment and foot position/orientation errors following Ground Contact Model Personalization

Treatment optimization process

After personalizing the subject’s neuromusculoskeletal model using all four Model Personalization tools, we ran all three Treatment Optimization tools to predict how underutilized paretic and non-paretic leg synergy activations should be upregulated to achieve a 0.8 m/s self-selected walking speed without increasing metabolic cost per unit time above the subject’s 0.5 m/s level.

First, the TO tool was used to generate a dynamically consistent full body walking motion at 0.8 m/s that closely reproduced the subject’s joint motion, joint moment, ground reaction, and muscle activation data simultaneously using muscle synergy controls for the lower body and joint torque controls for the back and upper body. Since synergy-driven TO requires an initial guess for synergy vectors and activations, the XML settings file for the TO tool was configured such that the initial guess directly was defined to be the NCP results directory while the tracked quantities directory was defined to be the preprocessed data directory created after running the JMP tool. Other information, including an NMSM Pipeline model file (.osimx), the post-JMP OpenSim model, and the specified output directory, were also included in the TO tool XML settings file. All lower body joints were controlled by synergy controllers except for the toes coordinate, which was left uncontrolled, and all other joints were controlled by torque controllers. Cost function terms tracked joint positions and velocities found by a post-JMP IK analysis, joint moments found by a post-JMP ID analysis (excluding pelvis residual loads), muscle activations found by MTP, and external loads consisting of processed experimental ground reaction data. Constraint terms included bounded kinetic consistency errors (0.01 Nm) for all controlled joints, root segment residual loads for the pelvis (1 N for forces and 0.1 Nm for moments), and joint position periodicity errors (0.05 rad for joints below the pelvis and 0.1 rad for joints above it). The solution produced by the TO tool matched tracked joint positions with an average RMSE of 2.1 degrees, joint moments with an average RMSE of 5.1 Nm, ground reaction forces with an average RMSE of 25.3 N, ground reaction moments with an average RMSE of 4.0 Nm, and muscle activations with an average RMSE of 0.038. Additional TO results are included in Table 5 and in the supplementary material. Convergence was achieved after 305 IPOPT iterations, which required approximately 770 min of CPU time.

Table 5 Lower body Tracking Optimization errors relative to experimental data

Second, the VO tool was used to verify that tracking the muscle synergy controls and upper body plus toes joint kinematics found by the TO tool reproduced the joint motions, joint moments, and ground reactions predicted by the TO tool without tracking those quantities explicitly. For the VO run, the results of the previous TO run were used for both the initial guess and tracked quantities data. Cost function terms tracked TO-found synergy controls and TO-found torque-controlled and uncontrolled joint positions. Constraint terms included bounded kinetic consistency errors (0.01 Nm) for all controlled joints, root segment residual loads for the pelvis (1 N for forces and 0.1 Nm for moments), and joint position periodicity errors (0.05 rad for joints below the pelvis and 0.1 rad for joints above it). The solution produced by the VO tool matched tracked joint positions with an average RMSE of 0.1 deg, joint moments with an average RMSE of 0.1 Nm, ground reaction forces with an average RMSE of 0.7 N, ground reaction moments with an average RMSE of 0.1 Nm, and muscle activations with an average RMSE of 0.01. Additional O results are included in Table 6 and in the supplementary material. Convergence was achieved after 101 IPOPT iterations, which required approximately 170 min of CPU time. Since the VO results closely matched the TO results, we were confident that the controller and joint position tracking cost function terms and associated constraint terms selected for the Verification Optimization could be used as the basis for the Design Optimization.

Table 6 Lower body Verification Optimization errors relative to Tracking Optimization results

Finally, the DO tool was used to explore whether upregulating underutilized synergy activations in the paretic leg and downregulating overutilized synergy activations in the non-paretic legs could produce a faster self-selected walking speed of 0.8 m/s while maintaining the same metabolic cost per unit time calculated for the subject’s original self-selected walking speed of 0.5 m/s. The hypothesis in this problem formulation is that large differences in synergy activation amplitudes between the two legs are indicative of impaired paretic leg synergy activations that are smaller than desired and compensatory non-paretic leg synergy activations that are larger than desired. When we evaluated the synergy activations found by VO for the two legs, we identified four pairs that varied significantly in amplitude between the left and right sides. Based on this observation, we formulated a DO problem that sought to bring these four synergy activations with disparate amplitudes in the two legs closer to each other. Each pair of asymmetric synergy activations contained a “weak” and “strong” side. First, scaled-up “weak” synergy activations were tracked for each pair to increase their amplitudes. Second, a custom user-defined cost function was written to reward synergy amplitude similarities between the “weak” and “strong” synergy activations. A user-defined static parameter was added for each of the four synergy pairs representing the amplitude scale factor to be applied to the “strong’’ synergy activation using a user-defined model modification function. These two cost function terms increased “weak” synergy activations and encouraged over-active “strong” synergy activations to reduce to more symmetric magnitudes while not overconstraining the control space. In addition, a cost function term was added to find a solution with a metabolic cost per unit time that matched the value calculated for the subject’s original 0.5 m/s self-selected walking speed. This value was 90% of the original value calculated for the subject’s 0.8 m/s self-selected walking speed, making a 10% reduction the target. Cost function terms also tracked VO-found synergy controls for the two symmetric control pairs and VO-found torque controls. Constraint terms included bounded kinetic consistency errors (0.01 Nm) for all controlled joints, root segment residual loads for the pelvis (1 N for forces and 0.1 Nm for moments), and joint position periodicity errors (0.05 rad for joints below the pelvis and 0.1 rad for joints above it). The VO results were used in both the initial guess and tracked quantities data directories, and the final time was fixed at the value determined experimentally for the selected 0.8 m/s walking trial.

The DO predictive simulation produced modified kinematics, synergy activations, and ground reactions that reduced the subject’s metabolic cost per unit time by 9.9% for the selected gait cycle and achieved closely matching maximum amplitudes for the four modified synergy activations. Final results are plotted in Fig. 11. Convergence was achieved after 489 IPOPT iterations, which required approximately 730 min of CPU time. These DO results suggest that upregulating weak paretic leg synergy activations while downregulating paired non-paretic leg synergy activations could potentially allow the modeled subject to increase his self-selected walking speed by 60% without increasing his metabolic cost per unit time.

Fig. 11

Plots showing the final predictive simulation results from the DO tool. Subfigures A) kinematics, B) synergy activations, C) muscle activations, and D) ground reaction forces and moments after VO (blue) and DO (red), respectively

Discussion

The primary objective of this work was to develop a novel set of easy-to-use computational tools that facilitate the creation of personalized neuromusculoskeletal models and the subsequent design of personalized clinical treatments by researchers and clinicians with little to no programming experience. These tools dramatically reduce the barrier to entry for researchers and clinicians working together who would like to explore using personalized computational neuromusculoskeletal models to design personalized treatment approaches. As shown in the example application, the NMSM Pipeline has the potential to change the way treatments are investigated and applied throughout the neuromusculoskeletal injury space.

Going forward, two key questions need to be addressed if the NMSM Pipeline is to be used to design actual clinical treatments for actual patients. The first question is, how should the patient’s personalized model and neural control strategy be changed in response to the simulated treatment? The second question is, how can one validate a patient’s predicted post-treatment movement function? The answers to these two questions will depend on the specific clinical application and which aspects of the patient’s neuromusculoskeletal system need to be modeled and personalized to investigate that application. The answers will also involve developing an understanding of how the planned intervention will affect a patient’s model parameter values and neural control strategy. Collection of post-intervention movement data will be required so that the entire Model Personalization and Treatment Optimization process can be evaluated, refined, and ultimately validated for one clinical application at a time.

To conceptualize how these two important questions could be addressed, one can consider two NMSM Pipeline clinical applications being explored by the authors—one involving knee osteoarthritis rehabilitation design, and one involving pelvic cancer surgical planning. Past computational research involving a subject with bilateral medial compartment knee osteoarthritis investigated whether personalized gait modifications could be designed to reduce the subject’s peak adduction moment in both knees to an extent comparable to invasive high tibial osteotomy surgery [137]. For this treatment optimization problem, no muscles or neural control needed to be modeled, which simplified the treatment design process. The computational treatment design process involved: a) creating a personalized full-body dynamic skeletal model of the subject using pre-rehabilitation gait data, b) selecting a post-treatment control approach to be used when predicting post-treatment walking function (in this case, minimize changes in lower body joint moments relative to pre-rehabilitation gait data), c) performing a treatment optimization process with the personalized model and selected control approach to identify gait modifications that would offload the medial compartment of both knees, d) training the subject to perform and consolidate the modified walking motion, and finally e) re-testing the subject in a gait lab to evaluate how well the predicted gait modifications could reduce both peak adduction moment peaks. After informal training, the subject was able to perform the predicted gait modifications and achieve bilateral reductions in the peak knee adduction moment that were comparable to the treatment optimization predictions. A similar computational treatment design process is currently being repeated using the NMSM Pipeline except with greater flexibility in the range of potential gait modifications that can be predicted.

Ongoing computational research involving subjects with pelvic cancer is investigating whether surgical decisions can be planned to achieve a post-surgery gait motion and self-selected walking speed that is as close to normal as possible for each patient [68, 138]. For this treatment optimization problem, both muscles and neural control need to be modeled, which complicates the treatment design process. The computational treatment design process being explored involves: a) creating a personalized full-body dynamic neuromusculoskeletal model of the patient using pre-surgery gait, EMG, and imaging data (to personalize pelvis bony geometry), b) selecting a post-surgery control approach to be used when predicting post-surgery walking function (in this case, minimize changes in muscle synergy quantities relative to pre-surgery EMG data) [138], c) performing a treatment optimization process with the personalized model and selected control approach to identify surgical decisions that achieve clear margins while also maximizing recovery of normal walking function [68], d) having the orthopedic oncologist perform the patient’s surgery using the optimized surgical plan, and finally e) re-testing the patient in a gait lab following plateau in recovery to evaluate how well the predicted surgical plan allowed the patient to achieve the predicted post-surgery gait motion and speed. Since trial and error treatment is not possible with orthopedic surgeries, extensive testing, evaluation, and refinement of this treatment optimization process must be performed before actual surgical plans can be designed for actual pelvic cancer patients. To this end, the authors are collecting extensive pre- and post-surgery gait, EMG, and imaging data sets from this patient population, along with the implemented surgical decisions. Only after the NMSM Pipeline can predict post-surgery walking function reliably for multiple patients when post-surgery walking function and the implemented surgical decisions are known can use of the NMSM Pipeline be considered for this clinical application. A similar process would need to be followed for any other clinical problem to which the NMSM Pipeline would be applied.

Although the decision to use MATLAB for the NMSM Pipeline facilitates development and ease of access for users, the decision creates several limitations. C++ and other similar compiled programming languages generally run significantly faster than interpreted programming languages like MATLAB. However, they also introduce challenges since few engineering researchers are strong C++ programmers, and system-specific build procedures can often be complicated [139, 140]. Due to these drawbacks, most biomechanics software written in compiled languages can be modified only by reporting a bug or requesting a feature and waiting for experienced software developers to release a new version of the software. The NMSM Pipeline actively encourages modifications, bug fixes, and improvements through open-source easy-to-modify MATLAB toolsets. Other interpreted programming languages like Python can struggle with package management, virtual environments, and/or dependency version mismatches that increase barriers to use. MATLAB’s native debugging tools, English-readable error messages, and robust, stable core functions allow for development with confidence, making the NMSM Pipeline robust to issues common to other interpreted languages.

Our hope is that the NMSM Pipeline’s open-source, easy-to-use framework will encourage investigations by the research community into “best practices” for each of the included tools. With the availability of seven new tools, each with many optional features, the research community is encouraged to explore novel approaches for achieving the best results with each tool. Investigations into changes in task order, data processing/organization, cost function error centers and maximum allowable errors, among others, will inform the entire research community on the best ways to leverage these unique tools for successful future research and clinical outcomes.

We are already planning future enhancements to improve the versatility and applicability of the NMSM Pipeline to a broader range of clinical and research problems. Potential enhancements would allow greater flexibility and customization during use of the Model Personalization toolset. Support for compliant tendon muscle models and muscle models that differentiate between fast and slow twitch fibers in the MTP tool could facilitate muscle–tendon model personalization and simulation for “fast” activities like running or jumping [61, 111]. Support for time-delayed sensory feedback models in the NCP tool could enable investigation of how spasticity, hyperreflexia, and hypertonia interact with incoordination and weakness to produce functional impairment in individuals with an upper motor neuron pathology. Support for multi-segment foot models in the GCP tool could facilitate higher fidelity simulations of foot biomechanics and better reproduction of ground reaction moments. Other potential enhancements could facilitate novel use of the Treatment Optimization toolset. User-defined controls, time-delayed sensory feedback controls, and individual muscle controls could facilitate computational design of personalized functional electrical stimulation protocols [67] and prediction of how downregulating overactive reflexes could improve an individual’s movement function. Support for user-defined loads as a function of model states and controls and making it easy to add new states and controls to a DO problem could facilitate computational design of personalized passive and/or active exoskeletons, which would be added to a subject’s personalized model after TO and VO had been performed. An enhancement supporting multi-phase optimal control problems, where each phase simulates a different movement cycle with control information shared across phases, could facilitate simulating neural feedback mechanisms, neural adaptation, and cycle-to-cycle motion and control variability. Enhancements supporting the use of closed-chain kinematic models would facilitate simulation of shoulder models with improved kinematic fidelity for surgical, exoskeleton, and exercise equipment applications. Support for connecting NMSM Pipeline controls to any type of OpenSim actuator could allow exploration of alternative optimal control problem formulations that significantly improve IPOPT convergence rate, especially for TO. Finally, for both toolsets, development of easy-to-use MATLAB utility functions that automatically process experimental movement data and put it into a format that can be used directly by NMSM Pipeline tools, without requiring any manual data manipulation, would make the toolsets accessible to a wider user base.

With all of these potential future enhancements, several limitations and challenges in future NMSM Pipeline development exist. First, the quality of the personalization results found by the Model Personalization tools is limited by the quality and quantity of collected experimental data. At this time, it can be challenging to collect all of the experimental data needed to exercise all seven NMSM Pipeline tools with a synergy-driven neuromusculoskeletal model. Leveraging markerless motion capture or inertial measurement units could potentially increase the ability of researchers to collect the movement data required to perform a high-quality NMSM Pipeline run-through. For Model Personalization, because gradient-based optimization techniques are used, solutions may become entrapped in a local minimum, though the final solution will still be better than the initial guess. The use of muscle activations rather than excitations for the NCP tool and the Treatment Optimization toolset is another limitation, though back calculation of muscle excitations from muscle activations could potentially be performed given estimates of activation dynamics parameter values. Using the OpenSim framework through the OpenSim MATLAB API prevents the NMSM Pipeline from using automatic differentiation to speed up optimizations performed during the Model Personalization or Treatment Optimization process. Computations leveraging the OpenSim MATLAB API are inefficient, and speed limitations inherent to this framework cannot be overcome.

Several research teams have released novel software tools for the musculoskeletal modeling research community [51,52,53,54,55,56,57,58], and these tools could be used synergistically with the NMSM Pipeline. Though none of these tools possess the broad range of model personalization and treatment optimization capabilities available in the NMSM Pipeline (see Table 7), all of these tools facilitate some aspect of creating scaled musculoskeletal models that incorporate personalized bone geometry and muscle attachment locations. Consequently, these tools can provide the scaled musculoskeletal model required as a starting point for the NMSM Pipeline Model Personalization process. For example, NMSBuilder, MAP Client, Torsion Tool, Bone Deformation Tool, or OpenSim Creator could be used to create personalized bone models from imaging data or physical measurements [52,53,54,55, 57, 141, 142]. SimCP, NMSBuilder, MAP Client, or OpenSim Creator could then be used to personalize muscle attachment locations [51,52,53, 57]. Next, SimCP, NMSBuilder, MAP Client, AddBiomechanics, or Scale Tool could be used to create a scaled OpenSim musculoskeletal model possessing the personalized bone geometry and muscle attachments [51,52,53, 56]. In addition, AddBiomechanics could potentially be used to perform limited initial personalization of joint functional axes [56], while SimCP could be used to perform initial personalization of optimal muscle fiber lengths, tendon slack lengths, and synergy control properties [51].

Table 7 Comparison of functionality available in the NMSM Pipeline and other musculoskeletal modeling software tools

As shown in Table 7 above, a current limitation of the NMSM Pipeline is the lack of functionality to personalize a patient’s bone morphology to clinical imaging data. The NMSM Pipeline was specifically designed to avoid the use of imaging data and require only experimental data typically collected in a human movement lab. However, personalized bone geometry generated using existing methods and software tools [53, 138, 143] can easily be incorporated into scaled generic OpenSim models with personalized joint functional axes [138] to produce more accurate muscle-tendon lengths, moment arms, and joint and bone loading scenarios.

Treatment Optimization also does not currently include a muscle fatigue model. To date, few muscle fatigue models for neuromusculoskeletal modeling applications have been published in the literature [144, 145]. Thus, selecting an appropriate existing muscle fatigue model would be challenging. Availability of a muscle fatigue model could improve the utility of predictive simulations for athletic, other high-intensity, or long duration movements, making this research area an important one for future investigation.

The NMSM Pipeline’s Treatment Optimization toolset has key differences in design and implementation from the existing OpenSim Moco tool. While Moco is generally the faster software due to its C++-based implementation, Treatment Optimization leverages the flexibility of MATLAB and GPOPS-II to facilitate novel features. Current unique features include user-defined functions and the use of .osimx personalized model components, while future unique features could include support for multi-phase optimal control problems, as discussed above. Treatment Optimization was designed to be accessible to both novice and advanced users through intuitive XML settings files, whereas Moco users typically set up problems through MATLAB or Python scripts. Even so, Moco’s implementation facilitates potential use of C++ automatic differentiation libraries to improve computational speed for gradient calculations, resulting in more efficient optimizations. Lastly, Moco is free while both MATLAB and GPOPS-II are commercial products with significant academic discounts. By incorporating NMSM Pipeline model components into native OpenSim in the future, we plan to give users the option of using their NMSM Pipeline personalized models in either our Treatment Optimization toolset or OpenSim Moco, allowing the research community to explore the tradeoffs between these two approaches for developing predictive simulations.

In conclusion, the NMSM Pipeline represents a significant step forward in the ability to personalize neuromusculoskeletal computer models to patient movement data and then optimize parameters of potential treatments using the personalized models. The software assembles over two decades of biomechanics research into a well-designed and accessible framework that makes previously “high end” modeling, simulation, and optimization techniques accessible to the entire research community. Our hope is that by empowering researchers and clinicians working together to explore clinically relevant questions and to design novel personalized interventions, the NMSM Pipeline will bring neuromusculoskeletal computer modeling not only to the doorstep of clinical utility but also across it.