(close to) minimizing metabolic cost

[Nicholas Bianco]

Thanks for the details Carlos! This is very helpful.

It sounds like you're on the right track with your approach and are getting good preliminary results. To clarify, you are seeing high metabolic cost both normal and impaired gait?

Metabolic cost is tricky to predict accurately. If you're seeing the expected kinematic changes for pathological gait, you might not need to worry about metabolic cost so much (unless it is more important to your study than I'm understanding). That being said, I don't think it would be incorrect to add energy expenditure to the cost function, as it is common in predictive simulations of gait.

This paper by Carmichael Ong from our group, where he simulated pathological gait by weakening muscle groups, may be interesting to you: https://journals.plos.org/ploscompbiol/ ... bi.1006993.

[Carlos Gonçalves]

Great! Great to see that I'm on the right track. Thanks a lot.

I will check the paper. I reckon that I'm familiar with it, but since I decided to invest in Moco instead of SCONE (no hard feelings), I skip it. It will be of great help right now since it is also a 2D simulation.

I've managed to get something very close to the "heel walking pattern." But my colleagues (physical therapists) don't see this pattern in our patients. That's the main issue. What we were expecting was something close to a crouch gait.

Maybe I will focus on the contracture case.

My other hypothesis is that a crouch gait will only emerge if I use a 3D model. So the hip rotations will make it even harder to sustain the load during middle stance.

About the metabolic cost, I'm still implementing the Umberger2010MetabolicsProbe in my code. Soon I will be testing it with the predictions. The plan was to increase the control weight of the muscles that have the highest metabolic costs.

My goal is to test different AFO designs in the pathologic gait prediction and evaluate the results. We are trying to customize and optimize the AFO design for each patient in our Orthopedic Workshop.

Thanks again for the help.

Best regards.

[Aaron Fox]

Hi All,

Given the interest on this thread and there being some Python discussion I thought I'd chip in with this example. I gave the 0.5 early release Chris posted a go in Python adapting the 2Dgait C++ metabolics example currently on GitHub. Here is the Python adaptation I used, note that this doesn't have GRF tracking included and is pretty lazily thrown together in parts (i.e. code commenting is pretty lazy):

Code: Select all

# -*- coding: utf-8 -*-
"""
Created on Wed Oct  7 14:58:48 2020

@author:
    Aaron Fox
    Centre for Sport Research
    Deakin University
    aaron.f@deakin.edu.au
 
     Test of example2d walking with metabolics, based off:
         https://github.com/opensim-org/opensim-moco/blob/master/Moco/Examples/C++/example2DWalking/example2DWalkingMetabolics.cpp
 
"""

import opensim as osim
import math

# %%

#Setup tracking object
track = osim.MocoTrack()

#Get model
baseModel = osim.Model('2d_gait.osim')

#Setup metabolics
metabolics = osim.Bhargava2004Metabolics()
metabolics.setName('metabolics')
metabolics.set_use_smoothing(True)
metabolics.addMuscle('hamstrings_r',baseModel.getComponent('hamstrings_r'))
metabolics.addMuscle('hamstrings_l',baseModel.getComponent('hamstrings_l'))
metabolics.addMuscle('bifemsh_r',baseModel.getComponent('bifemsh_r'))
metabolics.addMuscle('bifemsh_l',baseModel.getComponent('bifemsh_l'))
metabolics.addMuscle('glut_max_r',baseModel.getComponent('glut_max_r'))
metabolics.addMuscle('glut_max_l',baseModel.getComponent('glut_max_l'))
metabolics.addMuscle('iliopsoas_r',baseModel.getComponent('iliopsoas_r'))
metabolics.addMuscle('iliopsoas_l',baseModel.getComponent('iliopsoas_l'))
metabolics.addMuscle('rect_fem_r',baseModel.getComponent('rect_fem_r'))
metabolics.addMuscle('rect_fem_l',baseModel.getComponent('rect_fem_l'))
metabolics.addMuscle('vasti_r',baseModel.getComponent('vasti_r'))
metabolics.addMuscle('vasti_l',baseModel.getComponent('vasti_l'))
metabolics.addMuscle('gastroc_r',baseModel.getComponent('gastroc_r'))
metabolics.addMuscle('gastroc_l',baseModel.getComponent('gastroc_l'))
metabolics.addMuscle('soleus_r',baseModel.getComponent('soleus_r'))
metabolics.addMuscle('soleus_l',baseModel.getComponent('soleus_l'))
metabolics.addMuscle('tib_ant_r',baseModel.getComponent('tib_ant_r'))
metabolics.addMuscle('tib_ant_l',baseModel.getComponent('tib_ant_l'))
baseModel.addComponent(metabolics)
baseModel.finalizeConnections()

