Changes in the Neural and Non-neural Related Properties of the Spastic Wrist Flexors After Treatment With Botulinum Toxin A in Post-stroke Subjects: An Optimization Study

Introduction

Spasticity is a motor symptom commonly seen after a lesion of the upper motor neuron, e.g., after stroke, spinal cord injury, or cerebral palsy. Clinically, spasticity is defined as a velocity-dependent increase in tonic stretch reflex with exaggerated tendon jerks as a component of upper motor neuron lesions (Lance, 1980). In current practice, it is commonly measured subjectively by rotating a joint and estimating the resistance according to an ordinal scale, such as the Modified Ashworth Score (MAS) (Ashworth, 1964). MAS is usually employed to evaluate the clinical effectiveness of Botulinum toxin type A (BoNT-A), one of the most popular anti-spastic treatments aimed at reducing focal spasticity. The effects of BoNT-A include the decreased resistance to muscle stretch and increased the passive range of movement (Ozcakir and Sivrioglu, 2007). However, it is generally agreed that MAS has limited reliability and reproducibility (Fleuren et al., 2010), and is not able to discriminate between the underlying neural (stretch reflex) and non-neural contributions of resistant joint torque (Bakheit et al., 2001). Thus, alternative methods have been developed to objectively quantify resistant joint torque and its components (Harlaar et al., 2000; Mirbagheri et al., 2001; van der Helm et al., 2002). Several laboratories have developed custom-built apparatuses to examine joint stiffness in terms of underlying neural and mechanical properties, mostly at the ankle joint (Mirbagheri et al., 2001; Lorentzen et al., 2010; de Vlugt et al., 2013; Blanchette et al., 2015) and elbow joint (Ju et al., 2000; Calota et al., 2008). Recently, Lindberg et al. (2011) have developed a mechanical device that functions as an instrumented version of the MAS test to quantify spasticity of wrist and finger flexors and the test is convenient for daily clinical practice. The instrument (NeuroFlexor, Aggero MedTech AB, Solna Sweden) passively extends the wrist at two different velocities. A force transducer measures the resistance throughout the movements. The real-time analysis method, i.e., the NeuroFlexor Method (NF-Method), can separate the contributions to the measured resistance into a passive elastic component, passive viscous component, and a neural component. The validity of the NF-Method was demonstrated in chronic stroke patients through the strong correlation between the neural component and the electromyography (EMG) activity of flexor carpi radialis (FCR), both across all subjects and in each subject during the nerve block test (Lindberg et al., 2011). Recent investigations have also demonstrated that the NF-Method is reliable (Gäverth et al., 2013) and sufficiently sensitive to detect changes in spasticity by intramuscular injections of BoNT-A in wrist and finger muscles (Gäverth et al., 2014). However, the NF-Method does not allow for identifying effects on detailed mechanical properties of the joint nor the transformation of the motoneuron pool and muscle spindles.

Several theories have been proposed to explain the cause of spasticity, including decreased threshold or increased gain of the stretch reflex, hyper-excitability of the α-motor neuron and hypersensitivity of muscle spindles (Powers et al., 1988; Katz and Rymer, 1989; Sheean, 2002). The most common techniques of assessing neurophysiological outcomes after BoNT-A treatment include surface electromyography (EMG) (Hesse et al., 1996; Pandyan et al., 2002; Albani et al., 2010), electrical nerve stimulation (ENS) (Pauri et al., 2000; Sheykh et al., 2005; Manca et al., 2010) and stretch-reflex threshold (SRT) (Stampacchia et al., 2004; Chen et al., 2005; Lee et al., 2008). In some studies, EMG activity is recorded during the passive or active movement, and the magnitude or timing of the EMG activity is used to quantify the effect of the BoTN-A (Burgen et al., 1949). However, EMG magnitude can be highly variable between days, and the lack of standardized normalization can potentially lead to a misleading interpretation of treatment efficacy (Phadke et al., 2012). In the ENS, electrical stimulation of a peripheral nerve is used to evoke the Hoffmann's reflex (H-reflex), M-wave, or F-wave. The responses of the stimulations have been used as a potential outcome of the treatment (Sheykh et al., 2005). Although the ENS is widely used, these electrically evoked responses often do not correlate with the responses evoked with mechanical activation of the reflex, i.e., muscle stretching (Leis et al., 1996; Morita et al., 1998). The SRT is usually defined as the smallest angle or the lowest velocity that elicits muscle activity during passive stretching, which offers another alternative to measure the neurophysiological effect of BoNT-A. The SRT can provide meaningful information on the sensitivity of the neurophysiological outcomes after BoNT-A as long as a standard movement velocity and intensity can be consistently delivered (Phadke et al., 2012). Regardless of the varied limitations of the abovementioned techniques, direct measurement of motoneuron, and biophysical properties of human is not possible (Hidler and Rymer, 1999). In order to investigate the potential pathophysiological mechanisms of spasticity, a forward neuromusculoskeletal model that accounts for dynamics of muscle spindles, motoneuron pools, muscle activation and musculotendon of the wrist flexors has therefore been developed recently (Wang et al., 2017). The model can describe the overall resistant behavior of the wrist joint during passive extension test of spasticity. By tuning the passive and stretch-related parameters to fit the measured torque, neural and non-neural related properties of the spastic wrist flexors were estimated (Wang et al., 2017). Furthermore, the model explicitly accounts for motoneuron pool and muscle spindle behavior unlike the aforementioned NF-Method. The validity of the proposed optimization scheme was demonstrated through good fit to data, overall robustness to fluctuations in parameters and agreement of available EMG recording and the model estimation in the chronic stroke survivors and healthy controls (Wang et al., 2017). In line with a previous study by Koo et al. (Koo and Mak, 2006), two motoneuron pool parameters (threshold and gain) in the model have been shown to be physiologically meaningful and capable of characterizing the overall transmission of the motoneuron pool of spastic muscles. The aim of this work is to demonstrate the feasibility of the proposed modeling and optimization framework for investigating the potential neurophysiological effects of BoTN-A treatment. It is established that intramuscular administration of BoTN-A results in a decreased output from neuromuscular junction by inhibiting the release of acetylcholine from presynaptic motor neurons (Jankovic, 2004). However, how the BoTN-A injection influences the transmission of the overall motoneuron pool has not yet been investigated, even though it is believed that the injection causes decreased neural excitations that can paralyze the spastic muscles (Jankovic, 2004).

