In-vivo Identification of Skeletal Muscle Dynamics with Nonlinear Kalman Filter -Comparison between EKF and SPKF

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. In-vivo Identification of Skeletal Muscle Dynamics with Nonlinear Kalman Filter -Comparison between EKF and SPKF Mitsuhiro Hayashibe, David Guiraud, Philippe Poignet


Introduction
1.1.Muscle Modeling and FES.Functional electrical stimulation (FES) is an efective technique to evoke artiicial contractions of paralyzed skeletal muscles.It has been employed as a general method in modern rehabilitation to partially restore motor function for patients with upper neural lesions [1,2].Recently, the rapid progress in microprocessor technology provided the means for computer-controlled FES systems [3][4][5].A fundamental problem concerning FES is how to handle the high complexity and nonlinearity of the neuromusculoskeletal system [6,7].In addition, there is a large variety of patient situations depending on the type of neurological disorder.To improve the performance of motor neuroprosthetics beyond the current limited use of such system, subject-speciic modeling is essential.he use of a mathematical model can improve the development of neuroprosthetics by optimizing their functionality for individual patients.
A mathematical model makes it possible to describe the relevant characteristics of the patient's skeletal muscle a n dt oa c c u r a t e l yp r e d i c tt h ef o r c ea saf u n c t i o no ft h e stimulation parameters.Indeed, synthesis of stimulation sequences or control strategies to achieve movement can be eiciently computed and optimized using numerical models [8].herefore, it can contribute to enhancing the design and function of controls applied to FES.Over the years, a great variety of muscle models have been proposed, difering in their intended application, mathematical complexity, level of structure considered, and idelity to the biological facts.Some of them have attempted to exhibit the microscopic or macroscopic functional behavior, for instance Huxley [9]and Hill [10]. he distribution-moment model [11]c o n s t i t u t e sa bridge between the microscopic and macroscopic levels.It is a model for sarcomeres or whole muscle which has been extracted via a formal mathematical approximation from Huxley crossbridge models.Models integrating geometry of the tendon and other macroscopic consideration can be found in [12].A study by Bestel and Sorine [13], based on both microscopic Huxley and macroscopic Hill type model, proposed an explanation of how the beating of cardiac muscle is achieved through a chemical control input.It integrated the calcium dynamics in muscle cells that stimulate the contractile element of the model.Starting with this concept, we adapted it for striated muscle [14].We proposed a musculotendinous model which considered the muscular masses and viscous frictions in muscle-tendon complex.his model is represented by diferential equations where the outputs are the muscle's active stifness and force.he model input represents the actual electrical signal as provided by the stimulator in FES.
1.2.In Vivo Identiication for Subject-Speciic Parameters.In actual FES system, the appropriate tuning is achieved empirically by intensively stimulating the patient's muscle for each task.If this adjustment could be calculated in the simulation, and if we could ind the best signal pattern using virtual skeletal muscle, such method would be very helpful for movement synthesis for spinal cord injured (SCI) patients.However, in order to perform the simulation, an accurate skeletal muscle model is required to reproduce a well-predicted force for each muscle corresponding to the patient-speciic characteristics.
For any biological systems, identiication is a diicult problemduetothefactthat(i)measurementsmustbeasnoninvasive, particularly on humans, (ii) some entities cannot be directly measured, (iii) an experimental setup and protocol have to be designed and certiied, (iv) intersubject variations can be large, and (v) the large nonlinearity and complexity of the models cause some optimization algorithms to fail.hus, few papers address biomechanical parameters identiication in FES context, and they used macroscopic model for global force production [15].Consequently, we described an approach for coupling the model with in-vivo measurements, that is, using a multiscale muscle model in an estimation procedure in order to perform the identiication of the parameters, hence, giving access to physiologically meaningful parameters of the muscle.Preliminary result was reported in [16].
he skeletal muscle dynamics are in particular highly nonlinear, and we need to identify many unknown physiological parameters if a multiscale model is applied.he main objective of this paper is to develop an experimental computational method to identify unknown internal parameters from the limited information.For such identiication of nonlinear system, extended Kalman ilter (EKF) h a so t e nb e e nu s e df o rs k e l e t a lm u s c l e [ 17]a n dc a r d i a c muscle identiication [18].It can work when the model is not complicated and not highly nonlinear.In this paper, a sigma-point Kalman ilter (SPKF) was applied to the in vivo experimental data to identify internal states in the nonlinear dynamics of multiscale skeletal muscle model.SPKF has higher accuracy and consistency for nonlinear estimation than EKF theoretically. he feasibility of both identiications is veriied by comparison. he computational performance is discussed.he outline of the paper is as follows.he next section presents the formulation of the skeletal muscle model controlled by FES. he following section is devoted to the experimental identiication of the model for isometric contraction, including the identiication protocol.he experimental measurement was performed in-vivo on rabbits.hen, we present detailed results of the parameter estimation, the comparison with EKF, and cross-validation which illustrates the pertinence of the identiication.Finally, we present some discussions, conclusions, and perspectives in the last sections.

