Neuro-musculoskeletal simulation of instrumented contracture and spasticity assessment in children with cerebral palsy

Cerebral palsy (CP) is the most common neurological disorder in children and is attributed to non-progressive disturbances occurring in the developing fetal or infant brain [1]. It is characterized primarily by neural deficits (caused by the brain anomalies) and secondary by muscular and bone deformities [2, 3]. These deficits adversely affect normal development of functional activities, such as gait. Muscle hyper-resistance, i.e. increased tension in a passively stretched muscle, caused by neural impairments (spasticity and increased background activation) as well as by passive tissue stiffness, is the most prevalent problem. Spasticity is commonly defined as a velocity-dependent increase in tonic stretch reflexes due to hyper‐excitability [4], while stiffness is a mechanical resistance of the myotendinous tissue as it is passively lengthened [5].

The exact pathophysiology of hyper-resistance is complex and only partially understood. Specifically, there is lack of consensus on the mechanisms of, and interaction between, the neural and non-neural components of hyper-resistance. For example, it is generally thought that muscle stiffness occurs secondary to spasticity, with prolonged activation causing short and stiff fibers [6]. Compared with muscles in typically developing (TD) children, muscles in CP are smaller in cross sectional area, volume and muscle belly length, while Achilles tendons may be longer [7, 8]. However, studies on passive fascicle length and pennation angles have been inconclusive [9, 10], with some results indicating normal fascicle lengths [11] and long sarcomeres in spastic muscles [12]. Recent studies also suggest that muscle stiffness in CP already appears at an early age [13] and that maladaptation to growth, rather than spasticity, plays a crucial role in developing contractures [14].

The medial hamstrings are frequently affected by hyper-resistance in CP. In comparison to TD muscles, spastic hamstrings have been found to have lower activation thresholds during passive stretch [15] as well as altered muscle properties, including increased muscle stiffness due to larger amounts of extracellular matrix [16]. Recently, instrumented clinical tests to assess hyper-resistance in the hamstrings have been developed, with which both the biomechanical and electrophysiological responses are recorded during a passive stretch [17, 18]. Such tests yield a vast array of information about joint resistance, joint angles, angular velocity, and muscle electromyography (EMG). However, (instrumented) clinical assessment provides limited information at the underlying muscle level. This makes it difficult to distinguish between the neural and non-neural components of hyper-resistance, based on experimental data alone.

(Neuro-)musculoskeletal modeling is a powerful tool to gain insight into the underlying pathology at a more detailed level [19]. Thus far, most modeling studies of spasticity during passive stretch have focused on the ankle joint and have applied either simplified [20, 21], or highly complex [22–25] muscle models. Other modeling approaches have attempted to directly [26, 27] or indirectly [28] assess the effect of stretch reflexes on gait. However, the complexity of pathological muscle behavior during voluntary activation makes it very challenging to validate such models with measured data. A detailed analysis of spastic hamstring behavior during passive instrumented clinical assessment using neuro-musculoskeletal modeling could yield valuable information on tissue as well as reflective muscle behavior, but is so far lacking in literature.

Therefore, the aim of this paper was to simulate instrumented, clinical spasticity assessment of the hamstring muscles in CP using a conceptual model of contracture and spasticity. Specifically, we investigated whether we could explain 1) increased resistance during slow passive stretch (i.e. contracture) by altered passive muscle stiffness; and 2) additional increased resistance during fast passive stretch (i.e. spasticity) by a purely velocity-dependent stretch reflex.

Eleven children with spastic CP (6 male, age 11.5 ± 3.4y; weight 33.7 ± 12.8 kg) and 9 TD children (5 male; age 11.0 ± 3.2 y; 34.8 ± 13.5 kg) were included in the study. Of the CP subjects, 8 were bilaterally affected and 3 unilaterally; 8 were classified as Gross Motor Function Classification System (GMFCS) level I and 3 as II. All data were collected as part of a larger study [29]. Patients were selected from the available data set if the Modified Ashworth Scale score of the hamstring muscles was 1 or higher [30], i.e. a clinical sense of mild hyper resistance (1 subject had a score of 1; 6 of 1+, and 4 of 2; with an average Modified Tardieu Score of −78.5 ± 6.7°). Only data on left legs were included for practical reasons. To exclude patients with high background activation or who were not relaxed, measurements with a root mean square EMG (RMS-EMG) activity during slow passive stretch higher than 10 % of that measured during maximum voluntary contraction (about 2 % of the stretches) were excluded. Available data of 9 TD children were used for comparison [17]. All subjects older than 11 years and all parents signed an informed consent form. The data collection protocol was approved by the medical ethics committee of the KU Leuven University Hospital.

Instrumented spasticity assessment of the left leg hamstrings was performed on all subjects (Fig. 1a; [17]). The knee was stretched at a slow (>5 s for the entire range of motion) and a fast (<1 s) stretch velocity, while the subject was laying supine with the left hip in 90°, and the pelvis and thigh fixed by a second examiner. Knee angular displacement and velocity were recorded with inertial measurement units. The force applied to the shank by the examiner was recorded with a hand-held 6 degrees of freedom (DOF) force/torque sensor. The location of force application relative to the knee joint axis was measured manually parallel and perpendicular to the tibia. Surface EMG was collected from the medial hamstrings and rectus femoris (see [17] for measurement details).