All in all, the primary objective of this study is to evaluate the effect of the BoNT-A treatment on the potential spasticity-related manifestations in persons with stroke using the proposed modeling framework. The data for fitting the model were previously collected using NeuroFlexor apparatus (Gäverth et al., 2014). Although they were already used to investigate the sensitivity of the NF-Method to the changes in spasticity after the BoTN-A injection (Gäverth et al., 2014), the approach adopted in the current study allows for separating motoneuron pool and muscle spindle effects as they are explicitly accounted for in the model, and thus provides additional insights in the mechanisms of BoTN-A efficiency. The key finding was that stroke survivors exhibited decreased neural components at 4 weeks post BoNT-A injection, which returned to baseline levels after 12 weeks. Based on the model, the decreased neural component was predicted to stem from the increased motoneuron pool threshold.

Method

Subjects

A sample of 21 persons with chronic stroke (ST) (18 males/3 females, mean ± SD, age: 50 ± 12 yrs, body weight: 76 ± 13 kg) was selected from previous studies. All participants were recruited from a local rehabilitation clinic (Danderyd University Hospital, Department of Rehabilitation Medicine, Stockholm, Sweden) and received treatment of BoNT-A (Botox^®, Allergan, Irvine, CA, USA) in one or more of their wrist and/or finger flexors after assessing by physicians (Table 1). Inclusion criteria were the following: (1) stroke > 6 months prior to inclusion, (2) spasticity in wrist and finger flexors, with a MAS ≥ 1, (2) at least 50° passive range of motion (ROM) at the wrist joint (20°Flexion to 30° extension), and (3) ability to understand and comply with instruction. The physician selected the muscles and the dose based on clinical experience and goal-setting. No additional treatment or hand training was prescribed. The participants were examined before and after intramuscular injections of BoNT-A in the spastic wrist and finger flexor muscles. The NeuroFlexor data was collected on 3 occasions: at baseline (maximum 2 weeks before treatment), 4 weeks after treatment, and 12 weeks after treatment. Four of the participants dropped out from the assessment at 12 weeks, due to illness, problems with transportation, or time constraints. The passive ROM of the wrist was also measured using a goniometer at each occasion. All participants gave written informed consent according to the Declaration of Helsinki. The study was approved by the Regional Ethics Committee, Karolinska Institutet, Stockholm, Sweden. The neuromusculoskeletal model and optimization framework was implemented using Matlab 2014b (The MathWorks, Inc., Natick, US).

TABLE 1

Table 1. Clinical description of participants (FDS, flexor digitorum; FCR, flexor carpi radialis; FDP, flexor digitorum profundus; FCU, flexor carpi ulnaris).

Experimental Setup and Protocol

The resistance torque was measured using the NeuroFlexor that produced constant movements at a slow (5°s⁻¹) and a fast (236°s⁻¹) velocity (see details in Lindberg et al., 2011). The participants were seated comfortably with elbow in 90° of flexion and the hand securely fastened into the device with the wrist axis of rotation aligned with device (Figure 1A). In order to eliminate measurement variance, a clam and relax environment was carefully maintained during testing. The hand rested on the platform with the metacarpophalangeal joints in slight flexion and the fingers fully extended. Both the hand and arm were fastened using non-elastic velcro straps in order to ensure that movement could only occur at the wrist joint. The instrument performed slow and fast movements in a 50° ROM (20°Flexion to 30° extension) with a starting angle of 20° of palmar flexion. For each participant and test occasion, 5 slow movements and 10 fast movements were recorded. Resistance profiles were also obtained when the device ran empty (without hand; see an example data in Figures 1B,C) to enable the biomechanical model to isolate forces originating from the hand. As input to the neuromusculoskeletal model and optimization framework, the measured data of raw torque and angle were low pass filtered with a cut-off frequency of 10 Hz (second-order Butterworth filter).

FIGURE 1

Figure 1. (A) The NeuroFlexor instrument (Photo printed with permission from Aggero MedTech AB, Stockholm Sweden). (B,C) Example measurements of resistant force profiles during slow (B) and fast movements (C) in one representative stroke patient. Red trace represents the angle of wrist movement (from flexion to extension). Black trace represents one resisting force trial and black dashed trace represents one resistant force trial when device runs without hand. Three time points are identified by the NeuroFlexor Method (NF-Method): P3 1 s after slow passive stretch; P1 the first peak and P2 the peak toward the end of the fast movement.

The NF-Method

The NF-Method has been previously presented (Lindberg et al., 2011; Pennati et al., 2016), and it can separate the active muscle force component induced by the stretch reflex from the passive mechanical components. The method is briefly described below.

Three time points at the measured resistant force-time curve were defined specifically to estimate the inertial (IC), viscosity (VC), elasticity (EC), and neural components (NC): one point at the initial peak (P1) and one point at the later peak (P2) during the fast movement, and one point at fully stretched position (P3) during the slow movement.

Inertia component (IC) was defined as the force resisting the acceleration of the hand as:

$\begin{array}{rc}\text{IC}=\text{ma}& (1)\end{array}$

with m the mass of the hand and platform and a the acceleration.

Elasticity component (EC): The length –dependent elasticity was recorded at P3, 1s after the end of the slow movement, which included both the linear elasticity and nonlinear end range stiffness.

Viscosity component (VC): The early viscosity component was defined as the resisting force that remained after the inertia component had been subtracted from the initial peak of the total resistant force at P1.

$\begin{array}{rc}\text{V}{\text{C}}_{\text{p}1}=\text{Total forc}{\text{e}}_{\text{p}1}-\text{IC}& (2)\end{array}$

The late viscosity was estimated as:

$\begin{array}{rc}\text{VC}=(\text{Total forc}{\text{e}}_{\text{p}1}-\text{IC})\times 0.2& (3)\end{array}$

The neural component (NC) was estimated at P2 by subtracting the elasticity and viscosity component from the total force as:

$\begin{array}{rc}\text{NC}=\text{Total forc}{\text{e}}_{\text{p}2}-(\text{EC}+\text{VC})& (4)\end{array}$

Neuromusculoskeletal Model

The full description of the neuromusculoskeletal model can be found in our previous study (Wang et al., 2015, 2017). The model is only briefly described here.

Musculotendon Modeling

The musculotendon model was modified based on an lumped Hill - model described previously by Winters (Winters and Stark, 1985), which consisted of one lumped muscle actuator that represents synergistic wrist flexors and one passive unit. Wrist extensors were not included in the model due to their negligible activation during the passive wrist extension movement.

Equation of Motion

The modeled wrist-resistant torque was determined according to dynamic equilibrium during the passive wrist extension movement:

$\begin{array}{rc}{T}_m=(J_h+J_{hp})\ddot{\theta }+T_{MA}+T_{PP}+T_{Gh}+T_{Ghp}& (5)\end{array}$

with T_m the modeled resistant moment, J_h and J_hp the moments of inertia of the hand segment and hand-plate, $\ddot{\theta }$ the angular acceleration, T_MA the torque generated by the lumped muscle actuator, T_PP the torque generated by the passive unit and T_Gh and T_Ghp the torques due to the gravity of the hand segment and hand-plate, respectively.

Passive Unit

The passive unit in the model was considered as the passive viscoelastic elements of the muscle and connective tissues of the joint (Winters and Stark, 1985). Thereby, the elasticity component (EC) was an angle-dependent torque component of T_PP, and the viscosity component (VC) was an angular velocity-dependent torque component.

$\begin{array}{rc}EC=K_p(\theta -{\theta }_0)+k_1{e}^{k_2(\theta -{\theta }_0)-1}& (6)\end{array}$

and

$\begin{array}{rc}VC=B_p\stackrel{·}{\theta }& (7)\end{array}$

with K_p denoting the stiffness coefficient, B_p the viscosity coefficient, k₁ and k₂ the non-linear exponential coefficients, and θ₀ the resting joint angle at the muscle-tendon slack length, such that T_PP is equal to zero.

Lumped Muscle Actuator

The lumped muscle actuator was modeled as a contractile element in series with a series elastic element. The joint moment-angle and moment-angular velocity relationships were implemented in the actuator model (Winters and Stark, 1985).

$\begin{array}{rc}{T}_{MA}=a{f}_v(\stackrel{·}{\theta })f_l(\theta )& (8)\end{array}$

with a denoting the muscle activation, f_v the moment-angular velocity relationship, and f_l the moment-angle relationship. In passive wrist extension, muscle activation a was considered entirely induced by the stretch reflex of the wrist flexors. Therefore, the neural component (NC) was equivalent to T_MA.

Muscle Spindle Modeling

The modified hybrid v^0.6 model (Prochazka and Gorassini, 1998) was implemented to describe the firing characteristics of the muscle spindle in lengthening as a function of joint angle and angular velocity:

$\begin{array}{rc}r(t)=G_v{(\stackrel{·}{\theta }(t-{\tau }_{SRD}))}^{0.6}+G_l(\theta (t-{\tau }_{SRD})-{\theta }_0)+r_0& (9)\end{array}$

with r denoting the Ia afferent (the primary afferent) firing rate, θ₀ the rest joint angle, τ_SRD the stretch reflex delay, G_v and G_l the dynamic and static gains of the spindle model respectively, and r₀ the background discharge rate of the Ia afferent neuron. τ_SRD and r₀ were set to 30 ms and 6.45 imps⁻¹ (Wilson et al., 1999).

Motoneuron Pool Modeling

The input-output relationship for the α-motoneuron pool was represented by a Gaussian cumulative distribution function (Fuglevand et al., 1993). The motoneuron pool transformation, i.e., the motoneuron neural excitation as a function of the muscle spindle firing rate, was modeled as a sigmoid function which ranges from 0 to 1,

$\begin{array}{rc}u(t)=\frac{1}{1+e^{\frac{\mu -r(t)}{\sigma }}}& (10)\end{array}$

with u denoting the neural excitation of the motoneuron pool, μ the muscle spindle firing rate at 50% motoneuron recruitment, and σ the standard deviation of cumulative distribution function.

The threshold and gain of the motoneuron pool of wrist flexors were further calculated based on the definition from a previous study (Koo and Mak, 2006). A threshold μ₀ corresponds to the spindle-firing rate that results in neural excitation of 0.005, and gain G₀ is defined in the following manner:

$\begin{array}{rc}{G}_0=\frac{0.6}{r_{0.8}-r_{0.2}}& (11)\end{array}$

with r_0.8 and r_0.2 denoting the muscle spindle firing rate at 80 and 20% of full motoneuron recruitment.

Muscle Activation Dynamics

Muscle activation dynamics was described as the process of converting neural excitation u(t) to muscle activation a(t) as a first order differential equation (Thelen, 2003).

Optimization

A previous study indicated that the EMG activity of FCR was rarely observed in chronic stroke patients during the slow movement (Lindberg et al., 2011). Therefore, it was assumed that the resistant torque was entirely caused by the passive unit. Nonlinear least-squares (NLS) optimization was used to approximate the passive parameters ([K_p, B_p, k₁, k₂]) based on the experimental measurements of the resistant torque in the slow movement. A genetic algorithm (GA) (Golberg, 1989) was used to identify the optimal stretch reflex-related parameters ([G_v, G_l, μ, σ]) by minimizing the root-mean-square error (RMSE) between the measured and modeled resistant torque in part of the fast movement. Constraints were imposed on G_v and G_l to vary ±10% from the reported values (4.3 and 2.0 respectively).

Data Analysis and Statistics

A set of trials (one slow movement, one fast movement, one slow movement without the hand and one fast movement without the hand) from each participant's NeuroFlexor measurements were randomly selected to calculate the EC, VC, and NC at each assessment occasion using both NF-Method and optimization, respectively. The same trials were also used for computing passive and stretch reflex related parameters using the neuromusculoskeletal model and optimization. The validity of the model was determined for each trial by the Variance Accounted For (VAF) (de Vlugt et al., 2010):

$\begin{array}{rc}VAF=\left(1-\frac{\sum {\epsilon }^2(t)}{\sum {T_{meas}(t)}^2}\right)·100\%& (12)\end{array}$