Skeletal Muscle Model
Muscle modeling is complex, in particular when the model is based on biomechanics and physiological realities.Most of the muscle models have been based on phenomenological models derived from Hill's classic work [10]a n dw e l ls u mmarized by Zajac [12].Hill macroscopic model is a standard muscle model for practical use.Recent work [19]performed the validation of Hill model during functionally relevant conditions.hey concluded that model errors are large for diferent iring frequencies and largest at the low motor unit iring rates relevant to normal movement.hey pointed out that the reason may come from the Hill model assumption to consider muscle activation, force-length, and forcevelocity properties independently.It was suggested that more physiological coupling between activation and force-velocity properties can be demonstrated in microscopic crossbridge models incorporating a dependence between physiologically based activation and crossbridge attachment [19].hus, our approach is to provide a multiscale physiology-based model on the both micro-and macroscale fact to obtain meaningful internal parameters.
Our muscle model is composed of two elements of diferent nature: (i) an activation model which describes how an electrical stimulus generates an action potential (AP) and initiates the contraction and (ii) a mechanical model which describes the dynamics in force (Figure 1).For details of the muscle model, refer to the previously published article [14].Here, a brief summary of the model and the information necessary for the identiication are given.
2.1.Activation Model. he activation model describes the electrical activity of muscle which represents the excitationcontraction phenomena.In muscle physiology, it is known that two input elements dominate muscle contraction: iber recruitment rates and temporal activation [20]. he recruitment rates determine the percentage of recruited motor units.In Hill model, it has only one input which corresponds to this recruitment rate.Hill model is mainly applied for voluntary contraction where temporal activation is not dominant.However, in FES, since all muscle ibers are synchronously activated, this temporal activation which occurs by every stimulation pulse, is important.he recruitment rates can be determinedonthepulsewidth(PW)andpulseamplitudeof signal , generated by the stimulator [21]. he recruitment rate (PW,)canbeassumedtobeaconstantvaluewhenPWand remain constant.he temporal activation can be considered as the underlying physiological processes which describes the chemical input signal, , that brings muscle cells into contraction as shown in Figure 2. Muscle contraction is initiated by an AP along the muscle iber membrane, which goes deep into the cell through the T-tubules.It causes calcium releases that induce the contraction process when its concentration rises aboveathresholdandissustainedtilltheconcentrationdrops back below this threshold once again.Hatze [22]g i v e sa n example of calcium dynamics (Ca 2+ ) modeling.Since we focus on the mechanical response in this paper, we choose a simple model that renders the main characteristics of the dynamics.he contraction-relaxation cycle is triggered by the (Ca 2+ ) associated with two phases: (i) contraction and, (ii) relaxation as in Figure 2.W eu s ead e l a y e d( )m o d e l to take into account the propagation time of the AP and an average time delay due to the calcium dynamics.he frequency of the chemical input, , can be deined from the stimulation frequency.he time delay and the contraction timecanbeobtainedfromatwitchtestbysinglestimulation pulse.A contraction takes place with a kinetics then, if no other AP has been received in the mean time, an active relaxation follows indeinitely with a kinetics . is linked t ot h er a t eo fa c t i n e -m y o s i n ec y c l ew h e r e a s is related to the rate of crossbridge breakage.Finally, can be written as follows, where Π () is a trapezoidal switching function which connects relaxation and contraction state from 0 to 1: 2.2.Mechanical Model. he model is based on the macroscopic Hill-Maxwell type model and the microscopic description of Huxley [9].For crossbridge modeling, Huxley [9] proposed an explanation of the crossbridge interaction in a sarcomere.A sarcomere model can be used to represent a wholemusclewhichisassumedtobeahomogeneousassembly of identical sarcomeres. he distribution-moment model of Zahalak [11] is a model for sarcomeres or whole muscle which is extracted via a formal mathematical approximation from Huxley crossbridge models.his model constitutes a bridge between the microscopic and macroscopic levels.B a s e do nH u x l e ya n dH i l l -t y p em o d e l s ,B e s t e la n dS o r i n e [13] proposed an explanation of how the beating of cardiac muscle may be performed, through a chemical control input, c o n n e c t e dt ot h ec a l c i u md y n a m i c si nm u s c l ec e l l ,t h a t stimulates the contractile element of the model.

