Adjust Rajagopal model tendon stiffness

[mrrezaie]

Hi, in the study by Luis et al. (2022), triceps surae and quadriceps tendons were modeled more compliant than others (stiffness: 15, 20, and 35, respectively; please see the provided link), and I was thinking that it would be great if you update the stiffness of these seven tendons in Rajagopal model accordingly.

I know that tendon stiffness is being handled by strainAtOneNormForce in Millard2012EquilibriumMuscle; so according to the tendon force-length formulation in De Groote et al. (2016) (S1), the stiffness constants can be converted to the percentage of tendon slack length (TSL) as follows:

DeGroote (kT)	Millard (% TSL)
35 [Default]	4.73 [Default=4.9]
20 [Quadriceps]	8.59
15 [Achilles]	11.59

This would help physiologically based musculoskeletal modeling. Thank you in advance.

[nickbianco]

Hi @mrrezaie, good idea! Would you be interested in opening a pull request for this?

We should make sure that convert between the Millard and DeGroote muscles produces the expected values from your table above.

[mrrezaie]

Hi, thanks for your response. Sure, I will open a pull request.

There are a couple of points about this conversion:

1. It's been said that the standard stiffness of 35 in the De Groote muscle has a shift of zero [REF]. But it leads to negative force multipliers at the beginning. I calculated the exact shift for this stiffness and it would be 1.17507567e-02. Now, the curve exactly starts from zero at zero strain.
2. The default strainAtOneNormForce in Millard muscle is 0.049 while the converted value from De Groote muscle with standard stiffness of 35 would be 0.04735946664256785. I think it's better not to change the default strainAtOneNormForce in Millard muscle as the values are relatively close. Moreover, as I said previously, the mathematical functions in De Groote paper do not seem to be so accurate.

Here is how I converted the values:

import numpy as np
import opensim as osim
import matplotlib.pyplot as plt

# normalized tendon length
NTL  = np.linspace(1., 1.5, 10000)

# default Millard tendon curve
MillardTendonCurve = osim.TendonForceLengthCurve()
# MillardTendonCurve.getStrainAtOneNormForce() # default value = 0.049 
plt.plot(NTL, [MillardTendonCurve.calcValue(i) for i in NTL], color='k', label='Millard default\n(0.049)')

# Tendon force-length multiplier De Groote et al. (2016) (S1)
c1=0.200; c2=0.995; c3=0.250

def calcShift(kT):
    # the force multiplier standard stiffness at zero strain is negative 
    # TFLM35 = c1*np.exp(35*(1.-c2))-c3+0.
    TFLM35 = c1*np.exp(35*(1.-c2))-c3+1.17507567e-02
    TFLMkT = c1*np.exp(kT*(1.-c2))-c3
    return TFLM35-TFLMkT

for kT in [35,20,15]:
    shift = calcShift(kT)
    TFLM = c1*np.exp(kT*(NTL-c2))-c3+shift

    # intersect with y=1 (Strain At One Norm Force)
    x = np.interp(1, TFLM, NTL)
    xdiff = x - 1

    MillardTendonCurve = osim.TendonForceLengthCurve()
    MillardTendonCurve.setStrainAtOneNormForce(xdiff)

    plt.plot(NTL, [MillardTendonCurve.calcValue(i) for i in NTL], label=f'Millard ({round(xdiff,3)})')
    plt.plot(NTL, TFLM, label=f'DeGroote ({kT})', linestyle='--')
    plt.scatter(x, 1, s=30, color='r')

    # percentage of TSL change
    print(f'kT: {kT} == strainAtOneNormForce: {xdiff}')

plt.axhline(1, color='k', linestyle='dotted', alpha=0.35)
plt.xlim(0.999,1.125)
plt.ylim(-0.01,1.2)
plt.title('Tendon force-length curves: Millard vs. De Groote')
plt.xlabel('Tendon Strain')
plt.ylabel('Normalized force')
plt.legend()
plt.show(block=False)

kT: 35 == strainAtOneNormForce: 0.047089602027437394
kT: 20 == strainAtOneNormForce: 0.0854568043378019
kT: 15 == strainAtOneNormForce: 0.11531044186005412

Please let me know your ideas. Thank you.

[nickbianco]

@mrrezaie, could you make these comparisons using the current DeGrooteFregly2016Muscle, which already has some modifications to the published curves? I would see if your modifications are necessary based on what's already in our implementation.

Also, this class uses the property tendon_strain_at_one_norm_force to set the stiffness, rather than setting kT directly. So you could just set strain values directly rather than converting from kT. You should pull the strain values directly from the papes cited by Luis et al. (2022), if they report them.

[mrrezaie]

I studied the OpenSim implementation of De Groote et al. (2016) muscle and saw the improvements. Thanks for the link.

The study by Luis et al. (2022) and the articles cited for triceps surae tendon stiffness, Swinnen et al. (2019) and Delabastita et al. (2020), all were conducted by Dr. De Groote's research team at KU Leuven and they only reported tendon stiffness, kT, rather than tendon strain. The value was determined according to agreement between simulated fiber length and in vivo muscle fascicle length changes measured by ultrasound.

In the references cited for quadriceps tendon stiffness, kT was also determined according to in vivo muscle fascicle/sarcomere length changes, not tendon/aponeurosis.

Please let me know what you think about it.

[nickbianco]

@mrrezaie, thanks for the update. Good to know that the kT values are based on ultrasound measurements. In that case, I think we could update the tendon_strain_at_one_norm_force for these muscles in the RajagopalLaiUhrlich2023.osim model. I think we should leave Rajagopal2016.osim with the default values.

I agree that we shouldn't change any values from 0.049 to 0.0474. But we should update values for muscles with kT = 20 or kT = 15. Based on the expression for kT in DeGrooteFregly2016Muscle, values of 0.089 and 0.119 would be reasonable for tendon_strain_at_one_norm_force.

If you make these changes in RajagopalLaiUhrlich2023.osim, could you also update the README file, including references?

[mrrezaie]

Sure, I'm going to create a pull request. Thank you.

[mrrezaie]

Hi, I just read something about tendon mechanics:

In the linear region of the stress–strain curve, where the tendon is stretched less than 4%, collagen fibers lose their crimp pattern. ... If the tendon is stretched over 4%, microscopic tearing of tendon fibers occurs. Beyond 8–10% strain, macroscopic failure occurs. And further stretch causes tendon rupture (Butler et al., 1978)

Wang J. H. (2006). Mechanobiology of tendon. Journal of biomechanics, 39(9), 1563–1582. https://doi.org/10.1016/j.jbiomech.2005.05.011

So, these changes might not be necessarily helpful.