Extra excitation of biceps femoris during neuromuscular electrical stimulation reduces knee medial loading

Medial knee joint osteoarthritis (OA) is a debilitating and prevalent condition. Surgical treatment consists of redistributing the forces from the medial to the lateral compartment through osteotomy, or replacing the joint surfaces. As the mediolateral load distribution is related to the action of the musculature around the knee, the aim of this study was to devise a technique to redistribute these forces non-surgically through changes in muscle excitation. Eight healthy subjects participated in the experiment, and neuromuscular electrical stimulation was used to change the muscle forces around the knee. A musculoskeletal model was used to quantify the loading on the medial compartment of the knee, and a novel algorithm devised and implemented to simulate neuromuscular electrical stimulation. The forces and moments at the knee, ground reaction forces, walking velocity and step length were quantified before and after stimulation. Stimulation of the biceps femoris resulted in a significant decrease in the second peak of the medial knee joint loading by up to 0.17 body weight (p = 0.016). Kinematic parameters were not significantly affected. Neuromuscular electrical stimulation can decrease the peak loads on the medial compartment of the knee, and thus offers a promising therapy for medial knee joint OA.


Introduction
Osteoarthritis (OA) is a major degenerative disease, the prevalence of which is predicted to increase significantly due to a growing & 2019 The Authors. Published by the Royal Society under the terms of the Creative ageing population [1]. The highest prevalence of OA is in the medial compartment of the knee [2], and excessive contact force is known to be a risk factor for the onset and progression of medial knee OA [3]. Therefore, reducing the contact force on the medial compartment of the knee is a tantalizing opportunity to treat medial knee OA.
It has been shown that gait modification by varying the external forces on the knee through valgus bracing and lateral wedges can decrease the knee adduction moment (KAM) and/or shift the centre of pressure of the ground reaction force (GRF) laterally [4][5][6][7]. These measures are surrogates of the medial knee contact forces (KMF) and the efficacy of these approaches is highly dependent on the individual subject [8,9].
As the forces acting upon the knee comprise external forces that are counterbalanced by the internal forces (reaction forces between the bony segments, muscle forces and ligament forces) [10], and internal forces are a redundant system with multiple different muscle contraction patterns able to provide this counterbalance, so it is conceivable that by changing the internal muscle forces, the mediolateral load distribution can be modified without changing the external forces. This distribution is highly related to the action of the hamstrings, quadriceps and gastrocnemius muscles [11]. For example, it is known that an elevated KAM is directly correlated with increased medial loading, and that quadriceps and gastrocnemius offer the most resistance to KAM [12]. Therefore, an alternative to orthotic gait modifications is to directly change the muscle forces during gait.
Muscles forces can be changed by altering the timing and strength of contraction through, for example, physiotherapy. An alternative approach that is used for quadriceps muscle strengthening for knee OA patients, when they present with chronic pain and joint stiffness [13,14], is neuromuscular electrical stimulation (NMES). NMES is a non-invasive technology that induces a voltage gradient along axons that can depolarize membranes that then induce action potential, which causes muscle contraction. NMES of the quadriceps has been shown to modify knee moments [15] and could be used to reduce KMF by selecting the appropriate muscles for stimulation. Muscles known to realign the distribution of the medial and lateral knee forces are the long head of biceps femoris (loBF), lateral gastrocnemius (latGAS) and vastus lateralis (VL) [11]. Therefore, in this study, it is hypothesized that the selective excitation of these muscles will redistribute the joint loading from the medial to the lateral side in healthy subjects, and thus, NMES may be used to ameliorate symptoms in sufferers of medial knee joint OA.
KMF can be measured directly through the use of implanted prostheses [16], partially inferred through the KAM [17 -21], are indirectly related to knee flexion moment (KFM) [21,22], and can be reliably estimated through the use of musculoskeletal modelling [23][24][25][26]. Of these, the only method suitable for non-invasive direct quantification of KMF in those without arthroplasty is musculoskeletal modelling. Some musculoskeletal models are driven by electromyography (EMG) measurements; however, NMES stimulus contaminates EMG [27] meaning that this approach cannot be used. Other methods are based on static optimization [28][29][30][31] that use cost functions that minimize muscle fatigue [32], defined as the sum of cubed muscle stresses. However, when NMES is used to externally drive muscle recruitment, the cost function is not representative. Therefore, in this study, a modified static optimization cost function is proposed to estimate KMF in NMES-assisted gait.
The aim of this study was to test the hypothesis that extra excitation of the lateral leg muscles (i.e. loBF, latGAS and VL) can reduce KMF during level walking, where KMF is quantified through the use of the newly proposed cost function in a validated musculoskeletal model [23]. This is then compared with the commonly used surrogate of KMF, the KAM, and finally, the input kinematic and kinetic data are analysed to explain the results.