Stim. PW Contraction Relaxation Contraction Relaxation
he model is composed of macroscopic passive elements and a contractile element controlled by input commands: the chemical input, , as suggested by Bestel and Sorine [13] for the cardiac muscle on the sarcomere scale and the recruitment rate, , on the iber scale as shown in Figure 3.In order to express isometric contractions, whereas the skeleton is not actuated, our muscle model is introduced with masses (kg) and linear viscous dampers 1 , 2 (Ns/m) to ensure energy dissipation.On both sides of ,t h e r ea r ee l a s t i c springs, 1 , 2 (N/m), and viscous dampers to express the viscoelasticity of the muscle-tendon complex. he parallel element, , mainly represents surrounding tissues, but it can be omitted in isometric contraction mode.We assume the s y m m e t r i cf o r mw i t ht h et w om a s s e sa n dp a s s i v ee l e m e n t s identical.Here, 0 and 0 (m) are the lengths of and in the rest state.Initially, we can deine the relative length variation as positive when the length increases, as in (2).In particular, in the case of isometric contraction, the following relationship exists (3): he dynamical equation of one of the masses is given by (4).and express the force and the stifness of ,respectively.
he force a tt h eo u t p u to ft h i sm u s c l em o d e li st h es u m of the spring force and the damping force .W h e n we measure the tension of skeletal muscle under in vivo conditions, the experimentally measured force is equivalent to .Whenwetaketheratioof and , 0 is ofset and can bewrittenasin( 6).( 7) is the relational equation in Laplace transform.From this relationship, the diferential equation ( 8)canbeobtained: Finally, in isometric contraction, the diferential equations of this model can be described as follows: he dynamics of the contractile element correspond to ( 9) and (10). he terms and are the maximum values for and ,respectively .From(3)and( 4), the diferential equation for is obtained as in (12). he internal state vector of this system should be set as x =[ ].