#Process model
modelProcessor = osim.ModelProcessor(baseModel)

#Setup of MocoTrack
track.setModel(modelProcessor)
tabProcessor = osim.TableProcessor('referenceCoordinates.sto')
tabProcessor.append(osim.TabOpLowPassFilter(6))
track.setStatesReference(tabProcessor)
track.set_states_global_tracking_weight(30.0)
track.set_allow_unused_references(True)
track.set_track_reference_position_derivatives(True)
track.set_apply_tracked_states_to_guess(True)
track.set_initial_time(0.0)
track.set_final_time(0.47008941)
study = track.initialize()
problem = study.updProblem()

#Goals

#Symmetry
periodicityGoal = osim.MocoPeriodicityGoal('symmetryGoal')
problem.addGoal(periodicityGoal)
model = modelProcessor.process()
model.initSystem()

#Set symmetric coordinate values
#Set symmetry pairs
symPairs = [['hip_flexion_r','hip_flexion_l'],
            ['knee_angle_r','knee_angle_l'],
            ['ankle_angle_r','ankle_angle_l']]
for ii in range(0,len(symPairs)):
    #Set the jointsets depending on current pair
    if 'hip' in symPairs[ii][0]:
        jointName = ['/jointset/hip_r/','/jointset/hip_l/']
    elif 'knee' in symPairs[ii][0]:
        jointName = ['/jointset/knee_r/','/jointset/knee_l/']
    elif 'ankle' in symPairs[ii][0]:
        jointName = ['/jointset/ankle_r/','/jointset/ankle_l/']
    
    #Set the pair for coordinate value
    pair = osim.MocoPeriodicityGoalPair(jointName[0]+symPairs[ii][0]+'/value',
                                        jointName[1]+symPairs[ii][1]+'/value')
    #Add to the goal
    periodicityGoal.addStatePair(pair)
    
    #Set the pair for coordinate speed
    pair = osim.MocoPeriodicityGoalPair(jointName[0]+symPairs[ii][0]+'/speed',
                                        jointName[1]+symPairs[ii][1]+'/speed')
    #Add to the goal
    periodicityGoal.addStatePair(pair)

