Experimental parameter identification of a multi-scale musculoskeletal model controlled by electrical stimulation: application to patients with spinal cord injury

Abstract

We investigated the parameter identification of a multi-scale physiological model of skeletal muscle, based on Huxley’s formulation. We focused particularly on the knee joint controlled by quadriceps muscles under electrical stimulation (ES) in subjects with a complete spinal cord injury. A noninvasive and in vivo identification protocol was thus applied through surface stimulation in nine subjects and through neural stimulation in one ES-implanted subject. The identification protocol included initial identification steps, which are adaptations of existing identification techniques to estimate most of the parameters of our model. Then we applied an original and safer identification protocol in dynamic conditions, which required resolution of a nonlinear programming (NLP) problem to identify the serial element stiffness of quadriceps. Each identification step and cross validation of the estimated model in dynamic condition were evaluated through a quadratic error criterion. The results highlighted good accuracy, the efficiency of the identification protocol and the ability of the estimated model to predict the subject-specific behavior of the musculoskeletal system. From the comparison of parameter values between subjects, we discussed and explored the inter-subject variability of parameters in order to select parameters that have to be identified in each patient.

Appendix: Multi-scale muscle model

The appendix focuses on the contractile element for which a multi-scale approach was used to state, from the microscopic scale using Huxley’s formulation to a macroscopic and controlled muscle model that fully fulfills the well-established muscle law. At the microscopic scale, the authors in [29] started from Huxley’s formulation, which describes the distribution of the fraction of actin–myosin pairs. The key point is the choice of rate functions of attachment f and detachment g of the cross-bridges that allow integration together with a meaningful f and g definitions. Let us begin with the inputs of the model. We designed a model with two separate inputs, i.e. the recruitment rate α and the activation u. The first one is linked to the number of fired motor units deduced from the level of stimulation intensity and pulse width. It is a static relationship described by a sigmoid function [27] and modulated here by the pulse width PW of the electrical stimulus:

$$ \alpha(\rm PW)=\frac{c_1}{1+\exp\{c_2(c_3-{\text{PW/PW}}_{\rm max})\}} $$

(10)

where c ₁, c ₂ and c ₃ are parameters which represent the plateau level, maximum slope and inflexion point, respectively. The u function is linked to the calcium dynamics and expresses the rate of the creation or breakage of bridges. It is composed of three phases which are the muscle contraction and relaxation, and the transition between both states. It is defined by the following function [29]:

$$ u(t) = \Uppi_{\rm c}(t) U_{\rm c} + (1-\Uppi_{\rm c}(t)) U_{\rm r} $$

(11)

where $\Uppi_{\rm c}$ is a function defined as follows:

$$ \Uppi_{\rm c}= \left\{ \begin{array}{lll} 1&\hbox{during the contraction phase }\tau_{\rm c}\\ \frac{\tau_{\rm r}-t_{\rm r}}{\tau_r} & \hbox{during the transition phase }\tau_{\rm r}\\ 0 & \hbox{otherwise} \end{array} \right. $$

(12)

where t _r is the relative time position from the beginning of the transition phase. U _c represents the rate of the actin–myosin cycle although U _r represents the rate when mainly breakage occurs during the calcium restoring phase in the sarcoplasmic reticulum.

For all subjects, the same values τ _c = 30 ms and τ _r = 20 ms were obtained through curve fitting of isometric muscle forces under stimulation twitches. The values of U _c and U _r were fixed for all subject as well to 30 s⁻¹ and 10 s⁻¹, respectively.

From the activation function above and when stimulus is under way, f (although f is always zero) can be defined as follows:

$$ f(\xi,t)= \left\{ \begin{array}{ll} \Uppi_{\rm c}(t) U_{\rm c} &\hbox{when } \xi \in [0,1], \\ 0 & \hbox{elsewhere} \end{array} \right. $$

(13)

where ξ is the relative position of the myosin head and the closest actin binding site. We consider g to be split in two parts, one linked to the chemical process from the cell at a rate u(t) that depends on the contraction phase and the other one linked to the relative shortening speed that may increase the detachment rate accordingly. It gives (compared to [29] the parameter a = 1):

$$ g(\xi,t) = u(t) + |\dot{\varepsilon}_{c}(t)| $$

(14)

Our definition is thus linked to the microscopic mechanisms involved within the muscle cell, but another assumption is now needed to change scale and to assess the macroscopic behavior. We do not compute the distribution n of the actin–myosin pairs fraction or express a specific distribution to solve the problem. Instead, we approximate the g function during the contraction phase as follows:

$$ g(\xi,t) = |\dot{\varepsilon}_{\rm c}(t)| $$

(15)

when f ≠ 0 and $\xi \in [0,1].$ It can be shown that a constant can be added to this formulation without any mathematical change if the definition of g has to be adjusted. This is in fact equivalent to changing the U _c and k ₀ values (see below). In any case it is consistent with Huxley’s theory, which expresses that g is lower within than outside of this interval. Then we can write the modified g function:

$$ g(\xi,t) = u(t) + |\dot{\varepsilon}_{\rm c}(t)| - f(\xi,t) $$

(16)

such that f + g no longer depends on the variable ξ. The dynamics of the first two moments of the distribution n is then possible. The detailed computation is available in [29], and the main steps are:

$$ \left\{ \begin{array}{ll} k(t) = k_0 \int_{-\infty}^{+\infty} n(y,t){\text{d}}y=k_0 M_0(t)\\ F(t) = k_0 h \int_{-\infty}^{+\infty} (y+y_0)n(y,t){\text{d}}y=k_0 h M_1(t) \end{array} \right. $$

(17)

with y linked to both the relative elongation of cross-bridges and relative position of thin and thick filaments, y ₀ is the initial position of the spring, k ₀ is the scaled stiffness at the sarcomere scale. The dynamics of these moments become:

$$ \left\{ \begin{array}{ll} \dot{M}_0(t) = \int_{-\infty}^{+\infty} f(y,t){\text{d}}y - \int_{-\infty}^{+\infty} (f+g)(y,t) n(y,t){\text{d}}y\\ \dot{M}_1(t) = \frac{S_0}{h} \dot{\varepsilon}_c(t) M_0(t) + \int_{-\infty}^{+\infty} (y+y_0) f(y,t){\text{d}}y - \int_{-\infty}^{+\infty} (y+y_0) (f+g) (y,t) n(y,t){\text{d}}y \end{array} \right. $$

(18)

Due to the wise definition of f and g, f + g does not depend on y, so integration is possible and we obtain:

$$ \left\{ \begin{array}{ll} \dot{M}_0 = \Uppi_{\rm c}(t) U_{\rm c} fl(\varepsilon_{\rm c}) - (u + |\dot{\varepsilon}_{\rm c}|) M_0\\ \dot{M}_1 = \frac{S_0}{h} \dot{\varepsilon}_{\rm c} M_0 + \frac{(1+2y_0) \Uppi_{\rm c}(t) U_{\rm c} FL(\varepsilon_{\rm c})}{2} - (u + |\dot{\varepsilon}_{\rm c}|) M_1 \end{array} \right. $$

(19)

fl is the force length relationship at the sarcomere scale, which represents the relation between the maximum number of actin–myosin pairs that can bind and the relative length of the contractile element. Then considering that all sarcomeres are identical, going to the fiber scale involves multiplying by a scale factor. And then going to the muscle scale involves scaling of the maximum available force and stiffness. This leads to Eq. (1) with macroscopic parameters which can be obtained by a combination of the microscopic parameters.