with ε(t) denoting the error, i.e., the difference between the measured resistant torque, T_meas(t) and the modeled resistant torque, T_m(t), over t = 1, …, N samples. An extensive validity and reliability analysis of the used method and the estimated model parameters was performed previously (Gäverth et al., 2013; Wang et al., 2017).

Statistical analyses were performed using IBM SPSS Statistics (IBM Corp., Armonk, NY, USA). Descriptive statistics were expressed as median and range (or mean and standard deviation as appropriate). Spearman's rank order correlation test was used to investigate the correlation between components identified by the NF-Method and optimization, respectively. Differences over three time points in passive and stretch reflexed-related parameters were calculated using a linear mixed model with one subject effect: time, consisting of three time points (compound symmetry). A post-hoc pairwise comparison with Bonferroni adjustment was applied if the omnibus test was significant. Differences were considered statistically significant if p-value ≤ 0.05.

Results

Passive Parameters

As an example, Figures 2A,C shows the modeled and measured total resistant torque of one participant during the slow movement. It illustrates a good fit of the model. Consistent findings were observed in all participants. The VAF values of three assessment occasions were similar for optimization (pre BoNT-A: 95.2% ± 2.9%, 4 weeks post BoNT-A: 95.1% ± 2.9%, 12 weeks post BoNT-A: 93.1% ± 3.5%), which validated the generalization properties of the model fitted to individual participants. Compared to the baseline, there was a significant increase of non-linear coefficient k₁ (p = 0.02) at 12 weeks post BoNT-A. In addition, a significant decrease of passive ROM (p = 0.01) was determined at 12 weeks post BoNT-A. No significant changes in other passive parameters were found between test occasions (Table 2).

FIGURE 2

Figure 2. Measurements obtained with the NeuroFlexor and model based estimates of resistant torque for one representative stroke patient at three test occasions. The blue solid line represents the measured joint angle. The black solid line represents the measured resistant torque. The line with red circles represents the modeled resistant torque for the same subject. (A–C) In the slow movement, the wrist joint was passively extended from 20° (0.35 rad) to 30° (0.52 rad) of extension at 5°/s using a NeuroFlexor. A nonlinear-least squares optimization was used to predict passive parameters by fitting to the torque measured in the slow movement. (D–F) In the fast movement, the wrist joint was passively extended from 20° (0.35 rad) of flexion to 30° (0.52 rad) of extension at 236°/s (4.08 rad/s) using a NeuroFlexor. A genetic algorithm (GA) was used to estimate stretch reflex-related parameters. The optimization was only used to match data recorded in the period between PII to PIII.

TABLE 2

Table 2. Individual components, the passive and reflex related parameters (median and range) determined using modeling and optimization at baseline, 4 weeks and 12 weeks post BoNT-A injection.

Stretch Reflex-Related Parameters

Due to measurement error during the period 0.05–0.09s of the recording and the latency of the reflex, the optimization was only performed for a selected part of the fast movement (from the second peak, PII, to the third peak, PIII, of the measured total resistant torque, Figure 2E), during which the angular velocity remained stable. As an example, Figures 2D,F shows the modeled and measured resistant torque of one participant during the fast movement. It also illustrates a good fit of the model to the data. Consistent findings were observed in all participants. The VAF values of the three assessment occasions were similar for optimization (pre BoNT-A: 90.1% ± 2.1%, 4 weeks post BoNT-A: 89.8% ± 2.2%, 12 weeks post BoNT-A: 89.6 ± 2.5%).

For two motoneuron pool parameters (μ and σ), a significant increase in μ between baseline and 4 weeks post BoTN-A (p = 0.01) was observed, and a reduction of μ at 12 weeks after treatment was found, but the changes were not significant (Table 2). No significant changes were found in σ between test occasions. The participants had a significantly higher motoneuron pool threshold μ₀ at 4 weeks post BoTN-A, but no differences were found in the motoneuron pool gain G₀. Significant and high correlations between μ₀ and G₀ were determined in three test occasions (pre BoTN-A: p<0.01, r_s = 0.73, 4 weeks post BoTN-A: p<0.01, r_s = 0.75, 12 weeks post BoTN-A: p<0.01, r_s = 0.81, Figure 3). The averaged motoneuron pool profiles extracted before and after treatment are illustrated in Figure 4.

FIGURE 3

Figure 3. Highly correlated motoneuron pool parameters, threshold, μ₀, and gain, G₀, were determined at three test occasions: before BoTN-A treatment (A), 4 (B), and 12 (C) weeks after BoNT-A treatment. The solid lines represent the linear regression fit to the relationship of threshold and gain (r_s denotes the correlation coefficient).

FIGURE 4

Figure 4. (A,B) Model based estimates of muscle spindle firing rate induced by wrist extension in the slow (A) and fast (B) movement. (C,D) Mean motoneuron pool profiles were estimated using optimization at three test occasions (before BoTN-A treatment, 4 and 12 weeks after BoNT-A treatment). The dash lines in green, red and blue represent the motoneuron pool profile in a broader range of the spindle firing rate (0 imps⁻¹ to 100 imps⁻¹ for the slow movement; 0 imps⁻¹ to 250 imps⁻¹ for the fast movement) and the solid lines correspond to the range of the neuron excitation estimated in the current study.

Correlation of Individual Components From Model Estimation and the NF-Method

The EC, VC, and NC were identified at the third peak of the resistant torque in the fast movement, as the same as the NF-Method. To be consistent in the unit, individual torque component from model estimation was converted to the force component with a constant moment arm (the measured distance from the force sensor to the joint axis of NeuroFlexor). There was a significant reduction in NC between baseline and 4 weeks post BoNT-A (p = 0.01) and a significant increase between 4 weeks and 12 weeks post BoNT-A (p = 0.03). However, no significant difference in NC was found between baseline and 12 weeks after the treatment. There were no significant changes in EC and VC between test occasions. The individual components calculated from the current model estimation correlated well with those based on the NF-Method in each test occasion. A significant and high correlation was determined in NC and EC, as well as moderate correlation in VC (Table 3).

TABLE 3