#Symmetric pelvis joint coordinates
periodicityGoal.addStatePair(osim.MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_ty/value'))
periodicityGoal.addStatePair(osim.MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_ty/speed'))

#Set symmetric muscle activations
#Set symmetry pairs
symPairs = [['hamstrings_r','hamstrings_l'],
            ['bifemsh_r','bifemsh_l'],
            ['glut_max_r','glut_max_l'],
            ['iliopsoas_r','iliopsoas_l'],
            ['rect_fem_r','rect_fem_l'],
            ['vasti_r','vasti_l'],
            ['gastroc_r','gastroc_l'],
            ['soleus_r','soleus_l'],
            ['tib_ant_r','tib_ant_l']]
for ii in range(0,len(symPairs)):
    
    #Set the pair for coordinate value
    pair = osim.MocoPeriodicityGoalPair('/'+symPairs[ii][0]+'/activation',
                                        '/'+symPairs[ii][1]+'/activation')
    #Add to the goal
    periodicityGoal.addStatePair(pair)

#Symmetric coordinate actuator controls
periodicityGoal.addControlPair(osim.MocoPeriodicityGoalPair('/lumbarAct'))

#Effort
effort = osim.MocoControlGoal.safeDownCast(problem.updGoal('control_effort'))
effort.setWeight(0.1)

#Metabolics
metGoal = osim.MocoOutputGoal('met',0.1)
metGoal.setOutputPath('/metabolics|total_metabolic_rate')
metGoal.setDivideByDisplacement(True)
metGoal.setDivideByMass(True)
problem.addGoal(metGoal)

#Bounds
problem.setStateInfo('/jointset/groundPelvis/pelvis_tilt/value', [-20*math.pi/180, -10*math.pi/180])
problem.setStateInfo('/jointset/groundPelvis/pelvis_tx/value', [0, 1])
problem.setStateInfo('/jointset/groundPelvis/pelvis_ty/value', [0.75, 1.25])
problem.setStateInfo('/jointset/hip_l/hip_flexion_l/value', [-10*math.pi/180, 60*math.pi/180])
problem.setStateInfo('/jointset/hip_r/hip_flexion_r/value', [-10*math.pi/180, 60*math.pi/180])
problem.setStateInfo('/jointset/knee_l/knee_angle_l/value', [-50*math.pi/180, 0])
problem.setStateInfo('/jointset/knee_r/knee_angle_r/value', [-50*math.pi/180, 0])
problem.setStateInfo('/jointset/ankle_l/ankle_angle_l/value', [-15*math.pi/180, 25*math.pi/180])
problem.setStateInfo('/jointset/ankle_r/ankle_angle_r/value', [-15*math.pi/180, 25*math.pi/180])
problem.setStateInfo('/jointset/lumbar/lumbar/value', [0, 20*math.pi/180])

#Configure the solver
solver = osim.MocoCasADiSolver.safeDownCast(study.updSolver())
solver.set_num_mesh_intervals(50)
solver.set_verbosity(2)
solver.set_optim_solver('ipopt')
solver.set_optim_convergence_tolerance(1e-4)
solver.set_optim_constraint_tolerance(1e-4)
solver.set_optim_max_iterations(10000)

#Solve
solution = study.solve()

#Create full stride
full = osim.createPeriodicTrajectory(solution)

#Print metabolics output
print(str(round(solution.getObjectiveTerm('met'),3))+' J kg-1 m-1')

#Visualise
study.visualize(full)

This seems to work pretty well. Solves in ~17 minutes on 8 threads, and produces a metabolics output of 0.253 J.kg-1.m-1.

Aaron

[Carlos Gonçalves]

Thanks, Aaron!!!

After that, I think I will also give it a go in the 0.5 early release!

Quick questions:

• metabolics = osim.Bhargava2004Metabolics() ---- this is a probe? Does it matter if you use model.addComponent or model.addProbe? I searched in 0.4 documentation and couldn't find it.

• "Solves in ~17 minutes on 8 threads" --- you had to enable some multithreading in your Python project?

Thanks a lot, again!

Best regards.

Carlos

[Karthick Ganesan]

Hi All,
I also tried the 0.5 release (in Matlab). I tried it with Ross Miller's 2D model. It gave a cost of transport of 2.3 J/kg/m. It is of the order of 1J/kg/m lesser than realistic. Umberger's probe gave a value which is almost twice of this. I read in the literature that Bhargava's model could give a low estimate of metabolic cost in tracking simulations. With compliant tendon dynamics it seems to take considerably more time.
I would like to confirm whether the output includes basal heat rate. I could not find a method for setting the basal heat rate on. However getBasalCoefficient method gave a value of 1.2. I also want to know whether this implementation includes the length dependence for the maintenance heat rate.
I am curious how Falisse could get a realistic cost with this model. May be the predictive nature of the simulation or differences in the muscloskeletal parameters or a 3D model?

Thanks,
Karthick.

[Aaron Fox]

Hi Carlos and Karthick,

In response to Carlos' questions:

• I believe it is a probe - so I suspect you could add it in a different way. I think the major difference with this would be how you pass the total metabolics calculation to the output goal. If in the probe set you would probably need to include this as part of the output path.

• I didn't enable anything specifically. In early versions of Moco parallel computing was disabled in Python - but I recall asking Chris about this recently and now it is meant to work normally

In response to Karthick's post:

• Now that you mention it the cost I reported is quite low - could be to do with the simplicity of the 2D model and simplified muscle set? It might be helpful to get Antonie's comments on this given he developed the C++ implementation for the 2D model, and see how realistic it is.

Aaron

[Aaron Fox]

Hang on - I may have realised something regarding the metabolics calculation in my provided code. I am extracting the output goal term to print the cost of transport at the end of the code, but the weight of this is goal is 0.1. Perhaps Chris or Nick could clarify that this would mean the output goal is actually 10% of the total metabolics calculation?

Aaron

[Nicholas Bianco]

Hey everyone,

The Bhargava2004Metabolics component is essentially a probe, but it is not the same thing as an OpenSim probe (i.e., it doesn't derive from the Probe base class). Therefore, you should add a Bhargava2004Metabolics component to your model via addComponent().

Aaron -- getObjectiveTerm() does include the cost weight, so you're correct that your output value is 10% of the total metabolics.

-Nick

[Nicholas Bianco]

Also, thanks for converting that example to Python, Aaron! Let's clean it up and put it into Moco.

[Karthick Ganesan]

I have converted that example into Matlab. May be that also could be added.