Use the table of contensts to easily navigate this README. Click on the three lines next to README.md just above this sentence.
Predictive Simulations of Human Movement
This repository contains code and data to generate three-dimensional muscle-driven predictive simulations of human gait. The simulation framework is developed by Falisse et al.. The implementation in this repo is aimed at letting you run simulations with your customized musculoskeletal models*.
If you want to reproduce specific published results, you are adviced to visit the repo corresponding to the paper.
This repository is a work in progress
- If you encounter an error or get an unexpected output, check the issues. Please look at the existing issues before creating a new one, you might not be the first person to have this problem. Include the error messages in the issue.
- Feel free to suggest improvements. Submit a pull request with the changes, or create an issue with a desciption of the proposed change.
Citations
Please cite the original authors if you use this simulation framework.
-
Falisse A, Serrancoli G, Dembia C, Gillis J, Jonkers J, De Groote F. 2019 Rapid predictive simulations with complex musculoskeletal models suggest that diverse healthy and pathological human gaits can emerge from similar control strategies. Journal of the Royal Society Interface 16: 20190402. http://dx.doi.org/10.1098/rsif.2019.0402. You can find the original repository here: https://github.com/antoinefalisse/3dpredictsim (MATLAB)
-
Falisse, A., Serrancolí, G., Dembia, C. L., Gillis, J., & De Groote, F. 2019 Algorithmic differentiation improves the computational efficiency of OpenSim-based trajectory optimization of human movement. PLOS ONE, 14(10), e0217730. https://doi.org/10.1371/journal.pone.0217730
-
Falisse A, Afschrift M, De Groote F. 2022 Modeling toes contributes to realistic stance knee mechanics in three-dimensional predictive simulations of walking. PLoS ONE 17(1): e0256311. https://doi.org/10.1371/journal.pone.0256311. You can find the original repository here: https://github.com/antoinefalisse/3dpredictsim_mtp (python)
-
D’Hondt, L., De Groote, F., & Afschrift, M. 2024 A dynamic foot model for predictive simulations of human gait reveals causal relations between foot structure and whole-body mechanics. PLOS Computational Biology, 20(6), e1012219. https://doi.org/10.1371/journal.pcbi.1012219. You can find the original repository here: https://github.com/Lars-DHondt-KUL/3dpredictsim/tree/four-segment_foot_model (MATLAB)
If you are familiar with a previous version of this simulation framework, here is a rundown of the changes
In general, hard coded variables were removed as much as possible and put in a settings structure the user can adapt. This results in:
- More easily use a different model (other number of joints or muscles)
- Make input-output a bit easier
- More easily change the (musculo-skeletal) parameters of the simulation
Besides that, the code was also adapted in accordance with recent (October 2020) adjustments in the original repository of Falisse et al.
- Adjusted version of the collocation scheme (now really being an orthogonal radau scheme)
- Support for parallel computing
- Formulation with opti
Lastly, seperate pieces of code were put together to streamline performing predictive simulations:
- Automatic conversion of an OpenSim model to the external function (executable called from the workflow)
- Perorming Muscle Analysis and polynomial fitting (integrated into the workflow)
Required software
To run this code you need to have the following softwares on your machine:
How to setup the code
- Fork this repository to your github account. If you want the fork to be private, follow these steps instead.
- Clone the fork to your machine. If you have a computer with restricted permissions, make sure you have permission to run executables from the selected folder (For computers with KU Leuven BioMed Group policies, this is C:\GBW_MyPrograms\ ).
Do not download the code as zip.
- Get the OpenSim API running on MATLAB. See Setting up your Matlab Scripting Environment.
- In main.m, change S.solver.CasADi_path to reflect the location where you installed CasADi.
- Make sure the opensimAD submodule is installed. If PredSim\opensimAD\ is empty, open git command prompt, go to ...\PredSim\ , and run
git submodule update --init
.
After perfoming these steps, run the main script. (Expected run time is 40 minutes, depending on hardware.) If you don't receive any errors, your results should be the same as https://github.com/KULeuvenNeuromechanics/PredSim/tree/master/Tests/Falisse_et_al_2022_Results. If that is the case, you have succesfully intalled and set up the code. You are ready to do your own simulations.
How to use the code
The code is written such that as a user you only have to interact with main.m.
All user-defined settings are stored in structure S. In main.m you have to specify the required settings and are free to change/add the optional settings.
Before running a simulations
This code can automatically convert an OpenSim model to the external function used in the simulations. This greatly simplifies the process of going from a subject-specific model to a predictive simulation. Nevertheless, you should take care of the model you use since not all OpenSim models are suported:
- Model should be 3D.
- Your model should not have locked joints. Locked joints would technically require having kinematic constraints, which is possible but makes the problem more complicated. Replace them with weld joints instead.
- Constraints on coordinates will be ignored (eg, coupling constraints).
- Using SimmSplines to describe coordinates (e.g. Yamaguchi knee model) is not supported as the implementation in OpenSim is not really compatible with algorithmic differentiation. Change them to Polynomials instead. GeometryPaths can contain SimmSplines. AdaptOpenSimModel.m takes care of changing present SimmSplines to polynomials.
- Your model needs to have contact elements that interact with the ground. Only SmoothSphereHalfSpaceForce contact forces are supported. You can use AdaptOpenSimModel.m to add contact geometries and forces to your model. You can also scale the radius, stiffness and dissipation of the contact spheres.
- Your model can have any Hill-type muscle model, but it will be implemented as a DeGroote-Fregly muscle.
- Torque/force actuators of the class ActivationCoordinateActuator are supported. You can add actuators by running AdaptOpenSimModel.m. Actuators are not required.
Required Settings
- S.subject.name:
- the name or code of the subject you are simulating.
- osim_path:
- path to the scaled opensim model of the subject.
- S.subject.save_folder:
- path to the folder where you want to store the results of the OCP. If the folder does not exist yet on your machine, it will be created automatically.
- S.subject.IG_selection:
- either choose 'quasi-random' or give the path to a .mot file you want to use as initial guess.
- S.subject.IG_selection_gaitCyclePercent:
- if S.subject.IG_selection is a .mot file, S.subject.IG_selection_gaitCyclePercent is required. Here, specify what percent of gait cycle does the .mot file contain. For example, if the .mot file has 2 gait cycles, S.subject.IG_selection_gaitCyclePercent is 200.
- S.solver.run_as_batch_job:
OptionalSettings
S.bounds
S.metabolicE - metabolic energy
- S.metabolicE.tanh_b:
- hyperbolic tangeant smoothing factor used in the metabolic cost calculation. Default is 10 [double]
- S.metabolicE.model:
- the name of the metabolic energy model used. Default is Bhargava2004 [char]. Currently only Bhargava2004 model has been implemented. Other options that will be added in the future are:
S.misc - miscellanious
-
S.misc.gaitmotion_type:
- type of gait simulation. Default is HalfGaitCycle [char]. Other option is FullGaitCycle [char]. Simulating a half gait cycle reduces computation time, but is limited to symmetric models. Post-processing will always reconstruct a full gait cycle starting at right heel strike.
-
S.misc.msk_geom_eq:
- type of equation to approximate musculo-skeletal geometry (moment arm and muscle-tendon lengths wrt. joint angle). Default is polynomials [char]
-
S.misc.threshold_lMT_fit:
- Threshold RMSE on muscle-tendon length to accept the polynomial fit. Default is 0.003 m [double]
-
S.misc.threshold_dM_fit:
- Threshold RMSE on muscle-tendon moment arm to accept the polynomial fit. Default is 0.003 m [double]
-
S.misc.poly_order.lower:
- minimal order of polynomial function. Default is 3 [double]
-
S.misc.poly_order.upper:
- maximal order of polynomial function. Default is 9 [double]
-
S.misc.default_msk_geom_bound:
- file with default values for upper and lower bounds for approximating musculoskeletal geometry. Rotations are assumed in degrees, translations in meters. Default is 'default_msk_geom_bounds.csv' [char].
The provided file should be compatible with
readtable
. The table should contain a column with coordinate names (header: name), a column with lower bounds, in degrees and meters (header: lower), and a column with upper bounds, in degrees and meters (header: upper). To set only the upper or lower bound of a coordinate, set the other one to nan
;
-
S.misc.msk_geom_bounds:
- Cell array where 1st entry is dof name(s) , 2nd entry is its lower bounds, and 3rd entry is its upper bounds. Insert nan to lower bounds to only overwrite upper bounds. For another bound, add 3 more entries. For example, {{'knee_angle_r','knee_angle_l'},-120,10,'lumbar_extension',nan,30} implements limit of -120° and 10° on knee angles, and default lower bound with 30° upper bound for lumbar_extension.
Order of priority for bounds:
- Individual bounds from settings (S.misc.msk_geom_bounds)
- Default bounds from table (S.misc.default_msk_geom_bound)
- Read from model file. Qs: min and max coordinate values
-
S.misc.msk_geom_n_samples:
- Number of samples for the dummy motion that is used to fit the approximated musculoskeletal geometry. Default is 5000 [double]
-
S.misc.visualize_bounds:
- specify if bounds and initial guess are visualized (0 or 1). Default is 0 [double]
-
S.misc.dampingCoefficient:
- damping coefficient of muscles. Default is 0.01 [double]. Used as damping value that is multiplied by the normalized muscle velocity, in the muscle velocity dependent term in calculation of normalized contractile element force of the muscle.
-
S.misc.constant_pennation_angle:
- specify if pennation angle of the muscles is supposed to stay constant (0 or 1). Default is 0 [double]
-
S.misc.default_scaling_NLP:
- Filename with table that contains scale factors for Qs, Qdots, Qdotdots, and Moments in its columns. The first column should contain the coordinate names for its corresponding row. Default is '', i.e. scale factors are derived from bounds. [char]
-
S.misc.scaling_Qs:
- cell array of name-value pairs of coordinate names and the scale factor for the optimisation variables corresponding to their position.
-
S.misc.scaling_Qdots:
- cell array of name-value pairs of coordinate names and the scale factor for the optimisation variables corresponding to their velocity.
-
S.misc.scaling_Qdotdots:
- cell array of name-value pairs of coordinate names and the scale factor for the optimisation variables corresponding to their acceleration.
-
S.misc.scaling_moments:
- cell array of name-value pairs of coordinate names and the scale factor for the constraint violation on their moment equilibrium.
-
S.misc.git.local_hash:
- hash of the local instance [char]. This is the identifier of the version of the code on your machine. You cannot change this setting.
-
S.misc.git.branch_name:
- current branch of the local instance [char]. You cannot change this setting.
-
S.misc.git.remote_hash:
- hash of the last commit on the remote [char]. This is the identifier of the latest version on the remote, i.e. GitHub. You cannot change this setting.
-
S.misc.computername:
- name of the computer on which the simulation was run [char]. You cannot change this setting.
S.post_process
- S.post_process.make_plot:
- boolean to plot post processing results (0 or 1). Default is 0.
- S.post_process.rerun:
- boolean to rerun post-processing without solving OCP (0 or 1). Default is 0. If this option is set to 1, one should specify the S.post_process.result_filename.
- S.post_process.result_filename:
- File name for results. Used for the name of .mat file that saves the results, diary of the OCP, and name of the .mot file of the output motion. When rerunning post-processing of an existing result, giving this file name is required. Default value is
- S.post_process.savename:
- Type of savename to use if S.post_process.result_filename is empty. Defaults is structured [char]. This uses the name of the .mat file of results is used as <S.subject.name>_v\. Where <S.subject.name> is defined in S.subject.name. n = 1 if <S.subject.name>_v1.mat does not exist. n is increased until n is found such that <S.subject.name>_v\.mat does not exist. To change this structuring process, change its implementation in run_pred_sim.m file. An alternative option is datetime [char], this uses <S.subject.name>_\. Where \ is the system date and time when creating the savename.
S.solver
- S.solver.linear_solver:
- solver algorithm used for the OCP. Default is mumps [char].
- S.solver.tol_ipopt:
- the power (10^-x) the tolerance ipopt has to reach before the OCP can be regarded as solved; a higher number gives a more precise answer. Default is 4 [double]
- S.solver.max_iter:
- maximal amount of iterations after wich the solver will stop. Default is 10000 [double]
- S.solver.parallel_mode:
- type of parallel computing. Default is thread [char].
- S.solver.N_threads:
- number of threads in parallel mode. Default is 4 [double]. When using batch computing, this value is overwritten with the number of threads assigned to each worker in you parallel cluster.
- S.solver.N_meshes:
- number of mesh intervals. Default is 50 [double] for S.misc.gaitmotion_type = HalfGaitCycle and 100 for FullGaitCycle
S.subject
- S.subject.save_folder:
- folder path to store the intermediate subject specific results (e.g. external function with skeletal dynamics, CasADi function with musculoskeletal geometry polynomials). If the folder does not exist, it is created automatically.
- S.subject.mass:
- mass of the subject in kilograms. Default is [] kilograms [double]. Default is empty, it will be overwritten by the mass extracted from the OpenSim model.
- S.subject.IG_pelvis_y:
- height from ground to pelvis, in meters. Default is [] m [double]. Default is empty, it will be overwritten by pelvis height extracted from the OpenSim model.
- always used for the quasi-random initial guess
- used for data-informed initial guess when
S.subject.adapt_IG_pelvis_y = 1;
- S.subject.IG_pelvis_y is also used to establish bounds on vertical pelvis position.
- S.subject.adapt_IG_pelvis_y:
- boolean to adjust the trajectory of height of pelvis from the ground for data-informed initial guess. Default is 0. 0 means the trajectory will not be changed. If 1, the trajectory will be changed such that the average value of the trajectory is equal to s.subject.IG_pelvis_y.
- S.subject.v_pelvis_x_trgt:
- average velocity you want the model to have, in meters per second. Default is 1.25 m/s [double]
- S.subject.muscle_strength:
- structure with scaling factors for muscle strength. This scales the max muscle force of the active muscle force. Default is [], that is, no scaling. Input as a cell array where 1st input is the muscle(s) name, 2nd is the scale factor. If more than one scaling is to be performed, add 2 more inputs. For example, S.subject.muscle_strength = {{'soleus_l','soleus_r'},0.9,{'tib_ant_l'},1.1} will scale both soleus by a factor of 0.9 and tibialis anterior left by a scale of 1.1.
- S.subject.muscle_pass_stiff_scale:
- structure with scaling factors for muscle passive stiffness. Default is [], that is, no scaling. Input as a cell array where 1st input is the muscle(s) name, 2nd is the scale factor. If more than one scaling is to be performed, add 2 more inputs. For example, S.subject.muscle_pass_stiff_scale = {{'soleus_l','soleus_r'},0.9,{'tib_ant_l'},1.1} will scale both soleus by a factor of 0.9 and tibialis anterior left by a scale of 1.1.
- S.subject.muscle_pass_stiff_shift:
- structure with scaling factors for muscle passive stiffness shift. This property shifts the start of passive muscle force from normalized muscle length = 1 to the valuse specified in S.subject.muscle_pass_stiff_shift. Default is [], that is, no scaling. Input as a cell array where 1st input is the muscle(s) name, 2nd is the scale factore. If more than one scaling is to be performed, add 2 more inputs. For example, S.subject.muscle_pass_stiff_shift = {{'soleus_l','soleus_r'},0.9,{'tib_ant_l'},1.1} will scale both soleus by a factor of 0.9 and tibialis anterior left by a scale of 1.1.
- S.subject.tendon_stiff_scale:
- structure with scaling factors for tendon stiffnesses. Default is [], that is, no scaling. Input as a cell array where 1st input is the muscle(s) name, 2nd is the scale factor. If more than one scaling is to be performed, add 2 more inputs. For example, S.subject.tendon_stiff_scale = {{'soleus_l','soleus_r'},0.9,{'tib_ant_l'},1.1} will scale both soleus by a factor of 0.9 and tibialis anterior left by a scale of 1.1.
- S.subject.mtp_type:
- type of mtp joint. Default is '' [char], which treats the mtp like any other joint. Select '2022paper' to use passive mtp joints whose kinematics do affect the crossing muscle-tendon units (Falisse et al., 2022).
- S.subject.scale_MT_params:
- scale muscle tendon properties that are read from opensim model. Default is [], that is, no scaling. Input as a cell array where 1st input is the muscle(s) name, 2nd is what property you want to scale (FMo, lMo, lTs, alphao or vMmax), 3rd is the scale factor itself. If more than one scaling is to be performed, add 3 more inputs. For example, S.subject.scale_MT_params = {{'soleus_l','soleus_r'},'FMo',0.9,{'tib_ant_l'},'lTs',1.1} will scale max isometric force of both soleus by a factor of 0.9 and tendon slack length of tibialis anterior left by a scale of 1.1.
- S.subject.damping_coefficient_all_dofs:
- damping coefficient for all coordinates (except coordinates connected to ground, generally pelvis (also called floating base)). Default is 0.1 Nms/rad [double]
- S.subject.set_damping_coefficient_selected_dofs:
- damping coefficient can be specified here for each coordinate individually. For example, S.subject.set_damping_coefficient_selected_dofs = {{'hip_flexion_l','hip_flexion_r'},0.12,{'knee_angle_l'},0.11} will put damping coefficient of both hip flexions to 0.12 and that of knee angle left to 0.11. If not defined here for a particular coordinateS.subject.damping_coefficient_all_dofs will be used for that coordinate. Default is empty.
- S.subject.stiffness_coefficient_all_dofs:
- stiffness coefficient for all coordinates (except coordinates connected to ground, generally pelvis (also called floating base)). Default in 0 Nm/rad [double]
- S.subject.set_stiffness_coefficient_selected_dofs:
- stiffness coefficient can be specified here for each coordinate individually. For example, S.subject.set_stiffness_coefficient_selected_dofs = {{'hip_flexion_l','hip_flexion_r'},0.012,{'knee_angle_l'},0.011} will put stiffness coefficient of both hip flexions to 0.012 Nm/rad and that of knee angle left to 0.011 Nm/rad. If not defined here for a particular coordinate, S.subject.stiffness_coefficient_all_dofs will be used for that coordinate. Default is empty.
- S.subject.set_stiffness_offset_selected_dofs:
- position where moment of linear stiffness is zero can be specified here for each coordinate individually. For example, S.subject.set_stiffness_offset_selected_dofs = {{'hip_flexion_l','hip_flexion_r'},-0.3,} will offset the hip flexion stiffness such that their moment is zero at -0.3 rad hip flexion. Default is empty.
- S.subject.default_coord_lim_torq_coeff:
- file with default coefficients for coordinate limit torques. Default is 'default_coord_lim_torq_coeff.csv' [char].
The provided file should be compatible with
readtable
. The table should contain a column with coordinate names (header: name), 4 columns with stiffness coefficients, (headers: K_1, K_2, K_3, K_4), and 2 columns with offset coefficients (headers: theta_1, theta_2).
Limit torques are calculated in function of coordinate value q as: Tau = K(1)*exp(K(2)*(q-theta(2))) + K(3)*exp(K(4)*(q-theta(1)))
.
Default coefficients are taken from Anderson III, Frank Clayton. A dynamic optimization solution for a complete cycle of normal gait. The University of Texas at Austin, 1999.
- S.subject.scale_default_coord_lim_torq:
- scale factor for the amplitude of all default coordinate limit torques. Default is empty [double].
All
K(1)
and K(3)
values read from the file with default coefficients are multiplied by this factor.
To scale the limit torques for individual degrees of freedom, use S.subject.set_limit_torque_coefficients_selected_dofs.
- S.subject.set_limit_torque_coefficients_selected_dofs:
- Set limit torque coefficients for a coordinate. Default is empty [cell array] with pattern {coordinate name(s) [char, cell array of chars], K [4x1 double], theta [2x1 double]}.
For the specified coordinates, the coefficients from this setting take priority over the coefficients from the file with defaults (S.subject.default_coord_lim_torq_coeff).
Note: Setting
K(1) = nan
will cause the limit torque for that coordinate to be excluded, even if it was given in the defaults.
- S.subject.base_joints_legs:
- Joint name that is the base of a leg, left and right. Default is 'hip' [char]. Inputs of the form 'hip_r', {'hip_l'}, {'hip_r','hip_l'} are equivalent.
- S.subject.base_joints_arms:
- Joint name that is the base of an arm, left and right. Default is 'acromial' [char]. Inputs of the form 'acromial_r', {'acromial_l'}, {'acromial_r','acromial_l'} are equivalent. Set to empty [] if the model does not have arms.
- S.subject.stiffness_all_ligaments:
- Default stiffness model (i.e. force-length) used for ligaments. Default is ligamentGefen2002 [char]
- S.subject.set_stiffness_selected_ligaments:
- Use name-value pairs to use different stiffness models for specific ligaments. Default is {'PlantarFascia','plantarFasciaNatali2010'}
S.weights
- S.weights.E:
- weight on metabolic energy rate. Default is 500 [double]
- S.weights.E_exp:
- exponent for the metabolic energy rate. Default is 2 [double]
- S.weights.q_dotdot:
- weight on joint accelerations. Default is 50000 [double]
- S.weights.e_arm:
- weight on arm excitations. Default is 10^6 [double]
- S.weights.pass_torq:
- weight on passive torques. Default is 1000 [double]
- S.weights.pass_torq_includes_damping:
- specify if damping torque = damping coefficient * coordinate velocity is to be included in the cost function (0 or 1). Default is 0 [double].
- S.weights.a:
- weight on muscle activations. Default is 2000 [double]
- S.weights.slack_ctrl:
- weight on slack controls. Default is 0.001 [double]
S.OpenSimADOptions
These settings are passed to OpenSimAD.
- S.OpenSimADOptions.compiler:
- command prompt argument for the compiler. [char]
By default, PredSim will look for the most recent version that is installed in either
C:/Program Files/Microsoft Visual Studio/
or C:/Program Files (x86)/Microsoft Visual Studio/
.
If you get an error about not finding a compiler, use this setting to specify your compiler:
- Visual studio 2015: 'Visual Studio 14 2015 Win64'
- Visual studio 2017: 'Visual Studio 15 2017 Win64'
- Visual studio 2017: 'Visual Studio 16 2019'
- Visual studio 2017: 'Visual Studio 17 2022'
- S.OpenSimADOptions.verbose_mode:
- print outputs from windows command prompt to matlab command window (and log file). Default is true [bool].
- S.OpenSimADOptions.verify_ID:
- verify the generated function versus the inverse dynamics tool in OpenSim. Default is false [bool].
- S.OpenSimADOptions.jointsOrder:
- If you want to choose the order of the joints outputs. Default is empty, which uses the joint order of the .osim file. [cell array of char]
- S.OpenSimADOptions.coordinatesOrder:
- If you want to choose the order of the coordinate outputs. Default is empty, which uses the coordinate order of the .osim file. [cell array of char]
S.OpenSimADOptions.jointsOrder and S.OpenSimADOptions.coordinatesOrder are included in the settings to aid backward compatibility.
- S.OpenSimADOptions.input3DBodyForces:
- add 3D force vectors that act on bodies. Default is empty. Needs further implementations before this can be used.
- S.OpenSimADOptions.input3DBodyMoments:
- add 3D moment vectors that act on bodies. Default is empty. Needs further implementations before this can be used.
- S.OpenSimADOptions.export3DPositions:
- export 3D position of points in bodies, in ground reference frame. Default is empty. Needs further implementations before this can be used.
- S.OpenSimADOptions.export3DVelocities:
- export 3D velocity of points in bodies, in ground reference frame. Default is empty. Needs further implementations before this can be used.
- S.OpenSimADOptions.exportGRFs:
- Export total ground reaction forces of left and right side. Default is true [bool]
- S.OpenSimADOptions.exportSeparateGRFs:
- Export ground reaction forces of each contact element. Default is true [bool]
- S.OpenSimADOptions.exportGRMs:
- Export total ground reaction moments of left and right side. Default is true [bool]
- S.OpenSimADOptions.exportContactPowers:
- Export power due to vertical compression of each contact element. Default is true [bool]