Experimental data
Eight healthy subjects (age, 27.3 + 2.4 years; mass, 62.7 + 14.4 kg; height, 1.69 + 0.11 m; three males and five females) participated in this study which consisted of performing four gait tasks, three trials each of: normal walking, and three walking tasks in which separate muscles were stimulated individually (loBF, latGAS and VL).
The skin overlying the muscles was prepared with 70% isopropyl alcohol skin wipes, and the electrode pads (PALS w Platinum, Axelgaard Manufacturing Co., Ltd, Fallbrook, CA, USA) were placed at the two motor points of the relevant muscles ( figure 1a). An asymmetrical biphasic waveform was used with royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181545 the maximum current of 60 mA (1 kV). The frequency of the pulse was set to 40 Hz as recommended by the NMES manufacturer (Odstock 2 Channel Stimulator, Odstock Medical Ltd, UK). The pulse width was adjusted from 132.5 to 306.5 ms to generate visible muscle contraction at the limit of the subjects' tolerance. The stimulation was started manually by the subject by pressing a hand-held switch during the left foot strike before stepping on the force plate. The stimulation lasted 2 s to cover the entire stance phase of the right limb.
Right lower limb motion was recorded at 200 Hz using a nine-camera VICON motion analysis system (Vicon Motion Systems Ltd, Oxford, UK). Eighteen reflective markers were used in this study: RASIS, LASIS (right and left anterior superior iliac spine); RPSIS, LPSIS (right and left posterior superior iliac spine); T1, T2, T3 (cluster markers on the thigh); FLE, FME (lateral and medial femoral epicondyle); C1, C2, C3 (cluster markers on the shank); FAM, TAM (apex of the lateral and medial malleolus); FCC (calcaneus); FMT (tuberosity of the fifth metatarsal); FM2 (head of the second metatarsal) and TF (the centre of the acrotarsium). The system was calibrated using a 5 marker L frame.
GRFs were recorded at 1000 Hz from a force plate (Kistler Type 9286BA, Kistler Instrument AG, Winterthur, Switzerland). During stimulated walking, only the right leg muscles were stimulated, and the force plate recorded the force applied by the right foot only.