Table 3. The correlation between components identified by the NF-Method and optimization method was determined using Spearman's rank order test at three assessment occasions.

Discussion

In this study, a previously developed neuromusculoskeletal model was applied to simulate the passive wrist extension test of spasticity before and 4 weeks and 12 weeks post BoNT-A injections in the chronic stroke patients. By modeling musculotendon, muscle spindle, and motoneuron pool parameter explicitly, the neural and non-neural related properties at the wrist joint were estimated using optimization techniques at each test occasion. Compared to the baseline, stroke survivors exhibited decreased neural component at 4 weeks post BoNT-A injection, which returned to the baseline levels after 12 weeks. The decreased neural component was predicted by the model to stem from the increased motoneuron pool threshold. In addition, the motoneuron pool threshold and gain were found to be highly correlated in all test occasions. To the best of our knowledge, this is the first study that has attempted to investigate the underlying neural and non-neural related changes in wrist muscles up to 12 weeks after BoNT-A injection by combining motorized measurements and neuromusculoskeletal modeling.

Individual Components and Comparison to the NF-Method

In order to examine whether there were systematic differences in component estimation using the NF-Method and the neuromusculoskeletal model, EC, VC, and NC were calculated at the same instance as the NF-Method (at the maximal resistance torque in the fast movement, PIII). At baseline and 12 weeks after injection, NC was the largest contributor to the total resistant torque followed by EC. At 4 weeks post BoNT-A, EC instead of NC became the largest because of the significant decrease in NC. This finding was in line with the action of botulinum toxin, where the effect was expected to be maximal at 4 weeks and to diminish after 12 weeks post the injection (Gäverth et al., 2014). No significant changes in EC and VC were found over the time as expected, since BoNT-A was not assumed to influence the mechanical properties of the muscles. Moderate to high correlations with the NF-Method were found in all three components over time (Table 3). Similar correlations were reported from our previous study, though only baseline data was analyzed (Wang et al., 2017).

Passive and Stretch Reflex Related Parameters

Passive Parameters

The passive wrist stiffness and viscosity at baseline in the current study were comparable to our previous report (Wang et al., 2017). No comparison data were found for stroke survivors after BoTN-A treatment in the literature. Compared to the baseline, although the linear stiffness (K_p) and viscosity coefficients (B_p) were similar, a tendency of the increased nonlinear exponential (k₁) term was observed over time and significance was determined at 12 weeks (Table 2). During the slow movement, the passive tissue demonstrated a quasi-linear resistant torque in response to the increased wrist angle. A higher non-linear coefficient (the exponential term) results in an earlier exponential increase in the resistant torque. It might indicate a reduced ROM of the wrist joint over the time, because larger resistance torque develops in the end of the range of motion, close to the participant' s maximal extension limits. Compared to the baseline, we found the passive ROM of the wrist was significantly decreased at 12 weeks, which verified our hypothesis. Shortening of the muscle and/or connective tissues often induced by longstanding spasticity is regarded as the major contributor to the passive movement limitation (de Bruin et al., 2014), which may cause an earlier start of the exponential rise of the stiffness torque. Studies on musculoskeletal stiffness after BoTN-A injection were often conducted with very different protocols, thus resulting in somewhat contradictory conclusions. Moreover, the effect of BoTN-A on passive muscle stiffness measured biomechanically in humans has rarely been reported in the literature. To our knowledge, the only study was conducted on children with cerebral palsy after an injection of BoTN-A in triceps surae (Alhusaini et al., 2011). In their study, no significant change in either passive range of movement or passive stiffness of the ankle muscles was reported 6 weeks post the injection. However, in animal models, the passive elasticity module of muscle fascicles increases 1 month post injection (Thacker et al., 2012). A recent study on rats also found that the passive stiffness of the injected limb was increased at 2, 4, and 8 weeks (Pingel et al., 2016).

Stretch Reflex Related Parameters

In the current study, the modeling of the reflex loop was simplified and consisted of muscle spindle and motoneuron pool models. The static and dynamic gain (G_l, G_v) of the hybrid V^0.6 model reflected the level of γ-motoneuron activation. Recent research has shown that the firing rate of muscle spindles in the stroke survivors was similar to those of the controls (Dietz, 2003). Therefore, we assumed that the hyperactive stretch reflex was mainly due to the abnormal input-output relationship of the motoneuron pool. The changes of NC were further investigated by identifying the motoneuron pool profile for each stroke survivor before and after treatment (Figure 4). By definition, all of the motor units that innervate a certain muscle, i.e., the lumped wrist flexors, are considered as a motoneuron pool. Each motor unit consists of a single motor neuron and the muscle fibers it innervates. Therefore, the overall motoneuron pool behavior reflects qualitatively the suppress effect of the BoTN-A treatment at the neuromuscular junction where synapse with the motor end plates of the muscle fibers. Significant differences were observed between the baseline and 4 weeks post BoNT-A injection (Table 2). According to our modeling and analysis, reduced hyperactive reflex was mostly due to higher μ. In order to achieve a more meaningful interpretation of each participant's motoneuron pool profile, the threshold μ₀ (the minimal spindle firing rate) and gain G₀ (the mid-range slope) were defined. Physiologically, G₀ can be interpreted as an intrinsic motoneuron property, e.g., passive membrane properties and the voltage sensitive membrane conductance (Koo and Mak, 2006), and μ₀ as a net excitatory and inhibitory inputs to the motoneuron pool. Compared to the baseline, μ₀ had a significant increase 4 weeks after the injection and returned back to the baseline level at 12 weeks, while no significant changes in G₀ were noticed over time. Similar findings using an experimental approach were reported by Stampacchia et al. (2004). In their study, they quantified the stretch reflex threshold as the sudden increase in the recorded EMG activity at the lowest wrist extension velocity. They found that the stretch reflex threshold increased after the botulinum toxin treatment. They explained that the observed increased stretch reflex threshold was mainly due to the intrafusal action of the drug, e.g., the reduced afferent discharge and stretch sensitivity in the muscle spindles. Theoretically, neural excitation can be enhanced by reducing μ₀ without changing the shape of the profile (G₀), i.e., shifting the profile curve to the left. However, we found high correlations between μ₀ and G₀ in all test occasions (Figure 3), in line with our previous study (Wang et al., 2017). This finding agreed with our previous assumption that the abnormal motoneuron pool excitability was possibly modulated through a combined effect from the overall synapse input and the intrinsic properties of motoneurons. However, under the intervention of BoTN-A, the dominant effect was from the former factor, i.e., the overall synapse input represented as μ₀ in our modeling approach. A significant increase in μ₀ at 4 weeks implied that less motor neuron could be excited after the BoTN-A injection. This is also in line with the pharmacological mechanisms of BoTN-A, which inhibits the neural transmission at the neuromuscular junction.