Experimental Identification
In this study, we have developed a method to identify the parameters in the mechanical part of a skeletal muscle model.he input controls of the model are set as the static recruitment rate and the chemical control input from the activation model.hese two controls are computed from the FES input signal.Note that the identiication was performed with constant FES parameters for pulse width and intensity of electrical stimulation so that the recruitment rate is constant.he amplitude and pulse width were selected to recruit 90% of the maximum muscular force; then can be set as 0.9 for the iber recruitment.In addition, the calcium dynamics in our model induce a time delay and an on/of control which represents contraction/relaxation so that correct data processing can avoid the detailed calcium dynamics.he trigger for signal can be calculated from the timing of the electrical stimulation.and are set as 20 and 15, respectively, as in [23].
In isometric contraction, the diferential equations of skeletal muscle dynamics are straightly given in ( 9)- (12).In this case, , , ,a n d are unknown time-varying values and , ,a n d 0 are unknown static parameters to be estimated.For the identiication of this model, it is a nonlinear state space model, and many state variables are not measurable.Moreover, in-vivo experimental data include some noise.Hence, we need an eicient recursive ilter that estimates the state of a dynamic system from a series of noisy measurements.
3.1.Sigma-Point Kalman Filter.For this kind of nonlinear identiication, extended Kalman ilter (EKF) is the wellknown standard method.In EKF, the nonlinear equation should be linearized to irst order with partial derivatives (Jacobian matrix) around a mean of the state.he optimal Kalman iltering is then applied to the linearized system.When the model is highly nonlinear, EKF may give particularly poor performance and diverge easily.In skeletal muscle dynamics, its state-space is dramatically changed between the contraction and relaxation phases.At this time, partial derivatives would be incorrect due to the discontinuity.herefore, we introduced the sigma-point Kalman ilter (SPKF).he initial idea was proposed by Julier and Uhlmann [24] and has been well described by Merwe and Wan [25].SPKF uses a deterministic sampling technique known as the unscented transform to pick a minimal set of sample points (called sigma points) around the mean.hese sigma points arepropagatedthroughthetruenonlinearity .hisapproachresults in approximations that are accurate to at least second order in a Taylor series expansion.In contrast, EKF results in only irst-order accuracy.
An outline of the SPKF algorithm is described.For details, the reader should refer to [25,26]. he general Kalman framework involves estimation of the state of a discrete-time, nonlinear dynamic system, where x represents the internal state of the system to be estimated and y is the only observed signal.he process noise k drives the dynamic system, and the observation noise is given by n .he ilter starts by augmenting the state vector to dimensions, where is the sum of dimensions in the original state, model noise, and measurement noise. he corresponding covariance matrix is similarly augmented to an by matrix.In this form, the augmented state vector x and covariance matrix P canbedeinedas where P x is the state covariance, Q k is the process noise covariance, and R n is the observation noise covariance.
In the process update, the 2 + 1 sigma points are computed based on a square root decomposition of the prior covariance as in (15), where = √ +and is found using = 2 ( + ) − . is chosen in 0<<1which determines the spread of the sigma points around the prior mean and is usually chosen equal to 0. he augmented sigmapoint matrix is formed by the concatenation of the state sigma-point matrix, the process noise sigma-point matrix, and the measurement noise sigma-point matrix, such that .hesigma-pointweightsto be used for mean and covariance estimates are deined as in (16). he optimal value of 2 is usually assigned to : And these sigma points are propagated through the nonlinear function.Predicted mean and covariance are computed as in (18)a n d ( 19) and predicted observation is calculated as follows: he predictions are then updated with new measurements by irst calculating the measurement covariance and statemeasurement cross-correlation matrices, which are then used to determine the Kalman gain.Finally, the updated estimate a n dc o v a r i a n c ea r ed e t e r m i n e df r o mt h i sK a l m a ng a i na s below: hese process updates and measurement updates should be recursively calculated in = 1,...,∞until the end point of the measurement.

Experimental Measurement for Identiication.
Stimulation experiments were performed on two New Zealand white rabbits at Aarhus University Hospital in Aalborg, Denmark, as depicted in Figure 4. Anesthesia was induced and maintained with periodic intramuscular doses of a cocktail of 0.15 mg/kg Midazolam (DormicumR, Alpharma A/S), 0.03 mg/kg Fetanyl, and 1 mg/kg Fluranison (combined in HypnormR, Janssen Pharmaceutica) [14]. he let leg of the rabbit was anchored at the knee and ankle joints to a ixed mechanical frame using bone pins placed through the distal epiphyses of the femur and tibia.A tendon of the medial gastrocnemius (MG) muscle was attached to the arm of a motorized lever system (Dual mode system 310B Aurora Scientiic Inc.). he position and force of the lever arm were recorded.An initial muscletendon length was established by lexing the ankle to 90 ∘ .A bipolar cuf electrode was implanted around the sciatic nerve, allowing the MG muscle to be stimulated.Data acquisition was performed at a 4.8 kHz sampling rate. he muscle f o r cea ga in s ttheelectricals tim ula tio nwasmeas ur edunder isometric conditions.

