Automatic analysis of ultrasound shear-wave elastography in skeletal muscle without non-contractile tissue contamination

Abstract

Current analysis methods for obtaining mean shear modulus of skeletal muscles with ultrasound shear-wave elastography are limited by contamination with non-contractile tissues and manual operation of video processing. In this work, we develop new ultrasound image processing methods to assess muscle activities. We build upon previous research by using a 6-DOF robotic manipulator system and indirectly quantifying extensor carpi ulnaris (ECU) and triceps brachii longus (TRI) muscle elasticity during loading using ultrasound shear-wave shear modulus elastography. The purposes of this study were to (1) develop an automatic image-processing algorithm for removing non-contractile tissues from muscle elastography videos and (2) understand the effect of the removal on comparison of mean shear modulus of muscles across static motor tasks with variable muscle loadings in healthy humans. The developed algorithm with optimized clustering and thresholding identified and removed non-contractile tissues from muscle elastography videos with > 90% accuracy in arm muscles, causing reductions in the spatial variability of shear modulus data within the region of interest in healthy young adults. Removal of non-contractile tissues can alter the mean shear modulus of the muscles and influenced task comparisons by substantially altering the ranking of tasks by mean shear modulus.

Appendix

The computational muscle model used in this research is briefly presented below.

1.1 Inputs and outputs

The MATLAB (Mathworks, Natick, MA, USA) model predicts the force generation of 51 muscles across 12 joints of the human right arm and torso. The user specifies the baseline or nominal loading condition, target muscles, desired ratio γ_d of force change in the target muscles, a vector θ of joint angles, and the available degrees of freedom for the robot to induce γ_d. The model returns the realizable changes γ_a in muscle activation, and the experimental robotic loading F_d to induce γ_a. The chosen target muscles were triceps brachii longus (TRI L), and extensor carpi ulnaris (ECU). Muscles in the forearm, such as ECU, are smaller, closely packed, exhibit greater redundancy, and thus present a greater challenge for control and measurement of muscle activity (Beck et al. 2004; Madeleine et al. 2001). The value γ_d is the ratio of desired muscle force to nominal case muscle force (e.g. γ_d = 1.5 means desired muscle force is 1.5 times the nominal case muscle force). The loading conditions and γ_d for all experimental and nominal tasks performed during the study are listed in Table 1. Predicted nominal loading, as a percentage of maximum voluntary contraction, was 34.0% for TRI L and 22.4% for ECU. All γ_d values were realizable (i.e. γ_a = γ_d) based on the model.

1.2 Inverse solution for nominal endpoint loading

The computational prediction of human muscle force requires a mathematical description of how loads applied at the hand translate to proximal joint torques. For our study, a human subject grasps a robotic manipulator as the manipulator applies forces and torques to the hand. These loads are represented as a vector F. For a human musculoskeletal model of M muscles controlling movement about N joints, the torques at each joint are given by

$$ \varvec{\tau}_{h} = \varvec{J}\left(\varvec{\theta}\right)^{T} \varvec{F} + g\left(\varvec{\theta}\right) $$

(1)

where τ _h ϵ R^N is the human joint torque vector resulting from muscle forces. For our model, N = 12 and M = 51. Let θ be the joint angle vector (i.e., the posture of the arm), J be the Jacobian matrix relating the joint torques to end-point loading at the hand, and. g(θ) be the joint torque vector due to gravity. The dynamics of the body and robot were deemed negligible and omitted. The human joint torque vector is described by

$$ \tau_{h} = \varvec{A}\left( {\varvec{\uptheta}} \right)\varvec{f} $$

(2)

where A ϵ R^NxM is the moment-arm matrix of the muscles and f =  [f ₁ …f _M ] ϵ R^M is the human muscle force vector. Muscle force must be positive and has a physiological upper-bound such that 0 ≤ f _j ≤ j_MAXj, (j = 1,…,M). A has been obtained from other human musculoskeletal models and published biometric data (Ueda et al. 2010).

In the human body, the number of muscles exceeds the number of joints (Martini et al. 2013; Jacob et al. 1978), i.e., $M \gg N$. This muscle redundancy makes mapping of joint torques to muscle force an ill-problem. In order to describe the contribution of each muscle to each joint torque, an empirical optimization of Crowninshield’s cost function was used:

$$ u(\varvec{f}) = \sum\limits_{j = 1}^{M} {(c_{j} f_{j} )^{r} } \to \hbox{min} :{\text{subject to Eqs}}.\left( { 1- 2} \right) $$

(3)

where u (f) is a physiology-based cost function comprising the sum of muscle stress raised to a power. The summed muscle stress is minimized for a given task. c _j is a weighting factor of the j-th muscle, indirectly proportional to its physiological cross sectional area (PCSA). For our model, muscle stress is squared, i.e. r = 2, in accordance with (Crowninshield and Brand 1981). When τ_h is given, calculation of f to predict muscle force is straightforward.

1.3 Forward solution: Force changes and feasibility

In order to control muscle force, one must solve the problem of muscle force prediction given a desired f. The desired f is created by first selecting a nominal load F and predicting f. Vector f includes three categories of muscles: target muscles, non-target muscles (active during the nominal case but not be targeted for muscle control) and inactive muscles (inactive during the nominal case). This f is then modified such that the activity of target muscles are controlled (e.g. increased, decreased, or kept constant) based on γ_d.

Let f₀ be the nominal force vector when a subject is performing the nominal task. The desired forces are given by explicitly specifying the change ratios for each of the target muscles: f _td  = diag[γ₁…γ_Mt]^T f_t0. The subscript “d” denotes desired forces, and “0” denotes nominal forces. For a perfect control, the closed-form solution is as follows:

$$ {\varvec{\uptau}}_{ex} = \left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t}^{T} } & {{\mathbf{A}}_{n}^{T} } \\ \end{array} } \right]w^{ - 1} \left( {\left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t} } \\ {{\mathbf{A}}_{n} } \\ \end{array} } \right]\left( {{\mathbf{A}}_{t}^{ + } [w({\mathbf{f}}_{td} ) - w({\mathbf{f}}_{t0} )] + ({\mathbf{I}} - {\mathbf{A}}_{t}^{ + } {\mathbf{A}}_{t} ){\varvec{\upbeta}}} \right)} \right) $$

(4)

where w(*) is a function defined as $w(f_{i} ) = \partial u({\mathbf{f}})/\partial f = rc_{i} f_{i}^{r - 1}$. ${\varvec{\upbeta}}$ is a parameter vector that determines the influences on non-target muscles. The following feasibility conditions (FC) must hold:

• FC 1: The target muscles forces must be exactly realized

• FC 2: Inactive muscles stay inactive

• FC 3: Non-target muscles remain active, exert positive force

• FC 4: The external joint loading, τ _ex , is ideal, such that all joint torques τ _h are fully controllable by τ _ex , where τ _ex = J(θ)^TF.

According to FC1, the model must exactly realize the target muscle forces such that $\left| {f_{\text{td}}^{\text{T}} f_{\text{nd}}^{\text{T}} 0^{\text{T}} } \right|^{\text{T}} \, = \,{\text{arg min}}\;u(\varvec{f})$ with minimum changes in non-target muscles, i.e., $\left| {f_{\text{nd}} {-}f_{{{\text{n}}0}} } \right| \to { \hbox{min} } .$ The solution is not trivial since f_td must be physiologically feasible, as described in Eq. (3).

Our previous studies with exoskeletons showed that the physical interaction between the exoskeleton and the anatomy of the subject presented limitations in joint control such that the system was inherently non-ideal. This violation of FC4 prompted a switch toward a non-wearable endpoint loading system that trades greater but still imperfect control for greater adjustability and simplicity. The computational model was then modified for relaxed control and the use of a non-wearable robotic manipulator.

Under relaxed control, the initial formulation as shown in Eqs. 1 thru 3 remains the same. However, the aforementioned limitations prompted the use of a tolerance term and generalization of the solution to make the control scheme more widely applicable. This relaxed method therefore reframed the feasibility conditions as optimization criteria to minimize errors in the controlled muscle forces. For FC 1, it is assumed that requiring f = f_d, is necessary and possible. In reality, neither is true. Therefore, we introduce a tolerance term e such that f = f_d + e. FC 2 and 3 stipulate that inactive muscles must stay inactive, and non-target active muscles must stay active. For this reason, the revised model allows initially active non-target muscles to become inactive and initially inactive muscles to become minimally active. This is achieved by introducing new Lagrange multipliers in the KKT condition:

$$ \nabla u({\mathbf{f}}) + \sum\limits_{i = 1}^{M} {\mu_{i} } \nabla h_{i} (\varvec{f}) + \sum\limits_{j = 1}^{N} {(\lambda_{jl} } \nabla g_{jl} (\varvec{f}) + \lambda_{ju} \nabla g_{ju} (\varvec{f}) = 0 $$

(5)

where the following apply :$h_{i} \left( f \right) = \tau_{i} - a_{i}^{T} f,\;g_{ji} \left( f \right) = - f_{j} ,\;g_{ji} \left( f \right) = f_{j} - f_{\hbox{max} j} .$

We then defined ${\mathbf{q}} = \left[ {\begin{array}{*{20}c} {{\mathbf{q}}_{t} } & {{\mathbf{q}}_{n} } & {{\mathbf{q}}_{v} } \\ \end{array} } \right]^{T} = w(\left[ {\begin{array}{*{20}c} {{\mathbf{f}}_{t} } & {{\mathbf{f}}_{n} } & {f_{v} } \\ \end{array} } \right]^{T} )$. Equation (5) is then represented by:

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{q}}_{t} } \\ {{\mathbf{q}}_{n} } \\ {\mathbf{0}} \\ \end{array} } \right] = {\mathbf{A}}^{T} {\varvec{\upmu}} + {\varvec{\uplambda}} = \left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t} } \\ {{\mathbf{A}}_{n} } \\ {{\mathbf{A}}_{v} } \\ \end{array} } \right]{\varvec{\upmu}} + \left[ {\begin{array}{*{20}c} { - \lambda_{tu} } \\ {\lambda_{nl} - \lambda_{nu} } \\ {\lambda_{vl} } \\ \end{array} } \right]. $$

(6)

In Eq. (6), μ was a transformed torque vector used when the transformed force, q, was used instead of f. λ consisted of Lagrange multipliers. With the introduced error tolerance, where ${\varvec{\upvarepsilon}}$ was a transformed error, e, when q was used instead of f, the conditions could be given by:

Nominal case:

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{q}}_{t0} } \\ {{\mathbf{q}}_{n0} } \\ {\mathbf{0}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t} } \\ {{\mathbf{A}}_{n} } \\ {{\mathbf{A}}_{v} } \\ \end{array} } \right]{\varvec{\upmu}}_{0} + \left[ {\begin{array}{*{20}c} { - \lambda_{tu0} } \\ { - \lambda_{nu0} } \\ {\lambda_{vl0} } \\ \end{array} } \right] $$

(7)

Desired case:

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{q}}_{td} + {\varvec{\upvarepsilon}}_{t} } \\ {{\mathbf{q}}_{nd} } \\ {{\varvec{\upvarepsilon}}_{v} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t} } \\ {{\mathbf{A}}_{n} } \\ {{\mathbf{A}}_{v} } \\ \end{array} } \right]{\varvec{\upmu}}_{d} + \left[ {\begin{array}{*{20}c} {\mathbf{0}} \\ {\lambda_{nld} - \lambda_{nud} } \\ {\lambda_{vld} } \\ \end{array} } \right] $$

(8)

In perfect control, ${\varvec{\upvarepsilon}}_{t} = {\varvec{\upvarepsilon}}_{v} = \lambda_{nl} = \lambda_{nu} = \lambda_{tu} = 0$ so that the desired forces were exactly achieved. This relaxed representation allowed a broader range of muscle activation. Furthermore, the relaxed muscle control never failed and was always robust and feasible if an appropriate error tolerance was assigned. The relaxed version of Eq. 4 was

$$ {\varvec{\uptau}}_{ex} = \left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t}^{T} } & {{\mathbf{A}}_{n}^{T} } \\ \end{array} } \right]w^{ - 1} \left( {\left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t} } \\ {{\mathbf{A}}_{n} } \\ \end{array} } \right]\left( {{\mathbf{A}}_{t}^{ + } [w({\mathbf{f}}_{td} ) - w({\mathbf{f}}_{t0} ) + {\varvec{\upvarepsilon}}_{t} - {\varvec{\uplambda}}_{tu0} ] + ({\mathbf{I}} - {\mathbf{A}}_{t}^{ + } {\mathbf{A}}_{t} ){\varvec{\upbeta}}} \right)} \right) $$

(9)

To summarize, the relaxed conditions were as follows:

• RC 1: Target muscle forces were perfectly realized if:

  $$ {\text{rank}}[{\mathbf{A}}_{t} ] = {\text{rank}}[\begin{array}{*{20}c} {{\mathbf{A}}_{t} } & {w({\mathbf{f}}_{td} ) - w({\mathbf{f}}_{t0} ) + {\varvec{\upvarepsilon}}_{t} - {\varvec{\uplambda}}_{tu0} } \\ \end{array} ] $$

• RC 2: Inactive muscles were kept inactive if:

  $$ {\mathbf{A}}_{v} {\varvec{\upalpha}} + {\varvec{\uplambda}}_{vl0} > - \varepsilon_{v} ,\exists {\varvec{\upalpha}} $$

• RC 3: Non-target muscle forces remained positive if:

  $$ w({\mathbf{f}}_{\hbox{max} } ) > {\mathbf{A}}_{n} {\varvec{\upalpha}} + w({\mathbf{f}}_{n0} ) > 0,\exists {\varvec{\upalpha}} $$

In addition to implementing relaxed control, we also considered the limitations of the robotic system. The preceding treatment still assumes the robot is able to apply loads to each of the human’s joints and to the hand. With a robotic manipulator, the system can only apply torques to the hand and all other torques are set to 0. Therefore, an adequate robot loading F may not be available even if f_td is exactly realizable and a solution of τ_ex exists mathematically, the accompanying. This non-ideal robotic system may result in errors in target muscle forces. To overcome this issue, the realizability of τ_ex must be analyzed based on necessary rank conditions. The total torque may be rewritten as:

$$ {\varvec{\uptau}}_{ex} = {\mathbf{J}}^{T} {\mathbf{F}} - {\varvec{\uptau}}_{a} = \left[ {\begin{array}{*{20}c} {{\mathbf{J}}^{T} } & { - {\mathbf{E}}} \\ \end{array} } \right]\,\left[ {\begin{array}{*{20}c} {\mathbf{F}} \\ {{\varvec{\uptau}}_{ac} } \\ \end{array} } \right] = {\mathbf{J}}_{e}^{T} {\mathbf{F}}_{e} $$

(10)

where $\varvec{\tau}_{ac}$ is a vector that only contains the torque components controllable by the exoskeleton. E is an expansion matrix that maps these components onto the joint space. While previous studies have left r unspecified, this study explicitly used a cost function with r = 2 (Crowninshield and Brand 1981) to obtain a closed-form solution. The cost function is now quadratic: $u\left( \varvec{f} \right) = \varvec{f}^{T} \varvec{Sf}$ where $\varvec{S} = {\text{diag}}[c_{1} , \cdots c_{N} ]$, and $\varvec{q} = \varvec{w}\left( \varvec{f} \right) = 2\varvec{S}^{ - 1} \varvec{f}$. Substituting Eq. (6) for Eq. (4) yields the solution of $\varvec{ F}_{\varvec{e}}$, for relaxed control (Gallagher et al. 2013):

$$ {\mathbf{F}}_{e} = {\mathbf{G}}^{ + } (w({\mathbf{f}}_{td} ) - w({\mathbf{f}}_{t0} )) + ({\mathbf{I}} - {\mathbf{G}}^{ + } {\mathbf{G}})\psi $$

(11)

where ${\mathbf{G}} = \left( {2\left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t}^{T} } & {{\mathbf{A}}_{n}^{T} } \\ \end{array} } \right]\,{\mathbf{S}}\,\left[ {\begin{array}{*{20}c} {{\mathbf{A}}_{t} } \\ {{\mathbf{A}}_{n} } \\ \end{array} } \right]} \right)^{ - 1} {\mathbf{J}}_{e}^{T}$ and $\psi$ is a free vector parameter mapped onto the null space of G. This solution does not guarantee perfect realization, but minimizes errors in the target muscles in the sense of least mean square. If ${\text{rank}}[\varvec{G}] = {\text{rank}}[\varvec{G w}\left( {\varvec{f}_{td} } \right) - \varvec{w}(\varvec{f}_{t0} \varvec{ })]$ (rewritten FC1), perfect control is realized. Merging results given in Eqs. (9, 11) will directly determine a motor task from desired target muscles forces with the relaxed conditions. Investigating the association between the range of allowable error tolerance and induced muscle forces will provide quantitative assessment of the approach.