Modelling
The flowchart (figure 2) outlines the data flow and modelling process.
The open source musculoskeletal model Freebody 2.0, which is based on the original model of Cleather & Bull [33], was used [23]. The validation of this original model has been described previously [23,[33][34][35][36] through comparison between the model calculations and in vivo measurements of joint reaction forces or EMG. The revised model with NMES has been validated [36] through the measurement of gluteus maximus activation using EMG pre-and post-application of NMES on loBF;    the kinematics and kinetics input to the model were measured simultaneously with the EMG, and strong positive correlations were found in early stance peak (R ¼ 0.78, p ¼ 0.002) and impulse (R ¼ 0.63, p ¼ 0.021). The original anatomical model used is from Klein et al. [37] with 163 muscle elements. However, in order to improve the relevance of the dataset to a clinical population, and improve on the joint contact point data, the anatomical geometry were newly derived from MR imaging data of one younger male (height: 1.83 m; mass: 96 kg; age: 43) selected from a database of anatomical models [38]. The anatomical geometry included the muscle origins, via points, and insertion points, along with joint centres of rotation, and contact points between the femur and tibial plateau.
In the local coordinate system of each segment (foot, shank, thigh and pelvis), the x-axis is in the anterior-posterior direction with anterior being positive; the z-axis is in the medial-lateral direction with lateral being positive and the y-axis is in the longitudinal direction of the segment (figure 1b). The anatomical model consists of the muscle insertion points of 163 muscle elements. These are scaled to each subject using scaling factors, where a muscle point in the anatomical model has coordinates (x 0 , y 0 , z 0 ); this is multiplied by a scaling factor to obtain the coordinates of the muscle point in the subject (x, y, z) (equation (2.1)). Scaling factors, sf, of the segments ( pelvis, thigh, shank and foot) were determined from segmental ratios as in equation (2.2) The model has five segments in total: foot, shank, thigh, pelvis and patella. A segment-based inverse dynamics calculation was conducted to quantify the intersegmental forces and moments [33], and an optimization-based approach was used to quantify the muscle forces [28,34]. The inverse dynamics method of Dumas et al. [39] was used to formulate the Newton -Euler equations of motion that include the segmental weight, external forces, muscle forces and joint reaction forces. The muscle forces were constrained to be between 0 and a maximal force, which was calculated by multiplying the physiological cross-sectional area (PCSA) from Klein's research by an assumed maximum muscle stress of 3.139 Â 10 5 N m 22 [37]. Muscle forces of 163 muscle elements, one patella ligament force and five joint reaction forces (ankle, medial knee, lateral knee, hip and patellofemoral) were considered.
The cost function of the static optimization is the sum of the cubic relative muscle forces, which was proposed based on the force -endurance relationship [32] (below equation) where F j is the muscle force of muscle element j, and F jmax is the maximal muscle force of muscle element j. The joint reaction forces except the tibiofemoral forces were considered to act through the centres of rotation of the joints taken from MR imaging data: the hip joint centre was defined as the centre of the femoral head, which was determined by fitting a sphere to the femoral head; the knee joint centre was defined at the midpoint of the central axis of a cylinder that was fitted to the boundaries of the femoral condyles; and the ankle joint centre was defined at the centre of a sphere that was fitted to the talar dome. The centres of rotation were fixed to the distal segment once determined. However, as the tibiofemoral joint reaction force is compartmentalized into a medial and a lateral component, the contact points of these two compartments are different from the knee joint centre of rotation, and they were located according to the MR imaging data of the subject mentioned above. The contact points were fixed to the shank segment and defined as the points of action of the medial and lateral knee joint forces; these were the middle of the most distal ends of the femoral condyles in the axial slice. These positions will change on the thigh segment according to the translations between two adjacent segments. The subjects' contact points were obtained by scaling these two points in the MR imaging data to the subject's shank using the shank scaling factor.
royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181545 The outputs of the optimization are: magnitude of 163 muscle forces, one patella ligament force and five vectors of joint reaction forces for the ankle, lateral knee, medial knee, hip and patellofemoral joints.

NMES simulation
The optimization described above quantifies muscle forces according to an endurance function [32]. However, stimulated muscles are expected to have larger forces than in the normal case. In order to adapt this model to NMES-assisted walking, a new parameter c was added to obtain the new cost function as in the below equation [40] X 163 where c ¼ c s for stimulated muscles 1 others: ð2:5Þ The adjustable weighting factor, c s , enables the stimulated muscle force to increase relative to standard simulation for a value of c s of less than 1. Conversely, increasing c s to greater than 1 will result in under-activation of the muscle. As electrical stimulation recruits muscles directly, the combination of voluntary contractions with electrical stimulation produces stronger contraction [41] and so, clinically, c s is related to stimulation intensity and it is used to produce changes in muscle forces. In order to explore the sensitivity of c s on the muscle and joint forces, varying values of c s (1, 0.75, 0.5, 0.25, 0.1 and 0) were implemented. Finally, c s ¼ 0.1 for BF and c s ¼ 0.25 for latGAS and VL were selected to produce a significant increase in muscle force during the NMES simulation.
The new cost function (NMESsim) drives the optimization to produce larger forces for the stimulated muscles, and thus allows the simulation of NMES-assisted walking.
The joint reaction forces in normal gait were calculated with the original cost function in equation (2.3). As the imposition of muscle contraction using NMES could potentially alter the gait kinematics and kinetics, the effect of NMES on these variables only was assessed first using the original cost function. To investigate the effect of increasing muscle activation of the stimulated muscles, NMESsim in equations (2.4) and (2.5) was also used to quantify NMES-assisted walking. Joint kinematics were calculated as Euler angles based on the local coordinate systems of the adjacent segments [39]. The angles were referenced to the standing posture at which the joint angles were defined to be zero. The sequence of angle calculation is first flexion/extension ( plantar flexion/dorsiflexion for the ankle), then external/internal rotation and lastly abduction/adduction.

Statistical analysis
The Wilcoxon signed-rank test (Matlab, R2014a, The Mathworks, Inc., 2014) was used to determine if there was a significant difference between normal walking and NMES-assisted walking. Muscle forces were compared to validate the effectiveness of NMESsim. As two different cost functions were used, and three muscles were stimulated, the loBF, latGAS and VL muscle forces in normal gait were compared with those forces during stimulated walking using both the original cost function and NMESsim.
The magnitude (in body weight-BW) of KMF and KAM were compared between cases. The ankle, knee and hip joint angles were calculated and compared for normal and NMES-stimulated walking tasks.
The two peak values of GRF and KFM, mean walking velocity and step length during the stance phase were compared between normal gait and NMES-stimulated walking. In all analyses, each of the three trials per subject were analysed, and then mean parameters are presented.
Anonymized data files from the physical experiments and the computational simulations are available on request from the corresponding author.

Sensitivity analysis
The muscle forces and KMFs with varying c s for eight subjects are plotted in figure 3. The stimulated muscle forces (including loBF, latGAS and VL forces) increase from the turquoise line (c s ¼ 1) to the royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181545 blue line (c s ¼ 0) as c s decreases (figure 3a). The first peak of KMF decreases as the loBF force increases, and the second peak of KMF decreases as the latGAS force increases, from the turquoise line to the blue line. The VL has a similar effect as loBF on the first peak, but the decrease in KMF only became clearly apparent with c s ¼ 0. Figure 4 depicts the mean forces of the main lower limb muscles of eight subjects using NMESsim during four tasks. The nine muscles are loBF, latGAS, VL, medial head of GAS (medGAS), soleus (SL), tibialis anterior (TA), rectus femoris (RF), gluteus maximus (GlutMax) and gluteus medius (GlutMed).

Muscle forces
There was no significant difference between maximal muscle force values in normal gait and stimulated walking when using the original cost function. However, maximal muscle forces in the stimulated tasks were significantly larger than the maximal muscle forces in normal gait when using NEMSsim for stimulated walking (loBF: 0. 27

Knee forces and moments
The mean KMFs are shown in figure 5. Generally, the KMF has two peaks during the stance phase occurring at the mean values of 28 and 78% of the stance phase. However, two subjects do not have two explicit peaks. Therefore, the mean KMFs of 20-35 and 70 -85% of the stance phase are calculated to represent the two peak values. The second peak KMF in loBF-stimulated walking is decreased significantly using the original cost function (

Ground reaction forces, walking velocity and step length
GRF is presented in figure 7. There was no statistical difference in peak GRFs, walking velocity and step length between NMES stimulated and normal walking (table 5). The peak GRFs were defined as those at the same timing as the peak KMFs. Specifically, these were calculated as the mean values across 20-35 and 70-85% of the stance phase.

Discussion
This study has shown that the musculoskeletal system can be tuned by increasing muscle forces to reduce medial knee joint contact forces (decreased peak values of KMF) in healthy subjects through the excitation of loBF, and has tested this using NMES. This opens up the possibilities of using NMES to treat knee OA, a major clinical condition causing disability and pain. This study is not in agreement with previous computational studies that found that no muscles that cross the knee could reduce the KMF [24,42]. A key difference between these studies and the present study is that they assumed that the muscle activations were calculated without kinematic compensations, i.e. the kinematic and kinetic datasets remained the same for normal and stimulated  walking. They did not stimulate the muscles, and did not incorporate the changing kinematics and kinetics due to NMES. One of the advantages of the present study is that the kinematics and kinetics caused by increased muscle activations were used in the calculations. The different results of these studies highlight the importance of quantifying the effect of stimulation on kinematics and kinetics.

The use of c s
Voluntary contractions preferentially recruit type I fibres, and then type II fibres. This asynchronous activation of varied motor units is a physiological mechanism that decreases muscle fatigue [41]. However, electrically elicited contractions recruit motor units synchronously, and recruit motor units   [41]. It is deduced that stimulated muscles produce larger forces than voluntary conditions, and the other muscles contract in a way to decrease muscle fatigue. In view of this, NMESsim was used with a parameter c to realize a reasonable muscle force distribution: the stimulated muscle force was increased, and all the muscle and joint forces satisfied the motor function constraints. The c s values were chosen between 0 and 1 (0, 0.1, 0.25, 0.5, 0.75 and 1), where c s ¼ 1 represents the original cost function. Where the muscles are not saturated (i.e. limited by their PCSA), assigning a value of c s of less than 1 enforces a greater muscle force than when using the original cost function; a value of c s ¼ 0 assigns the greatest activation to the muscle as this removes the muscle from the cost function. Using c s ¼ 0.1 for BF and c s ¼ 0.25 for latGAS and VL increased the activation of the stimulated muscles (table 1), but did not, in this study, saturate the muscle (one of the constraints requires the muscle force to be less than the maximal muscle force). Therefore, theoretically, the simulated muscle force with c s ¼ 0.1 and c s ¼ 0.25 is expected to be practically achievable. The physical experiment used NMES to contract the muscles at the maximum comfortable level for each subject.  The new algorithm, NMESsim, succeeded in simulating the stimulated condition, because the modelpredicted muscle forces under stimulation were significantly larger than those of normal cases and all the muscle forces and joint reaction forces satisfied the equations of motion.
EMG signals are affected by the stimulation artefacts caused by the stimulation current when NMES is used. This has been partially addressed in the literature with, for example, use of an EMG-amplifier with shut-down control [43], advanced filtering procedures [44] and optimizing the positioning of the EMG electrodes in relation to the stimulation electrodes [45]. However, these approaches do not fully

KMF and KAM
The shape of the KMF in stance phase has an expected double hump pattern with slightly higher peak than measured in vivo [19 -21,26,46,47]. The literature has much information on the reasons for higher joint force predictions using musculoskeletal modelling; in this case, we hypothesize that such forces are appropriate due to our subjects being younger and more active than those subjects in the literature who had total joint replacements. The decrease in KMF due to NMES is not substantial compared to the reduction caused by other methods. Kinney et al. [47] quantified KMF reduction through gait modification using long hiking poles with wide pole placement, and found that KMF at 75% stance decreased significantly from 1.35 BW to 0.89 BW. Other modifications in Kinney's study, such as mild crouch and medial thrust, did not show any significant KMF reduction. Orthotic interventions using a lateral foot wedge and valgus knee brace resulted in a decrease in model-estimated peak tibiofemoral load with increasing wedge size and applied valgus brace moment [7]. A higher level of activation using NMES might be required to create a greater reduction in KMF.
Although some of variables were not significantly different between groups, this does not negate the fact that changes can be significant for an individual. Table 3 lists how peak KMF varies for eight subjects and although the KMF decrease for the group is not substantial, the simulation could be used to identify which patient is well suited to this intervention.
The KAM curves shown in figure 6 are similar to those of Manal et al. [21]. The KAM is reported to be highly correlated with the medial contact force for different gait patterns, and is thus often used as a surrogate measure [19]. The significant decrease in KAM by loBF stimulation was consistent with that of KMF. A study of 62 OA patients verified the clinical importance of KAM in that peak KAM and KAM impulse explained variance in knee cartilage thickness [48]. However, whether the change of KAM could contribute to the alleviation of OA symptoms needs clinical validation.

GRF, walking velocity and step length
Changes in walking velocity and step length have been found to affect knee joint contact forces [49]. The GRF, as the external force on the whole body, is related to the overall kinematics and kinetics. Here, there is no significant difference in GRF, velocity and step length, so although the NMES changed the subtle kinematics of the lower limb, the lack of overall changes in pace could explain the lack of variation in GRF.
Although normal gait is not perfectly symmetrical, stimulation of one side alone makes the gait even more asymmetrical. However, the results here show that these changes were small and not large enough to obviously disturb normal walking as the GRF, walking velocity and step length of stimulated walking are similar to that of normal walking.

Limitations and future work
This study simulated the situations where only one muscle was stimulated. However, it is possible that the stimulation of combined muscles will not change the net joint moment, but together have a greater impact on KMF than stimulating each muscle independently. This case should be considered further.
The muscle forces of loBF, latGAS and VL generate movement of the lower limb mainly through changing the flexion-extension moments of the knee and ankle. The extra muscle activation by NMES may influence the periodic characteristics of normal gait. In this way, different timing of NMES leads to different gait characteristics, and will affect the reduction in KMF. Therefore, NMES timing is an important parameter for precise control of NMES. This study used a very crude manually implemented timing mechanism that should be tuned for optimal outcome.
The decrease in the KMF was limited by the increase in the muscle force, which was determined by the c s value in the model. Lower c s values than were applied in this study could be used to generate a greater reduction in KMF. However, the clinical feasibility of including such high muscle forces is not known.
This study is a pilot study on the effect of NMES on the KMF, and only healthy subjects were considered. The study verified that the stimulation of loBF will redistribute the knee loading by royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181545 decreasing the KMF. The result provides a new way to control the knee loading distribution. These results suggest that it may also be possible to decrease the medial knee joint loading of OA patients using NMES. Therefore, the next steps in this work should involve the recruitment of OA patients in order to verify the effect of NMES on a pathological group.

Conclusion
A new cost function for use in musculoskeletal models, NMESsim, has been proposed to allow the quantification of muscle forces when they are stimulated by external means. This was tested in living subjects, and our results show that the stimulated muscle forces are significantly increased (loBF: D ¼ 0.07 BW, p ¼ 0.039; latGAS: D ¼ 0.27 BW, p ¼ 0.008; VL: D ¼ 0.40 BW, p ¼ 0.008) and the peak values of knee joint medial loading are significantly decreased by applying NMES to loBF (NMESsim: D ¼ 0.17 BW, p ¼ 0.016). This study demonstrates that it is possible to redistribute the knee loading and reduce the loading on the medial compartment by activating selected muscles across the healthy knee and this opens the door for prevention and alternative non-surgical interventions for knee OA.
Ethics. The experiment was approved by the ethics committee of Imperial College London (reference number: 14IC2134), and all subjects gave written informed consent.