Clinical Implications

In combination with the NeuroFlexor measurement, a computational model accounting for musculotendon, muscle spindle, and motoneruon pool effects can help identify potential neural and non-neural related parameter changes along the stretch reflex pathway in the spastic patients after treatment of BoTN-A. These parameters are physiologically meaningful, but not directly measurable. Moreover, this approach can also be applied in other non-pharmacological interventions, e.g., physical interventions and orthotics. Monitoring neural and non-neural related parameters prior to and post treatment can provide insights into different intervention strategies. In addition, our findings can be used to improve the measurement protocol of the NF-Method. The current measurement protocol operates in a fixed ROM, which excluded the severely affected patients with limited ROM. Consequently, the estimation of EC was not sufficiently accurate to describe passive stiffness in each participant's full ROM. An individualized ROM would be more informative and appropriate for this heterogeneous population.

Limitations

Several simplifications were made in the model. First, the optimization was only performed in part of the fast movement. There is a short “mechanical chaotic state” after the first acceleration peak in the fast movement, which occurs around 0.06–0.07 s. During that period, the force sensor is unloading due to a small gap between the handplate and the force sensor. Thereby, no valid force data was collected. The improvement of the mechanical design of the NeuroFlexor is valuable in the future. Second, the goal of the current study was to use a general but valid forward neuromusculoskeletal model and optimization framework that allows for reliable estimation of passive muscle-tendon and stretch reflex-related parameters before and after the BoTN-A treatment in each person. However, the authors are aware that neurophysiological measurement, e.g., EMG activity, can provide more information of spasticity when combined with a resistance measurement. In addition, NeuroFlexor is not able to separate finger flexors and wrist flexors and our analysis cannot identify parameters on an individual muscle basis. Finally, the sample size of the patients is small and there was several drop-outs at 12 weeks. Future studies with larger cohorts will improve the generalization of our findings to the entire stroke population.

Conclusions

Neural and non-neural related properties of wrist flexors were estimated in persons with chronic stroke before and after BoTN-A injection using a forward neuromusculoskeletal model and optimization. The model describes the overall resistant behavior of the wrist joint during a motorized extension test of spasticity. Compared to the baseline, the neural contribution to the decreased resistant torque at 4 weeks was predominated by the increased stretch reflex threshold in the model. Moreover, an increased non-linear stiffness coefficient indicated a reduced wrist ROM over time, even though the treatment effect of reducing spasticity was observed. We concluded that in combination with the NeuroFlexor measurement, the proposed neuromusculoskeletal model and optimization scheme served as suitable tools for investigating potential parameter changes along the stretch reflex pathway in persons with chronic stroke. These neural related parameters are physiologically meaningful but not directly measurable. By monitoring their values before and after interventions, (e.g., BoNT-A treatment and other therapies), valuable insights can be provided for effective treatment response.

Author Contributions

RW conceived the model, developed it, analyzed the data. PH developed optimization framework. JG collected the data. All authors have been involved in manuscript preparation.

Conflict of Interest Statement

Author JG owns part of the commercial rights of the measurement instrument described in this study as shareholders of the manufacturing company Aggero MedTech AB.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by Promobilia Foundation, Norrbacka-Eugeniastiftelsen, and Carl Tryggers Stiftelsen.

References

Albani, G., Cimolin, V., Galli, M., Vimercati, S, Bar, D., Campanelli, L., et al. (2010). Use of surface EMG for evaluation of upper limb spasticity during botulinum toxin therapy in stroke patients. Funct. Neurol. 25, 103–107.

PubMed Abstract | Google Scholar

Alhusaini, A. A., Crosbie, J., Shepherd, R. B., Dean, C. M., and Scheinberg, A. (2011). No change in calf muscle passive stiffness after botulinum toxin injection in children with cerebral palsy. Dev. Med. Child Neurol. 53, 553–558. doi: 10.1111/j.1469-8749.2011.03930.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashworth, B. (1964). Preliminary trial of carisoprodol in multiple sclerosis. Practitioner 192, 540–542.

PubMed Abstract | Google Scholar

Bakheit, A., Pittock, S., Moore, A., Wurker, M., Otto, S., Erbguth, F., et al. (2001). A randomized, double-blind, placebo-controlled study of the efficacy and safety of botulinum toxin type A in upper limb spasticity in patients with stroke. Eur. J. Neurol. 8, 559–565. doi: 10.1046/j.1468-1331.2001.00277.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanchette, A. K., Mullick, A. A., Moïn-Darbari, K., and Levin, M. F. (2015). Tonic stretch reflex threshold as a measure of ankle plantar-flexor spasticity after stroke. Phys. Ther. 96, 687–695. doi: 10.2522/ptj.20140243

PubMed Abstract | CrossRef Full Text | Google Scholar

Burgen, A., Dickens, F., and Zatman, L. (1949). The action of botulinum toxin on the neuro-muscular junction. J. Physiol. 109, 10–24. doi: 10.1113/jphysiol.1949.sp004364

PubMed Abstract | CrossRef Full Text | Google Scholar

Calota, A., Feldman, A., and Levin, M. (2008). Spasticity measurement based on tonic stretch reflex threshold in stroke using a portable device. Clin. Neurophysiol. 119, 2329–2937. doi: 10.1016/j.clinph.2008.07.215

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J.-J. J., Wu, Y.-N., Huang, S.-C., Lee, H.-M., and Wang, Y.-L. (2005). The use of a portable muscle tone measurement device to measure the effects of botulinum toxin type a on elbow flexor spasticity. Arch. Phys. Med. Rehabil. 86, 1655–1660. doi: 10.1016/j.apmr.2005.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