Result of Identification
In order to facilitate the convergence of the identiication, the estimation process has been split into two steps.In the irst step, we only estimate geometrical parameter 0 .Inthe second step, we estimate the dynamic parameters: and .he biomechanical muscle model to be identiied is presented as ( 9)- (12).We deine the state vectors for geometric and dynamic estimations in SPKF as follows: successive pulses (doublet) at 16 Hz and 20 Hz with amplitude 105 Aa n dp u l s ew i d t h3 0 0s. he stifness has been estimated separately, using the experiments performed on the isolated muscle.he isolated muscle was pulled with the motorized lever. he stifness is taken as equal to the slope of thestraightlineofthepassivelengthversusforcerelationship in the experimental measurement as shown in Figure 5.he obtained was 4500 N/m. and can be obtained knowing the force response of muscle to a stimulation pattern with maximum signal parameter values (frequency, amplitude, pulse, and width).In this case, = 1000 N/m and =15N were obtained from the measurement.he experimental muscle force against the doublet stimulation (20 Hz) was used for parameter estimation during measurement updates for in SPKF. he covariance matrix has been initialized as in Table 1 for both geometric and dynamic parameter identiications.he evolutions of estimated internal states for and are obtained as shown in Figure 6.From these behaviors, we can conirm that the contractile element of the model is successfully tracking the contraction represented by the dynamics of diferential equations under the estimation process.Figure 7 shows the estimated parameter 0 and the error covariance. he fact that the error covariace is being decreased during identiication process shows stable convergence of the estimation.Even with the diferent initial values setting, parameter estimation results converged into a particular value in SPKF as shown in Figure 7(a).We tested the estimation from 7 diferent initial values with an interval of 5mm from 55mm to 85mm. he color is changed from magenta to cyan in the plot.
As can be seen from the resultant computational behavior on the graphs, the internal state vectors of the muscle model converged well to stationary values even with diferent initial values.Ater the complete estimation process for geometric and dynamic parameters, the estimated values are summarized in Table 2. Normalized RMS errors between measured and estimated force were also computed for each case.Moreover, the estimated length of the contractile element with each stimulation showed values close to the actual measured lengths of the extracted muscle.he border between the contractile element and tendon is not so clear visually, but it was approximately 7 cm for both rabbit GM muscles. he sizes of the two rabbits were similar; in particular, the estimated length of 0 h a dg o o dc o r r e s p o n d e n c ew i t ht h e measured length for both rabbits and with diferent frequency stimulations.Physical parameters are impossible to make visual veriication, but the estimated intrasubject property is maintained among the parameters obtained by diferent frequency stimulations from Table 2.
Initial covariance matrix for (dynamic) 4.2.Comparison with the Extended Kalman Filter. he extended Kalman ilter is generally chosen for nonlinear system identiication.However, in EKF, irst-order partial derivatives are used for the computation, which means that a matrix of partial derivatives (Jacobian) is computed around the estimate for each step. he detail of EKF is summarized in the appendix.When the process and measurement functions f and h are highly nonlinear, EKF can give particularly poor performance and diverge easily [25]. he estimate can have a bias due to the linear approximation especially for discontinuous systems.
Using the same computational conditions, such as initial values for states, parameters, and covariance matrices, parameter estimation was also performed with EKF to compare the estimation quality for this nonlinear system. he estimation results with EKF are shown in Figures 8 and 9. he numerical comparison for SPKF and EKF in 7 diferent initial conditions is summarized in Table 3 both for rabbit1 and rabbit2.With the observation noise covariance =2.5 × 10 −3 ,thesameaswithSPKF ,theresultwithEKFconverged computationally as shown at the error covariance of EKF1 in thetable.However,thematchedresultinforceleveldoesnot mean directly that the estimated internal state and parameter are correct.he internal state in EKF estimation was not well estimated as the contracted strain variation.he diference can be observed by comparing the plots of Figures 9 and 6 for the EKF and SPKF estimates, respectively.In Figure 6, the estimated state of has two sharp peaks which correspond to the strain of the contractile element induced by muscle contraction by doublet stimulation.hus, this internal state relects the expected phenomenon in muscle dynamics.However, the estimated state in Figure 9 does not show sharp deformation induced by stimulation pulse.
he linear approximation in the EKF transformation m a t r i xc o u l db et h ec a u s eo ft h el o w e rq u a l i t yo ft h e state estimation.Finally, it resulted in a bias for parameter estimation.hus, the converged value for 0 under these conditions is not realistic.he EKF estimation was also performed with the observation noise covariance = 2.5 × 10 −2 , giving more uncertainty to the measurement update.In this case, the estimated value became more realistic as the SPKF estimate did.However, the convergence became poor as in Figure 8(b) and it was highly dependent on the initial values as in Figure 8(a).