A dedicated musculoskeletal model (Fig. 1b) was developed using OpenSim software [31]. The model was adapted from the generic gait model (GAIT2392) and included torso, pelvis, and left and right femur, tibia and foot (of which torso and right leg for visualization only). All joints were locked in 0°, except for the left hip, which was locked in 90° as imposed during all measurements, and the evaluated left knee, which was free to move until full extension. All muscles were removed except for those around the left knee, which were lumped to 5 major muscle groups: hamstrings (HAM, representing semitendinosus, semimembranosus and biceps femoris long head), vasti (VAS, representing vastus lateralis, medialis and intermedius), rectus femoris (RF), biceps femoris short head (BFS), and gastrocnemius (GAS, representing gastrocnemius lateralis and medialis). The path of the lumped muscles and the total muscle tendon length (optimal fiber length plus tendon slack length) were equal to that of one representative original muscle (biceps femoris long head for HAM, vastus intermedius for VAS, and gastrocnemius lateralis for GAS). Paths were selected based on other simplified models (e.g. Gait10dof18musc.osim) and previous literature [32]. To obtain myotendinous behavior representative of the whole muscle group, the fiber to tendon length ratio was averaged over all three muscles. The optimal muscle force was taken as the sum of the represented muscles’ forces. The ‘default activation’ was set to 0.01 for all muscles.

The model was scaled to individual subject sizes based on the subject’s height, leg length and tibia length, as measured during clinical examination. Muscle strength was scaled with the subject’s mass to the ²/₃-power. Experimentally measured knee angles and applied forces were imported to OpenSim. The measured 6-DOF force was applied to the shank segment at the same location relative to the knee as measured experimentally. An inverse dynamic (ID) analysis was run in OpenSim to obtain the net knee moment, representing the sum of all (passive plus active) internal muscle moments around the knee.

As a model for contracture, passive muscle properties for HAM and VAS were optimized to match the slow-velocity net knee moment versus knee angle curve. The other three muscles were found to have negligible effect, as they were not in their end range of motion in any part of the assessment. The two parameters that determine the passive force-length curve in the Thelen muscle model [33] were adjusted, i.e. the parameters FmaxMuscleStrain (‘S’, the passive muscle strain at maximum isometric muscle force) and KshapePassive (‘K’, the exponential shape factor for the passive force-length relationship). An fminsearch optimization algorithm was used (Matlab 2010, The Mathworks), with the default values in OpenSim (S = 0.6; K = 4) as starting values, and no boundaries. To test the robustness of the optimization, additional different starting values were used, and this was found to have negligible effect. A muscle analysis was performed for each set of S and K, to obtain the total muscle knee moments generated by all five muscles. The sum of all muscle moments was compared with the net knee moment from the ID analysis. In the optimization, the RMS error between the two curves was minimized. The optimal S and K values for HAM and VAS, as well as the associated RMS error were obtained for each subject.

To test how well the model could predict measured motion, forward dynamic (FD) simulations were performed for all assessments. The measured initial position and applied forces were used as input, while predicted knee angle and muscle activity were output of the simulations. The passive muscle properties S and K were set to the optimal values for each individual for all FD simulations. First, the slow stretch data were simulated for both TD and CP, to verify that the model was capable of replicating data using FD simulations. Second, the fast stretch data were simulated to predict the effect of stretch reflexes in CP, and to obtain a reference in TD.

with E _m the muscle excitation for muscle m, t _d a time delay factor representing the stretch reflex delay, G a gain factor representing the severity of the enhanced reflex, v _m the fiber velocity of muscle m, and T the threshold velocity. The controller was implemented as a plug-in in OpenSim and is freely available on https://simtk.org/home/spasticitymodel [34]. t _d was set to a fixed value of 30 ms, representing the shortest possible stretch reflex delay as found in the literature [35]. Based on the measured data, threshold values were determined as the stretch velocity 30 ms before onset of EMG of the muscle-tendon complex (thus assuming that muscle fiber velocity equals muscle-tendon velocity, which is reasonable considering the low forces). EMG onset was identified using the method of Staude and Wolf [36]. Gain values were set to 0, 1, 2 and 4 for each subject, which was found to cover the range of potential values. As no notable EMG activity was seen in any of the TD subjects’ fast stretches, the stretch reflex controller was only applied during CP fast stretch simulations.

To validate the FD results, predicted knee angles were compared with measured knee motion, and muscle excitation as evoked by the stretch reflex controller was compared with measured EMG. Both comparisons were done qualitatively by looking at the graphs, and for knee motion also quantitatively by calculating the RMS difference between the curves. A distinction was made between the stretch phase, in which the knee was moved quickly to its end range of motion, and the hold phase, in which the knee was held stationary at this end range of motion. In total 1 s of data was simulated for all FD simulations. The muscle fiber length and velocity during the CP fast stretch simulations were also analyzed, to gain further insight in the working mechanisms of the stretch reflex controller.

For statistical analysis, the optimal S and K parameters and the RMS error of the optimization were compared between TD and CP using a Mann–Whitney U test. This non-parametric test was chosen because of the small sample and non-normally distributed data. The RMS errors of the FD simulations for knee angle were also compared between TD and CP using a Mann–Whitney U test, and between different gain values in CP using a Kruskal-Wallis test. The level of significance was set to p < 0.05.

BFS, Biceps femoris short head; CP, Cerebral palsy; DOF, degrees of freedom; EMG, Electromyography; FD, Forward dynamics; GAS, Gastrocnemius lateralis and medialis; GMFCS, Gross Motor Function Classification System; HAM, Hamstrings muscle group; ID, Inverse dynamics; K, KshapePassive, exponential shape factor for the passive force-length relationship; PIC, persistent inward currents; RF, Rectus femoris; RMS, Root mean square; S, FmaxMuscleStrain, passive muscle strain at maximum isometric muscle force; TD, Typically developing; VAS, Vasti muscle group