de Bruin, M., Smeulders, M. J., Kreulen, M., Huijing, P. A., and Jaspers, R. T. (2014). Intramuscular connective tissue differences in spastic and control muscle: a mechanical and histological study. PLoS ONE 9:e101038. doi: 10.1371/journal.pone.0101038

PubMed Abstract | CrossRef Full Text | Google Scholar

de Vlugt, E., de Groot, J. H., Schenkeveld, K. E., Arendzen, J. H., van der Helm, F. C., and Meskers, C. G. (2010). The relation between neuromechanical parameters and Ashworth score in stroke patients. J. Neuroeng. Rehabil. 7:35. doi: 10.1186/1743-0003-7-35

PubMed Abstract | CrossRef Full Text | Google Scholar

de Vlugt, E., de Groot, J. H., van der Heijden-Maessen, H. C., Wielheesen, D. H., van Wijlen-Hempel, R. M. S., Arendzen, J. H., et al. (2013). Differentiation between non-neural and neural contributors to ankle joint stiffness in cerebral palsy. J. Neuroeng. Rehabil. 10:81. doi: 10.1186/1743-0003-10-81

PubMed Abstract | CrossRef Full Text | Google Scholar

Dietz, V. (2003). Spastic movement disorder: what is the impact of research on clinical practice? J. Neurol. Neurosurg. Psychiatr. 74, 820–821. doi: 10.1136/jnnp.74.6.820-a

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleuren, J. F., Voerman, G. E., Erren-Wolters, C. V., Snoek, G. J., Rietman, J. S., Hermens, H. J., et al. (2010). Stop using the Ashworth Scale for the assessment of spasticity. J. Neurol. Neurosurg. Psychiatry. 81, 46–52. doi: 10.1136/jnnp.2009.177071

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuglevand, A., Winter, D., and Patla, A. (1993). Models of recruitment and rate coding organization in motor-unit pools. J. Neurophysiol. 70, 2470–2488. doi: 10.1152/jn.1993.70.6.2470

PubMed Abstract | CrossRef Full Text | Google Scholar

Gäverth, J., Eliasson, A., Kullander, K., Borg, J., Lindberg, P., and Forssberg, H. (2014). Sensitivity of the NeuroFlexor method to measure change in spasticity after treatment with botulinum toxin A in wrist and finger muscles. J. Rehabil. Med. 46, 629–634. doi: 10.2340/16501977-1824

PubMed Abstract | CrossRef Full Text | Google Scholar

Gäverth, J., Sandgren, M., Lindberg, P., Forssberg, H., and Eliasson, A. (2013). Test-retest and inter-rater reliability of a method to measure wrist and finger spasticity. J. Rehabil. Med. 45, 630–636. doi: 10.2340/16501977-1160

PubMed Abstract | CrossRef Full Text | Google Scholar

Golberg, D. E. (1989). Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley Professional.

Harlaar, J., Becher, J., Snijders, C., and Lankhorst, G. (2000). Passive stiffness characteristics of ankle plantar flexors in hemiplegia. Clin. Biomech. 15, 261–270. doi: 10.1016/S0268-0033(99)00069-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Hesse, S., Krajnik, J., Luecke, D., Jahnke, M. T., Gregoric, M., and Mauritz, K. H. (1996). Ankle muscle activity before and after botulinum toxin therapy for lower limb extensor spasticity in chronic hemiparetic patients. Stroke 27, 455–460. doi: 10.1161/01.STR.27.3.455

PubMed Abstract | CrossRef Full Text | Google Scholar

Hidler, J. M., and Rymer, W. Z. (1999). A simulation study of reflex instability in spasticity: origins of clonus. IEEE Trans. Rehabil. Eng. 7, 327–340.

PubMed Abstract | Google Scholar

Jankovic, J. (2004). Botulinum toxin in clinical practice. J. Neurol. Neurosurg. Psychiatry 75, 951–957. doi: 10.1136/jnnp.2003.034702

PubMed Abstract | CrossRef Full Text | Google Scholar

Ju, M., Chen, J., Lee, H., Lin, T., Lin, C., and Huang, Y. (2000). Time-course analysis of stretch reflexes in hemiparetic subjects using an on-line spasticity measurement system. J. Electromyogr. Kinesiol. 10, 1–14. doi: 10.1016/S1050-6411(99)00018-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Katz, R., and Rymer, W. (1989). Spastic hypertonia: mechanisms and measurement. Arch. Phys. Med. Rehabil. 70, 144–155.

PubMed Abstract | Google Scholar

Koo, T. K., and Mak, A. F. (2006). A neuromusculoskeletal model to simulate the constant angular velocity elbow extension test of spasticity. Med. Eng. Phys. 28, 60–69. doi: 10.1016/j.medengphy.2005.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lance, J. (1980). The control of muscle tone, reflexes, and movement: Robert Wartenberg lecture. Neurology 30, 1303–1313. doi: 10.1212/WNL.30.12.1303

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H.-M., Chen, J.-J. J., Wu, Y.-N., Wang, Y.-L., Huang, S.-C., and Piotrkiewicz, M. (2008). Time course analysis of the effects of botulinum toxin type a on elbow spasticity based on biomechanic and electromyographic parameters. Arch. Phys. Med. Rehabil. 89, 692–699. doi: 10.1016/j.apmr.2007.08.166

PubMed Abstract | CrossRef Full Text | Google Scholar

Leis, A., Kronenberg, M., Stetkarova, I., Paske, W., and Stokic, D. (1996). Spinal motoneuron excitability after acute spinal cord injury in humans. Neurology 47, 231–237. doi: 10.1212/WNL.47.1.231

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindberg, P., Gäverth, J., Islam, M., Fagergren, A., Borg, J., and Forssberg, H. (2011). Validation of a new biomechanical model to measure muscle tone in spastic muscles. Neurorehabil. Neural Repair. 25, 617–625. doi: 10.1177/1545968311403494

PubMed Abstract | CrossRef Full Text | Google Scholar

Lorentzen, J., Grey, M., Crone, C., Mazevet, D., Biering-Sørensen, F., and Nielsen, J. (2010). Distinguishing active from passive components of ankle plantar flexor stiffness in stroke, spinal cord injury and multiple sclerosis. Clin. Neurophysiol. 121, 1939–1951. doi: 10.1016/j.clinph.2010.02.167