Model Cross-Validation.
A cross-validation of the identiied model was carried out to conirm the validity of this method.he resultant muscle force was calculated using the identiied models with the control input generated by the stimulation frequency.Figure 10 shows the measured force response of the MG muscle of the rabbit and the simulated force with the identiied model.he stimulation was three successive pulses in = 105 A, PW = 300 s, and Freq = 3 1 .2 5H z .her e dl i n ei n d i c a t e st h em e a s u r e dm u s c l ef o r c e , theblueandgreendottedlinesaretheplotsbytheidentiied model with SPKF and EKF1, respectively.he model identiied by SPKF could predict the nonlinear force properties of stimulated muscle quite well in this cross-validation.Figure 11 Average errors (%) 2.9 44.9 47.9 EKF1 and EKF2 represent EKF estimation with the observation noise covariance = 2.5 × 10 −3 and = 2.5 × 10 −2 ,r especti v el y .represents steady-state error covariance.
is the result in the case of stimulation trains with six pulses in 20 Hz. here is a good agreement between the measured and the predicted force; however, when muscle fatigue appears, the experimental force is lower than that from the simulation.Particualrly, the diference between the measured and the estimated forces at the end can be considered as the error coming from modeling without considering fatigue.Apart from the error coming from muscle fatigue, we could conirm the feasibility of the FES muscle modeling and the efectiveness of the identiication.Normally, the muscle property varies greatly, being dependent on the type of muscle, and the force response is highly subject-speciic.herefore, the muscle force prediction by cross-validation is quite diicult even for the approximate prediction especially when force production is predicted for each stimulation pulse. he identiication based on experimental response would contribute strongly to realistic force prediction in electrical stimulation and FES controller development.Further investigation is required for the expression of fatigue phenomena [27], for its reproducibility and identiication.