PubMed Abstract | CrossRef Full Text | Google Scholar

Manca, M., Merlo, A., Ferraresi, G., Cavazza, S., and Marchi, P. (2010). Botulinum toxin type A versus phenol. A clinical and neurophysiological study in the treatment of ankle clonus. Eur. J. Phys. Rehabil. Med. 46, 11–18.

PubMed Abstract | Google Scholar

Mirbagheri, M., Barbeau, H., Ladouceur, M., and Kearney, R. (2001). Intrinsic and reflex stiffness in normal and spastic, spinal cord injured subjects. Exp. Brain Res. 141, 446–459. doi: 10.1007/s00221-001-0901-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Morita, H., Petersen, N., Christensen, L., Sinkjær, T., and Nielsen, J. (1998). Sensitivity of H-reflexes and stretch reflexes to presynaptic inhibition in humans. J. Neurophysiol. 80, 610–620. doi: 10.1152/jn.1998.80.2.610

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozcakir, S., and Sivrioglu, K. (2007). Botulinum toxin in poststroke spasticity. Clin. Med. Res. 5, 132–138. doi: 10.3121/cmr.2007.716

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandyan, A. D., Vuadens, P., van Wijck, F. M., Stark, S., Johnson, G. R., and Barnes, M. P. (2002). Are we underestimating the clinical efficacy of botulinum toxin (type A)? Quantifying changes in spasticity, strength and upper limb function after injections of Botox® to the elbow flexors in a unilateral stroke population. Clin. Rehabil. 16, 654–660. doi: 10.1191/0269215502cr536oa

CrossRef Full Text | Google Scholar

Pauri, F., Boffa, L., Cassetta, E., Pasqualetti, P., and Rossini, P. M. (2000). Botulinum toxin type-A treatment in spastic paraparesis: a neurophysiological study. J. Neurol. Sci. 181, 89–97. doi: 10.1016/S0022-510X(00)00439-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Pennati, G. V., Plantin, J., Borg, J., and Lindberg, P. G. (2016). Normative NeuroFlexor data for detection of spasticity after stroke: a cross-sectional study. J. Neuroeng. Rehabil. 13, 30. doi: 10.1186/s12984-016-0133-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Phadke, C. P., Ismail, F., and Boulias, C. (2012). Assessing the neurophysiological effects of botulinum toxin treatment for adults with focal limb spasticity: a systematic review. Disabil. Rehabil. 34, 91–100. doi: 10.3109/09638288.2011.591882

PubMed Abstract | CrossRef Full Text | Google Scholar

Pingel, J., Wienecke, J., Lorentzen, J., and Nielsen, J. B. (2016). Botulinum toxin injection causes hyper-reflexia and increased muscle stiffness of the triceps surae muscle in the rat. J. Neurophysiol. 116, 2615–2623. doi: 10.1152/jn.00452.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Powers, R., Marder-Meyer, J., and Rymer, W. (1988). Quantitative relations between hypertonia and stretch reflex threshold in spastic hemiparesis. Ann. Neurol. 23, 115–24. doi: 10.1002/ana.410230203

PubMed Abstract | CrossRef Full Text | Google Scholar

Prochazka, A., and Gorassini, M. (1998). Models of ensemble firing of muscle spindle afferents recorded during normal locomotion in cats. J. Physiol. 507, 277–291. doi: 10.1111/j.1469-7793.1998.277bu.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheean, G. (2002). The pathophysiology of spasticity. Eur. J. Neurol. 9, 3–9. doi: 10.1046/j.1468-1331.2002.0090s1003.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheykh, M., Aliaei, G. R., Abou, A., and Asghar, H. H. (2005). Evaluation of plantar flexor muscles spasticity after the BTXA injection in hemiplegic patients using electrophysiological tests. Pejouhandeh 10, 191–195.

Google Scholar

Stampacchia, G., Bradaschia, E., and Rossi, B. (2004). Change of stretch reflex threshold in spasticity: effect of botulinum toxin injections. Arch. Ital. Biol. 142, 265–273. doi: 10.4449/aib.v142i3.375

PubMed Abstract | CrossRef Full Text | Google Scholar

Thacker, B. E., Tomiya, A., Hulst, J. B., Suzuki, K. P., Bremner, S. N., Gastwirt, R. F., et al. (2012). Passive mechanical properties and related proteins change with botulinum neurotoxin a injection of normal skeletal muscle. J. Orthop. Res. 30, 497–502. doi: 10.1002/jor.21533

PubMed Abstract | CrossRef Full Text | Google Scholar

Thelen, D. G. (2003). Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. J. Biomech. Eng. 125, 70–77. doi: 10.1115/1.1531112

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Helm, F., Schouten, A., de Vlugt, E., and Brouwn, G. (2002). Identification of intrinsic and reflexive components of human arm dynamics during postural control. J. Neurosci. Methods 119, 1–14. doi: 10.1016/S0165-0270(02)00147-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Herman, P., Ekeberg, Ö., Gäverth, J., Fagergren, A., and Forssberg, H. (2017). Neural and non-neural related properties in the spastic wrist flexors: an optimization study. Med. Eng. Phys. 47, 198–209. doi: 10.1016/j.medengphy.2017.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Herman, P., Ekeberg, Ö., Gäverth, J., Fagergren, A., Forssberg, H. A., et al. (2015). “A neuromusculoskeletal model to simulate passive joint extension test of spasiticy,” in Proceedings of the The XXV Congress of the International Society of Biomechanics 2015 (Glasgow).

Wilson, L., Gandevia, S., Inglis, J., Gracies, J., and Burke, D. (1999). Muscle spindle activity in the affected upper limb after a unilateral stroke. Brain 122, 2079–2088. doi: 10.1093/brain/122.11.2079

PubMed Abstract | CrossRef Full Text | Google Scholar

Winters, J., and Stark, L. (1985). Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models. IEEE Trans. Biomed. Eng. 32, 826–839. doi: 10.1109/TBME.1985.325498

PubMed Abstract | CrossRef Full Text | Google Scholar