Discussion
he skeletal muscle model used in this identiication protocol is based on both macroscopic and microscopic physiology which is unlike the black-box or other approaches using simple Hill muscle model.he structured model requires more parameters in highly nonlinear dynamics, which have to be estimated by experiment.However, the great advantage of our model is the insight which it can give to a muscle's biomechanical and physiological connections, where the parameters have a physical signiicance such as length and mass.his paper describes an identiication method which uses experimental response and a nonlinear, physiological muscle model to obtain subject-speciic parameters. he advantage of animal experiment is to have direct access to extracted muscle for conirmation of the obtained parameter ater the experiment.his work was performed for application of FES stimulated muscles; however, since many living systems have nonlinear dynamics and subject speciicity, this kind of identiication approach itself could be applied to other organ models and clinical situations.In this work, both SPKF a n dE K Fa l g o r i t h m sw e r ea p p l i e df o rm u s c l ed y n a m i c s identiication.EKF showed a high dependency on initial settings of the computation and it was diicult to ind the efective range of initial states and covariances.In SPKF, t h eo b t a i n e dr e s u l tw a sm o r ec o n s i s t e n ta n dr o b u s tw i t h respect to various initial conditions.Moreover, for SPKF,  the computation of a Jacobian is not necessary, so it can e a s i l yb ea p p l i e de v e nt oc o m p l e xd y n a m i c s .S P K Fr e s u l t s in approximations that are accurate to at least second order in Taylor series expansion.In contrast, EKF results in irstorder accuracy.Further, the identiication accuracy is clealy improved, especially for nonlinear systems, but the computation cost still remains the same as for EKF.Advanced and robust system identiication, including designing experimental protocols, has a very important role to improve the control issues in neuroprosthetics.In addition, the main diiculty in understanding human systems is caused by their timevarying properties. he systems are not static and change over time. he function of a human being is not always the same; for example, muscle fatigue can easily change the expected force response.In order to deal with the timevarying characteristics of a human system, robust biosignal processing and model-based control which corresponded to nonlinearity and time variance would provide a breakthrough in the development of neuroprosthetics [28].

Conclusion
We have developed an experimental identiication method for subject-speciic biomechanical parameters of a skeletal muscle model which can be employed to predict the nonlinear force properties of stimulated muscle.he mathematical muscle model accounts for the multiscale major efects occurring during electrical stimulation.hus, the identiication method was required to deal with the high nonlinearlity and discontinous states between muscle contraction a n dr e l a x a t i o np h a s e . he performance was compared and summarized under thesamecomputationalconditions.SPKFgivesmorestable performance than EKF. he internal state in EKF estimation was not well estimated as the contracted muscle strain.
In SPKF estimation, keeping the realistic state transition and independence from initial conditions, it could realize converged solution for each identiication trial.SPKF is a Bayesian estimation algorithm which recursively updates the posterior density of the system state as new observations arrive online.his framework can allow us to calculate any optimal estimate of the state using newly arriving information.We believe that the proposed identiication method has also the advantage for human muscle identiication while it provides not only better accuracy in nonlinear dynamics but also adaptability to time-varying systems. he preliminary result is reported as in [29].

Figure 1 :
Figure 1: Outline of skeletal muscle model and its identiication.

Figure 3 :
Figure 3: Macroscopic mechanical coniguration of the muscle model.

Figure 7 :
Figure 7: Estimated parameter for the original length of the contractile element 0 (a) and its error covariance (b) with SPKF. he plot color is changed from magenta to cyan with 7 diferent initial values from 55 mm to 85 mm.

Figure 8 :Figure 9 :
Figure 8: Estimated parameter for the original length of the contractile element 0 (a) and its error covariance (b) with EKF. he plot color is changed from magenta to cyan with 7 diferent initial values from 55 mm to 85 mm. he results with two diferent settings for the observation noise covariance are shown.

Figure 10 :
Figure 10: Measured and simulated isometric muscle force by the identiied models with three successive stimulation pulses (in = 105 A, PW = 300 s, and Freq = 31.25Hz).
h ei d e n t i i e dm o d e lw a se v a l u a t e d by comparison with experimental measurements in crossv a l i d a t i o n .he r ew a sag o o da g r e e m e n tb e t w e e nm e a s u r e d a n ds i m u l a t e dm u s c l ef o r c eo u t p u t s .her e s u l ts h o w e d the performance which can contribute to the prediction of the nonlinear force of stimulated muscle under FES. he feasibility of the identiication could be demonstrated by comparison between the estimated parameter and the meas u r e dv a l u e .I nt h i ss t u d y ,t h ei d e n t i i c a t i o nw a sp e r f o r m e d by sigma-point Kalman ilter and extended Kalman ilter.