MD M Schuur

PhD P Henneman

J C Swieten

PhD M C Zillikens

I Koning

PhD A C J W Janssens

PhD J C M Witteman

PhD Y S Aulchenko

PhD R R Frants

PhD B A Oostra

PhD J C Van Swieten

PhD I De Koning

PhD K Willems Van Dijk

Professor Cornelia M Van Duijn
email: c.vanduijn@erasmusmc.nl

Insulin-resistance and metabolic syndrome are related to

Keywords: cognition, executive function, metabolic syndrome, insulin-resistance, adiponectin

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Cognitive impairment is a common problem in the elderly and an important predictor of dementia (DSM-IV). There are multiple risk factors contributing to cognitive decline. For long, type 2 diabetes has been recognized as a major risk factor of dementia [START_REF] Ott | Diabetes mellitus and the risk of dementia: The Rotterdam Study[END_REF].

Epidemiological studies indicate that hyperinsulinemia and insulin-resistance, which characterize type 2 diabetes, may cause cognitive dysfunction [START_REF] Geroldi | Insulin resistance in cognitive impairment: the InCHIANTI study[END_REF][START_REF] Young | Hyperinsulinemia and cognitive decline in a middle-aged cohort[END_REF]. Insulin-resistance is also a feature of the metabolic syndrome (MetS) and there is increasing interest in the role of MetS in cognition. The syndrome is defined as the presence of at least 3 out of 5

factors, which include central obesity, hyperglycemia, hypertension, hypertriglyceridemia, and low high-density lipoprotein cholesterol (HDL-c) levels [START_REF] Grundy | Diagnosis and management of the metabolic syndrome: an American Heart Association/National Heart, Lung, and Blood Institute Scientific Statement[END_REF]. It has a high prevalence of over 22% in Western societies [START_REF] Ford | Prevalence of the metabolic syndrome among US adults: findings from the third National Health and Nutrition Examination Survey[END_REF], and is an important risk factor of cardiovascular disease [START_REF] Alberti | The metabolic syndrome--a new worldwide definition[END_REF][START_REF] Olijhoek | The metabolic syndrome is associated with advanced vascular damage in patients with coronary heart disease, stroke, peripheral arterial disease or abdominal aortic aneurysm[END_REF]. Potentially, MetS is an important, modifiable, determinant of cognitive pathology [START_REF] Kwon | Metabolic syndrome as an independent risk factor of silent brain infarction in healthy people[END_REF][START_REF] Kumari | Minireview: mechanisms by which the metabolic syndrome and diabetes impair memory[END_REF][START_REF] Roriz-Cruz | Cognitive impairment and frontalsubcortical geriatric syndrome are associated with metabolic syndrome in a stroke-free population[END_REF]. Previous studies have shown associations between cognitive dysfunction and MetS, but results are not consistent and predominantly concern elderly people [START_REF] Yaffe | The metabolic syndrome and development of cognitive impairment among older women[END_REF]. Of its individual components, hypertension and hyperglycemia have proven to be the most consistent determinants of cognitive function in a recent review of published studies [START_REF] Van Den Berg | Type 2 diabetes mellitus, hypertension, dyslipidemia and obesity: A systematic comparison of their impact on cognition[END_REF].

Adipose tissue also plays an important role in the etiology of MetS. Adipose tissue functions as endocrine tissue, among which secretion of adipocytokines and cytokines, which indirectly affect C-reactive protein (CRP) levels [START_REF] Devaraj | Adiponectin decreases Creactive protein synthesis and secretion from endothelial cells: evidence for an adipose tissue-vascular loop[END_REF][START_REF] Guerre-Millo | Adipose tissue and adipokines: for better or worse[END_REF][START_REF] Rader | Inflammatory markers of coronary risk[END_REF]. The visceral fat depot is thought to play an important role in adipocytokine secretion and MetS [START_REF] Bacha | Adiponectin in youth: relationship to visceral adiposity, insulin sensitivity, and beta-cell function[END_REF]. Several groups have suggested that in particular adiponectin might be a promising biomarker for MetS [START_REF] Gilardini | Adiponectin is a candidate marker of metabolic syndrome in obese children and adolescents[END_REF][START_REF] Luchsinger | Adiposity, hyperinsulinemia, diabetes and Alzheimer's disease An epidemiological perspective[END_REF][START_REF] Lara-Castro | Adiponectin and the metabolic syndrome: mechanisms mediating risk for metabolic and cardiovascular disease[END_REF]. Whether it is related to cognitive function in humans is not known, but animal studies indicate that adiponectin may influence brain metabolism [START_REF] Kubota | Adiponectin stimulates AMPactivated protein kinase in the hypothalamus and increases food intake[END_REF]. Regarding CRP levels, it has been reported that the association of MetS with cognitive function is stronger in those with high levels of CRP [START_REF] Yaffe | Metabolic syndrome and cognitive disorders: is the sum greater than its parts?[END_REF][START_REF] Dik | Contribution of metabolic syndrome components to cognition in older individuals[END_REF].

The MetS is an important health problem and may be a suitable target for early intervention. The aim of the current study was to examine whether the MetS and related factors, insulin-resistance, plasma adiponectin and CRP-levels, are related to cognitive function in a middle-aged population in the Netherlands.

Methods

Study population

The study was performed in the Erasmus Rucphen Family Study (ERF), which is an ongoing family-based cohort study in a genetically isolated population [START_REF] Aulchenko | Linkage disequilibrium in young genetically isolated Dutch population[END_REF]. The ERF population includes around 3000 individuals, all living descendents of 22 couples who, at the end of the 19 th century had at least six children baptized in the community church.

Extensive data on cardiovascular risk factors, body composition, cognitive functioning, and blood chemistry are available. These data were collected between 2002 and 2005. All participants gave informed consent and the study was approved by the Medical Ethics Committee at ErasmusMC.

Clinical and laboratory assessments

To define MetS, we used the definition proposed by the National Cholesterol Education

Program Adult Treatment Panel (ATPIII) [START_REF] Grundy | Diagnosis and management of the metabolic syndrome: an American Heart Association/National Heart, Lung, and Blood Institute Scientific Statement[END_REF]. Waist circumference was measured exactly halfway the distance between the lower costal margin and the iliac crest [START_REF] Santos | Heritability of fasting glucose levels in a young genetically isolated population[END_REF]. Blood pressure was measured twice with an automatic device in a sitting position [START_REF] Sleegers | Cerebrovascular risk factors do not contribute to genetic variance of cognitive function: the ERF study[END_REF].

Blood was taken after fasting of at least 8 hours. Insulin-resistance was determined as homeostasis model assessment insulin-resistance (HOMA-IR) [START_REF] Matthews | Homeostasis model assessment: insulin resistance and beta-cell function from fasting plasma glucose and insulin concentrations in man[END_REF], which was computed by multiplying fasting plasma insulin (micro-international units per milliliter) and fasting plasma glucose (millimoles per liter)/22. are asked to connect a randomly placed series of letters and numbers as quickly as possible [START_REF] Reitan | The relation of the trail making test to organic brain damage[END_REF]. In TMT-A only the numbers have to be connected, in TMT-B numbers and letters have to be connected alternately (from 1 to A, to 2 to B etc). For analysis, we used time in seconds depicted as the ratio of TMT-B/TMT-A. The Stroop Color and Word

Test is also a time demanding task consisting of three cards. In card I, participants have to name the right words as quickly as possible, in card II they have to name the right colors as quickly as possible and in card III they have to name the colors in which the words are printed as quickly as possible. On this card the meaning of the word is different than the color the word is printed in (e.g. blue is written whereas the word is printed in red) [START_REF] Hammes | Stroop Kleur-woord Test: Dutch Manual[END_REF].

For analyses, we used time in seconds depicted as the ratios of Stroop-III/Stroop-II. In the verbal fluency tests, participants had to name as many animals and as many words starting with the same letter (D, A and T), each in one minute. The score was defined as the total number of correct words on both tasks. Lastly, visuo-spatial ability was assessed with the WAIS-III block-design subtest. The test consists of reproducing two-dimensional patterns using cubes that have red, white, and half-red-half-white faces. The total number of correctly reproduced blocks within a time limit was used as a score.

In addition to analyzing individual test scores, three composite scores were used: 1.

memory (z-memory); 2. executive function (z-executive); 3. global cognition (z-global).

These scores were computed by taking the average of z-scores of the individual tests, which were computed by subtracting the mean from the test score divided by the standard deviation [START_REF] Prins | Cerebral small-vessel disease and decline in information processing speed, executive function and memory[END_REF]. For time-demanding tasks, the test score was subtracted from the mean.

The composite scores were only computed when at least 3 memory scores (z-memory), 2 executive scores (z-executive) and 6 total test scores (z-global) were available. In the current study, we excluded participants who had a history of stroke or dementia (N=21)

or were physically not able to perform a neuropsychological test (N=44). There were 10 participants with missing test scores due to technical problems.

Statistical analysis

One large pedigree was extracted from the genealogical database, which was cut into smaller pedigrees before analyses with the software PedCut [START_REF] Liu | An approach for cutting large and complex pedigrees for linkage analysis[END_REF]. In total, 1898 people with complete data could be linked to pedigrees with a maximum bit size of 18 and these were included in the analyses.

Baseline descriptive analyses were performed with SPSS version 15.0, using the independent-samples T-test (continuous variables) or Chi-square statistics (categorical data). C-reactive protein levels were transformed by taking the common logarithm.

Because high levels of CRP could have been caused by acute inflammation, all values of the log-transformed variable greater than three times the standard deviation above the mean were excluded.

Possible confounders were first analyzed as single covariates in two regression models with either MetS or cognitive test scores as dependent variables. Also, interaction terms of MetS, HOMA-IR, adiponectin and CRP with age and sex were analyzed in these models. Covariates that were significantly associated were then entered simultaneously into a model using backward regression. We removed all covariates with a p-value > 0.10. The final model included age, level of education, HADS-D score, alcohol use and inbreeding coefficient, which is a measure of relatedness of the parents of an individual.

The analyses were stratified by sex, because the interaction terms of MetS with sex were significant.

To identify associations between cognitive test scores and MetS, HOMA-IR, adiponectin, and CRP, we performed regression analyses using the SOLAR software version 4.1.0 [START_REF] Almasy | Multipoint quantitative-trait linkage analysis in general pedigrees[END_REF]. The MetS was included in the model as a binomial variable, and additionally as a continuous variable ranging from 0 to 5 to test the effect of increasing number of MetS components. We used quartiles of HOMA-IR, adiponectin and CRP with the lowest quartile as a reference category to estimate the association with adipose tissue endocrine function.

Results

The mean age of the study population was 47.7 (±14.3) years in women and 49.2 (±14.1) years in men, with a range of 18 to 86 years and less than 10% older than 67. Metabolic syndrome was present in 227 (21%) women and in 207 (25%) men. Of these, 155 women and 150 men fulfilled three of the MetS diagnostic criteria, 63 and 54 fulfilled four and 9 and 3 fulfilled all five. Null criteria were fulfilled by 282 women and 153 men. A comparison of the individual MetS components is given in Table 1.

Table 2 shows the baseline characteristics of men and women with and without MetS regarding determinants of cognitive function and vascular and endocrine factors.

We studied the associations between MetS and cognitive test scores. Women with MetS had significantly lower scores on Stroop en z-executive (p=0.03 and p=0.02) than women without MetS (Table 3). An increasing number of MetS components was associated to lower verbal fluency and z-executive scores (p = 0.04 and 0.02) in women. The same trend, however non-significant (p=0.10), was seen for Stroop.

The findings related to adipose tissue endocrine function are given in Figure 1 with analyses limited to executive function. Women with high HOMA-IR scores had lower scores on TMT, Stroop and z-executive than women with low HOMA-IR scores. Higher adiponectin levels were associated with higher scores on z-executive in both men and women and with Stroop in women only. The associations of HOMA-IR with Stroop and z-executive were more significant (p=0.004 and 0.007) when analyzing HOMA-IR as continuous variables instead of quartiles. The other associations, however, became nonsignificant. No association between CRP and executive function was found in this population.

Discussion

The main finding of the current study is that MetS and HOMA-IR is associated with poorer executive performance in women.

There a number of limitations of the study to be discussed before interpretation of the findings.The differences between persons with and without the MetS were relatively small, which makes the clinical relevance or our findings difficult to predict. Also, we performed many tests. Although statistically these tests are not independent, some of the significant p-values could be false-positive findings. When applying strict Bonferonni adjustment, our results are at best only borderline significant (p homa-ir =0.11). We did not a priori have an hypothesis which of the different cognitive outcomes were associated. The finding of a relation to executive function but not to memory remains to be interpreted with caution. Although some are more consistent that others, the consistent associations argue against false-positive findings. The cross-sectional design is another limitation of our study, making it difficult to assess causality. Several longitudinal studies, however, suggest that MetS might be a causative factor in the relationship with cognition [START_REF] Yaffe | The metabolic syndrome, inflammation, and risk of cognitive decline[END_REF].

Another point to be discussed is that we studied the relationship of MetS and cognition in a genetically isolated population. We have taken relatedness into account by adjusting for pedigree structure and relatedness, and since it is a recently isolated population, the genetic makeup does not differ largely from the general population and we do not expect large differences in the effects of cardiovascular risk factors on cognition [START_REF] Pardo | The effect of genetic drift in a young genetically isolated population[END_REF]. As to the generalizaibility, there is no reason to assume that the findings in terms of the etiological pathways uncovered should not aply to outbred populations. However, the associations may be more or less pronounced in this family-based study. The strength of this study is the population-based design in which participants were not selected on the disease of interest. Also, the inclusion of a large sample size and a range of cognitive tests allowed detailed study of cognition.

Our findings of association of MetS with cognitive performance confirms previous crosssectional and longitudinal studies [START_REF] Dik | Contribution of metabolic syndrome components to cognition in older individuals[END_REF][START_REF] Komulainen | Metabolic syndrome and cognitive function: a population-based follow-up study in elderly women[END_REF][START_REF] Yaffe | Metabolic syndrome and cognitive decline in elderly Latinos: findings from the Sacramento Area Latino Study of Aging study[END_REF]. Although we studied different cognitive domains, mainly executive function was associated with MetS. This is in line with two recent studies that found significant differences in processing speed, semantic fluency and executive function between persons with and without MetS [START_REF] Van Den Berg | Cognitive functioning in elderly persons with type 2 diabetes and metabolic syndrome: the Hoorn study[END_REF][START_REF] Segura | Mental slowness and executive dysfunctions in patients with metabolic syndrome[END_REF]. Associations of the MetS with fronto-subcortical syndrome and with frontal white matter changes were found by others [START_REF] Roriz-Cruz | Cognitive impairment and frontalsubcortical geriatric syndrome are associated with metabolic syndrome in a stroke-free population[END_REF][START_REF] Segura | Microstructural white matter changes in metabolic syndrome: a diffusion tensor imaging study[END_REF]. Damage to cerebral small vessels could be one of the explanations for the observed associations in the current study [START_REF] O'brien | Cognitive associations of subcortical white matter lesions in older people[END_REF][START_REF] Van Dijk | The association between blood pressure, hypertension, and cerebral white matter lesions: cardiovascular determinants of dementia study[END_REF][START_REF] Pa | Clinical-neuroimaging characteristics of dysexecutive mild cognitive impairment[END_REF] and also why we find an association with executive function and not memory. A number of studies find that an increasing number of MetS components is associated with lower cognitive scores, which is also confirmed by our observations [START_REF] Roriz-Cruz | Cognitive impairment and frontalsubcortical geriatric syndrome are associated with metabolic syndrome in a stroke-free population[END_REF][START_REF] Yaffe | Metabolic syndrome and cognitive decline in elderly Latinos: findings from the Sacramento Area Latino Study of Aging study[END_REF].

An effect of MetS components on the early development of Alzheimer's disease (AD) pathology may be another possible mechanism underlying our observations, since the insulin-degrading enzyme (IDE) functions less well in the presence of insulin-resistance, resulting in higher levels of amyloid β-peptide (Aβ) [START_REF] Farris | Partial loss-of-function mutations in insulin-degrading enzyme that induce diabetes also impair degradation of amyloid beta-protein[END_REF][START_REF] Craft | Insulin resistance syndrome and Alzheimer's disease: age-and obesityrelated effects on memory, amyloid, and inflammation[END_REF]. Long-term exposure to insulin due to insulin-resistance, may by itself also have a direct damaging effect on neurons [START_REF] Li | Hippocampal neuronal apoptosis in type 1 diabetes[END_REF]. Our observed association of HOMA-IR and executive function is in line with this.

Despite possible pathways through which adiponectin could be linked to cognition, such as improvement of insulin-sensitivity, anti-inflammatory effects and effects on brain metabolism and the vascular system, our results do not show a consistent association of adiponectin with cognition [START_REF] Han | Antiatherosclerotic and anti-insulin resistance effects of adiponectin: basic and clinical studies[END_REF][START_REF] Kubota | Adiponectin stimulates AMPactivated protein kinase in the hypothalamus and increases food intake[END_REF] . Adiponectin levels are however, lower in individuals with MetS in our study and a trend was seen with executive function in women, which supports further study on the role of adiponectin in cognition.

Levels of CRP were not related to cognition in the present study, however, data on other markers, such as interleukin-6 (Il-6) and α1-antichymotrypsin (ACT) were not available.

These markers have been related to cognitive decline and increased systemic inflammation has been proposed as a possible modifier between MetS and cognition [START_REF] Dik | Contribution of metabolic syndrome components to cognition in older individuals[END_REF].

We could not confirm this, which may also be due to a more pronounced effect of inflammation on cognition in the elderly than in middle aged individuals.

Our study shows that there are gender differences in the association of MetS with cognitive function. Gender-specific effects of metabolic components on cognitive function have not extensively been addressed in previous studies, while there is evidence that MetS has greater effect on cardiovascular and cerebrovascular disease in women than in men, indirectly supporting our observed differences [START_REF] Pischon | Inflammation, the metabolic syndrome, and risk of coronary heart disease in women and men[END_REF][START_REF] Mcneill | The metabolic syndrome and 11-year risk of incident cardiovascular disease in the atherosclerosis risk in communities study[END_REF]. Furthermore, there is evidence that men are treated more aggressively for cardiovascular risk factors [START_REF] Gouni-Berthold | Sex disparities in the treatment and control of cardiovascular risk factors in type 2 diabetes[END_REF]. As a consequence, men may be less likely to experience the adverse effects of MetS. Another factor that could explain gender differences is depression, which is highly correlated to cognition and more common in women. Since our analyses were adjusted for depressive symptoms, however, we do not think that depression explains the observed differences in our study.

In summation, the most important findings of our study are that MetS and HOMA-IR are associated with executive function in women. We interpret these results with caution due to the cross-sectional design and the potential problems related to multiple testing.

However, our study yields an important starting point to replicate these findings and take them further. Longitudinal studies will be needed to gain insight in the causality of our reported findings and may result in more conclusive findings.

Table 1 : Metabolic Syndrome in the Erasmus Rucphen Family study

 1

	(n=1898)

Table 2 : Metabolic Syndrome in relation to vascular and endocrine factors and determinants of cognitive function

 2

		Men		Women
		MetS absent (n=615) MetS present (n=207)	MetS absent (n=849)	MetS present
	Vascular and endocrine factors
	Smoking (%)	32.7	33.5	46.5	45.1
	Insulin (μU/ml, FPI)	11.8 (5.6)	19.1 (11.1) **	11.5 (4.7)	18.0 (9.2) **
	HOMA-IR	2.4 (1.4)	4.5 (2.9) **	2.2 (1.0)	4.1 (2.4) **
	Adiponectin (mg/L)	8.5 (4.2)	6.7 (3.7) **	13.1 (5.8)	9.5 (4.6) **
	CRP (mg/L)	3.1 (6.2)	3.5 (5.8) **	3.4 (5.4)	5.1 (7.4) **
	Determinants of cognitive function
	Age (year)	47.4 (14.3)	54.4 (11.9) **	45.3 (13.7)	56.7 (12.6) **
	Education level (%)		**		**
	Low	30.0	37.7	26.2	46.9
	Low-intermediate	37.9	44.4	42.0	40.7
	High-intermediate	24.1	14.5	27.3	11.5
	High	8.0	3.4	4.5	3.7

Table 3 : Association metabolic syndrome and cognitive function in men and women

 3 Values presented as mean (standard deviation). AVLT=Adult Verbal Learning Test, TMT=Trail Making Test, Stroop=Stroop Color and Word Test; † time demanding task: high values represent low test scores. p=p-value derived from regression analysis adjusted for age, education, depressive score, alcohol use and family-relationship. * p < 0.05 Figure 1 Associations of HOMA-IR (A), adiponectin (B), and CRP (C) with cognitive function in men and women Results are presented as mean differences compared to quartile 1. The differences are adjusted for age, education, depressive score, alcohol use and family-relationship. HOMA-IR=homeostatis model assessment insulin-resistance, CRP=C-reactive protein * p-value < 0.05; ** p-value <0.01

	Cognitive test	Men			Women
		n	MetS absent	MetS present	n	MetS absent	MetS present
	General cognitive ability
	Dutch Adult Reading Test	750	62.1 (19.0)	58.3 (18.5)	993	60.4 (18.2)	52.6 (18.1)
	Memory
	AVLT -Working memory	757	4.3 (1.7)	3.9 (1.5)	999	4.6 (1.7)	4.0 (1.6)
	AVLT -Learning	757	32.1 (8.9)	29.2 (8.3)	999	35.5 (8.9)	31.8 (8.9)
	AVLT -Recognition	759	27.7 (2.2)	27.1 (2.5)	1000	28.3 (2.2)	27.7 (2.4)
	AVLT -Delayed recall	759	7.1 (2.9)	6.4 (2.8)	1000	8.3 (2.9)	7.2 (2.9)
	Executive
	Ratio TMT-B / TMT-A †	752	2.7 (1.0)	2.8 (1.0)	996	2.6 (1.0)	2.9 (1.1)
	Ratio Stroop III / Stroop II †	734	1.7 (0.4)	1.8 (0.4)	996	1.7 (0.3)	1.9 (0.6) *
	Verbal fluency	758	63.7 (18.0)	60.0 (17.1)	998	63.0 (17.9)	53.5 (18.7)
	Visuospatial
	Block design	756	30.8 (15.6)	25.7 (13.9)	998	28.7 (15.1)	20.8 (11.2)
	Composite scores
	Memory function	757	0.1 (0.8)	-0.2 (0.8)	999	0.1 (0.8)	-0.3 (0.9)
	Executive function	759	0.1 (0.7)	-0.1 (0.7)	998	0.1 (0.6)	-0.4 (0.8) *
	Global cognitive function	758	0.1 (0.7)	-0.2 (0.6)	998	0.1 (0.7)	-0.3 (0.7)

Acknowledgments

This study is financially supported by the Netherlands Organization for Scientific Research (NWO), the Nutrigenomics Consortium (NGC) and the Centre for Medical Systems Biology (CMSB) in the framework of the Netherlands Genomics Initiative (NGI). We thank Marijke Frolich of the dept. Clinical Chemistry, Leiden University Medical Center, for measuring plasma adiponectin, CRP and insulin levels.


Student Member, IEEE Qin Zhang
email: qin.zhang@lirmm.fr

Member, IEEE Mitsuhiro Hayashibe
email: hayashibe@lirmm.fr

Member, IEEE Philippe Fraisse
email: fraisse@lirmm.fr

Member, IEEE David Guiraud
email: guiraud@lirmm.fr

FES-Induced Torque Prediction with Evoked EMG Sensing for Muscle Fatigue Tracking

Keywords: Evoked Electromyography (eEMG), Torque Prediction, Muscle Fatigue Tracking, Kalman Filter with Forgetting Factor, Functional Electrical Stimulation (FES)

This study investigates a torque estimation method for muscle fatigue tracking, using stimulus evoked electromyography (eEMG) in the context of a functional electrical stimulation (FES) rehabilitation system. Although FES is able to effectively restore motor function in spinal cord injured (SCI) individuals, its application is inevitably restricted by muscle fatigue. In addition, the sensory feedback indicating fatigue is missing in such patients. Therefore, torque estimation is essential to provide feedback or feedforward signal for adaptive FES control. In this work, a fatigue-inducing protocol is conducted on five SCI subjects via transcutaneous electrodes under isometric condition, and eEMG signals are collected by surface electrodes. A myoelectrical mechanical muscle model based on the Hammerstein structure with eEMG as model input is employed to capture muscle contraction dynamics. It is demonstrated that the correlation between eEMG and torque is time-varying during muscle fatigue. Compared to conventional fixed-parameter models, the adaptedparameter model shows better torque prediction performance in fatiguing muscles. It motivates us to use a Kalman filter with forgetting factor for estimating the time-varying parameters and for tracking muscle fatigue. The assessment with experimental data reveals that the identified eEMG-to-torque model properly predicts fatiguing muscle behavior. Furthermore, the performance of the time-varying parameter estimation is efficient, suggesting that real-time tracking is feasible with a Kalman filter and driven by eEMG sensing in the application of FES.

I. INTRODUCTION

A. Background F UNCTIONAL electrical stimulation (FES) is one of the existing solutions to partly restore lost motor function in persons with spinal cord injury (SCI). The electrical stimulus can artificially generate action potential on the axons of the alpha motor neurons to drive muscle contraction in place of the central nervous system (CNS). FES has been used in a wide range of rehabilitation applications, including the FES-aided support of standing, gait, grasping [START_REF] Triolo | Lower Extremity Applications of Functional Neuromuscular Stimulation After Spinal Cord Injury[END_REF] [START_REF] Mangold | Transcutaneous Functional Electrical Stimulation for Grasping in Subjects with Cervical Spinal Cord Injury[END_REF], drop foot correction, tremor compensation and bladder/bowel management.

Although FES has potential advantages for the improvement of functional restoration in SCI subjects, it has not so far gained widespread clinical use because of its limitations, such as rapid muscle fatigue and imprecise force/torque control. In order to obtain the desired movement or trajectory for a certain task, the electrical stimulation (ES) must produce strong, consistent muscle force. It is well known, however, that the fatigue resistance decreases in the paralyzed muscle after SCI [START_REF] Pelletier | Muscle Fatigue Characteristics in Paralyzed Muscle after Spinal Cord Injury[END_REF], so that muscle becomes fatigued more rapidly when artificially stimulated than when excited by the CNS, because of the way that the motor units are recruited: inverse size principle, synchronized activation of motor units and constant order of recruitment [START_REF] Blnder-Macleo | Muscle Fatigue: Clinical Implications for Fatigue Assessment and Neuromuscular Electrical Stimulation[END_REF]. In addition, since SCI individuals have also lost their sensory pathways, they cannot perceive muscle fatigue as stimulation proceeds, which leads to movement failures and sub-optimal FES parameter-tuning in practical FES control. Furthermore, the physiological and mechanical complexity-and nonlinearity-of the neuromuscular system increase the difficulty for FES to precisely control muscle force/torque output and to perform functional movements [START_REF] Chou | The Effects of Stimulation Frequency and Fatigue on the Force-Intensity Relationship for Human Skeletal Muscle[END_REF].

Some research has addressed the attenuation or delay of FES-induced fatigue. Optimal stimulation patterns, such as Nlets [START_REF] Karu | Reducing Muscle Fatigue in FES Applications by Stimulating with N-Let Pulse Trains[END_REF] or catch-like stimulation [START_REF] Shimada | Reduction of Muscle Fatigue by Catchlike-Inducing Intermittent Electrical Stimulation in Rat Skeletal Muscle[END_REF], were found to maximize muscle performance and minimize fatigue. Random modulation of FES parameters was proved to have no effect on muscle fatigue [START_REF] Graham | The Effect of Random Modulation of Functional Electrical Stimulation Parameters on Muscle Fatigue[END_REF]. Low frequency and long pulse duration was reported to be capable of producing less fatigue with constant frequency and intensity [START_REF] Kesar | Effect of Frequency and Pulse Duration on Human Muscle Fatigue During Repetitive Electrical Stimulation[END_REF]. Despite these efforts, no consensus was reached to date due to the complex, multifactorial and task-specific properties of muscle fatigue [START_REF] Russ | Factors in Fatigue During Intermittent Electrical Stimulation of Human Skeletal Muscle[END_REF]. Therefore, muscle fatigue is still a major limiting factor for the widespread application of FES.

Accurate torque/force signal is hence important to produce the desired feedback or feedforward signal for adaptive FES control in the the presence of muscle fatigue. However, the torque/force measurement via external equipments is not convenient for daily use, and generally do not directly measure the muscle output induced by stimulation. The implanted sensor is one possible technique [START_REF] Hoffer | Closed-Loop, Implanted-Sensor, Functional Electrical Stimulation System for Partial Restoration of Motor Functions[END_REF] but is not yet available for practical use. Thus, the inadequacy of torque/force sensors is another problem in rehabilitation application. This motivates the development of methods to estimate muscle torque from biosignals that can be measured or well estimated.

B. Related Work

As described above, on one hand, muscle fatigue may result in torque/force decline and eventually movement failure without any sensation in SCI individuals, on the other hand, convenient and precise sensor is not available for measuring muscle output generated by stimulation. Therefore, estimating the FES-induced torque/force is important for adaptive closedloop FES control, considering the compensation of fatigue. Several techniques have been developed for force/torque prediction in fatiguing muscle. Some mathematical fatigue models were developed from physiological knowledge or experimental analysis. A fatigue function was introduced into a biomechanical model to predict FES-induced shank motion [START_REF] Riener | Biomechanical Model of the Human Knee Evaluated by Neuromuscular Stimulation[END_REF]. A fatigue recovery function, based on metabolic profiles, was introduced into a musculotendon model [START_REF] Mizrahi | Predicted and Measured Muscle Forces After Recoveries of Differing Durations Following Fatigue in Functional Electrical Stimulation[END_REF]. A four-parameter fatigue model, coupled with a mathematic isometric force model, that predicts the fatigue induced by different stimulation patterns under isometric contractions, was reported in [START_REF] Ding | A Predictive Model of Fatigue in Human Skeletal Muscles[END_REF]. However, in the above works, the complex model parameters remain difficult to be identified in application, because of high nonlinearity and dynamic complexities. Moreover, the fatigue models can work only when the stimulation scheme can be predetermined. But in practice, the stimulation pattern is unknown in advance, so the fatigue model cannot work.

Some researchers proposed to use FES-evoked electromyography (eEMG) for torque/force prediction, as ES-induced muscle contraction represents both electrical and mechanical behavior, respectively manifested by eEMG and torque/force. Moreover, the eEMG signal permits noninvasive and reliable measurement of muscle activity. Under continuous stimulation, eEMG was found to be highly correlated with FES-induced torque in different muscle conditions, suggesting its use as a fatigue indicator [START_REF] Chesler | Surface EMG as a Fatigue Indicator During FES-induced Isometric Muscle Contractions[END_REF] [START_REF] Erfanian | Excitation-Contraction Fatigue During Sustained Electrical Stimulation of Paralyzed Muscle[END_REF]. An exponential function was used to express the electrical and mechanical relationship of the stimulated muscle in SCI patients [START_REF] Mizrahi | EMG as an indicator of fatigue in isometrically FES-activated paralyzed muscles[END_REF] [START_REF] Tepavac | Detection and Prediction of FES-Induced Fatigue[END_REF]. A predictive model based on eEMG instead of stimulation was developed, allowing the use of eEMG as a synthetic torque sensor [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically Stimulated Muscle[END_REF]. One possible explanation is that eEMG can capture most time-variations due to the effects of fatigue, which cannot be captured by a stimulation-to-torque model. Combining timedomain and frequency-domain variables was suggested for better prediction of force [START_REF] Tepavac | Detection and Prediction of FES-Induced Fatigue[END_REF]. In these works, a fixed relationship between eEMG and muscle force/torque during sustained stimulation was assumed. Nevertheless, this correlation was presented to be time-variant during different muscle-fatigue levels [START_REF] Mizrahi | Myoelectric and Force Characteristics in Transcutanteous Isometric FES[END_REF] and during recovery process [START_REF] Yu | Mechanical and Electromyographic Response to Stimulated Contractions in Paralyzed Tibialis Anterior Post Fatiguing Stimulations[END_REF]. Such evidence indicates that the relationship between eEMG and muscle force/torque is not constant but time-varying, and also implies the limitations of previous force/torque prediction strategies based on fixed-parameter models. [START_REF] Yu | Mechanical and Electromyographic Response to Stimulated Contractions in Paralyzed Tibialis Anterior Post Fatiguing Stimulations[END_REF] suggested an adaptive tuning of the parameters of eEMG for predicting the stimulated force. However, they did not propose a feasible and effective method to do this.

The present work aims to develop an online estimation method of FES-induced torque based on eEMG in the presence of muscle fatigue. This paper is organized as follows. The relevant experiment is introduced in Section II. The muscle modeling is presented in Section III. The proposed estimation method, Kalman filter with forgetting factor, is presented in Section IV. In Section V, we show the identification and prediction results with the ordinary least squares method and the Kalman filter technique. Discussions and conclusions can be found in Sections VI and VII, respectively.

II. EXPERIMENTAL METHODS

A. Experimental set-up

The experiments were conducted on five SCI subjects (see TABLE I) in the PROPARA rehabilitation center, Montpellier, France. All subjects were classified as the American Spinal Injury Association (ASIA) A, where no motor or sensory function is preserved. The experimental set-up is depicted in Fig. 1. This study was approved by the ethical committee of France and all subjects signed informed consent forms. The subjects were seated on the chair with the ankle at 90 o , while the joint center was aligned with the axis of a calibrated dynamometer (Biodex 3, Shirley corp., NY, USA). The shank was adjusted to be horizontal to the ground with the knee joint at 40 o . The foot was strapped to the pedal to transmit ankle torque to the dynamometer and to allow the optimal recording of isometric ankle torque. Electrical current pulses were delivered to the right triceps surae muscle group via surface electrodes (10cm×3cm) to induce muscle contractions and to plantarflex the ankle joint. One electrode was placed 5cm beneath the popliteal cavity and the other beneath the insertion point of the medial and lateral gastrocnemius on the Achilles tendon. The muscle group was stimulated with amplitude modulation at a constant frequency (30Hz) and constant pulse-width (450μs), under isometric conditions, by a portable stimulator (Cefar physio 4, Cefar Medical, Lund, Sweden).

EEMG activity of soleus in triceps surae muscle group was recorded, amplified (gain 1000) and sampled at 4KHz by an acquisition system (Biopac MP100, Biopac Systems Inc., Santa Barbara, CA, USA). The bipolar AgCl EMG electrodes were positioned over the muscle belly in the direction of muscle fiber with 20mm interelectrode spacing. The reference electrode was placed on the patella of another leg. The skin under the electrodes was shaved to minimize the impedance. Isometric ankle plantar-flexion torque was measured using the dynamometer (Biodex 3), sampled at 2KHz, and interfaced with the acquisition system (Biopac MP100).

B. Experimental Protocol

For each subject, the experiment consists of three test sessions: a fatigue-inducing test, a fatigue-recovery test and a random test, as shown in Fig. 2. The maximum stimulation amplitude was found for each subject at the beginning, by gradually increasing stimulation amplitude until torque saturation. The fatigue-inducing test includes several sequences (named as fatigue1-fatigue5). Each sequence contains five trapezoidal trains with each trapezoidal train consisting of 4s stimulation (1s ramp-up, 2s plateau and 1s ramp-down) and 2s rest. The stimulation amplitude during plateau is chosen at 50% of the maximum stimulation amplitude with constant stimulation frequency 30Hz and constant stimulation pulse width 450μs. In order to induce muscle fatigue, three such stimulation sequences were applied to subjects S1 and S2, four sequences to S3, and five sequences to S4 and S5. In fatigue-recovery test, one trapezoidal train (indicated by pre, postA, ..., postE and recovery1, ... recovery3) at maximum stimulation amplitude was delivered to the muscle just before and after each fatigue-inducing sequence. After the stimulation train postE, the same stimulation train was applied every 5 minutes, up to 15 minutes. At the end, a sequence including several trapezoidal trains was applied. In this case, the stimulation amplitude in each train was increased from zero to a randomly determined value below the maximum value, and then symmetrically decreased, during two minutes in total.

C. Torque and eEMG Data Processing

During electrically elicited contractions, the detected surface eEMG signal is contaminated by stimulation artifacts due to the stimulation current which is recorded by the eEMG amplifier. For the preparation of model identification and torque prediction based on eEMG, the measured eEMG and torque were processed in the following steps. An example of the processed results is illustrated in Fig. 3.

1) The blanking window method [START_REF] Frigo | EMG Signals Detection and Processing for On-Line Control of Functional Electrical Stimulation[END_REF] is used to remove stimulation artifacts from the raw eEMG signal and extract muscle response (Mwave).

2) A lowpass filter is applied to measured ankle torque (6thorder, cutoff frequency 100Hz) and measured eEMG (6th-order, cutoff frequency 300Hz).

3) The filtered eEMG signal is divided into epochs with each epoch containing one Mwave, and the mean absolute value (MAV) of eEMG is calculated every five epochs. The average torque is calculated within the same time window. 4) The MAV and average torque are normalized with respect to their maximum values. The normalized MAV and normalized torque are prepared to be the system input and output for model identification. An example of processed eEMG and torque. The raw eEMG signal and Mwave were zoomed in (0.5s ∼ 1.1s) to show the details during stimulation increase. The raw eEMG was contaminated by stimulation artifacts. Blanking window was applied to remove artifacts so that Mwave was effectively extracted.

III. CONTRACTION DYNAMICS MODEL OF ELECTRICALLY STIMULATED MUSCLE

A. Model Structure

The discrete-time Hammerstein structure was used to model the muscle contraction dynamics as shown in Fig. 4. This model consists of a memoryless nonlinear part followed by a dynamic linear part. It is popularly used to represent highly nonlinear systems, and particularly useful for modeling biomechanical systems, such as stretch reflex EMG signal [START_REF] Dempsey | Identification of Hammerstein Mod-els With Cubic Spline Nonlinearities[END_REF] and electrically-stimulated muscles, relating stimulation to muscle force under isometric conditions [START_REF] Chia | Recursive Parameter Identification of Constrained Systems: An Applicaiton to Electrically Stimulated Muscle[END_REF]. It has been shown to extend to dynamic conditions [START_REF] Farahat | A Method Identification of Electrically Stimulated Muscle[END_REF], which is essential for developing stable adaptive controllers for applications in FES. In this study, the MAV of eEMG and FES-induced ankle torque are system input u(t) and output y(t), respectively. The memoryless nonlinear function maps the system input u(t), to the intermediate variable h(t), which represents the activation level of the stimulated muscle. It is traditionally modeled by an nth-order polynomial of u(t) [START_REF] Dempsey | Identification of Hammerstein Mod-els With Cubic Spline Nonlinearities[END_REF] as follows:

h(t) = n i=1 γ i u i (t) (1)
where γ i is a model coefficient.

The linear time-variant system is described by an autoregressive model with exogenous input (ARX), which has been shown experimentally to yield good prediction of output torque/force in isometric situation [START_REF] Bernotas | A Discrete-Time Model of Electrically Stimulated Muscle[END_REF]. It can be described as

A(z)y(t) = B(z)h(t) + e(t) (2)
with transfer function G(z) = B(z)/A(z) and

A(z) = 1 + a 1 z -1 + a 2 z -2 + • • • + a l z -l B(z) = b 1 z -1 + b 2 z -2 + • • • + b m z -m (3)
where z -1 is the backward shift operator which makes z -1 y(t) = y(t -1). In (2), y(t) is the system output at time t, e(t) is zero mean and Gaussian white noise affecting the system. The h(t) is the output of the nonlinear element, and the input of the linear element as seen in Fig. 4. Substituting (1) and ( 3) into (2) and expanding, the output of a polynomial Hammerstein model (PHM) at a given time t can be parameterized as:

y(t, θ) = l i=1 a i y(t -i) + m i=1 n j=1 μ ij (u(t -i)) j ( 4
)
where

μ ij = b i γ j , θ = [a 1 , • • • , a l , μ 11 , • • • , μ mn ] T is a
parameter vector containing the model coefficients. The size of θ depends on model complexity. Therefore, the selection of model order (l, m, n) is a key step in the estimation of the unknown parameters in θ. We chose to model the recruitment curve of the muscle as a 3rd-order polynomial of instantaneous MAV of eEMG, (n = 3), as in [START_REF] Chia | Recursive Parameter Identification of Constrained Systems: An Applicaiton to Electrically Stimulated Muscle[END_REF]. Linear model order determination was determined by comparing the Rissanen's minimum description length (MDL) [START_REF] Rissanen | Modeling by Shortest Data Description[END_REF][28] obtained for different model orders, since the MDL principle provides a criterion for tradeoff between the simplicity of the model and the model's applicability to the data. Model-order parameters ranging from 2 ≤ l ≤ 6, and from 2 ≤ m ≤ 6 were considered. Finally, model order (l = 3, m = 4) was chosen with relatively less MDL value and a simpler model as shown in Fig. 5. The stimulated muscle model has l+m×n unknown parameters in all. The elements of θ as well as the eEMG-totorque features are time-varying due to the effects of fatigue and the associated biochemistry. At a given time t, the model estimates are predicted using (4) by assuming the system is stationary, or slowly time-varying, during the prediction horizon. The computation of past torque y(t -i) in ( 4) has two possible versions: (i) past measured torque y m (t -i) and (ii) past predicted torque y p (t -i). For model identification, the former approach is preferred, as measured torque is favorable for model fitting. For torque prediction, both can be employed, while from a practical point of view, the former is inappropriate, since measured torque is not available. Accordingly, the corresponding predicted output y(t) based on identified model has two versions, y m (t) and y p (t), which can be computed as a function of past measured MAV of eEMG, past measured torque or past predicted torque in this way:

y m (t) =f unc[u(t -1), u(t -2), ...u(t -m), y m (t -1), y m (t -2), ...y m (t -l)] ( 5
)
y p (t) =f unc[u(t -1), u(t -2), ...u(t -m), y p (t -1), y p (t -2), ... y p (t -l)] (6)
When torque measurement is not available, the second approach as shown in ( 6) makes it possible to use eEMG as a synthetic torque sensor. In this case, we can initialize the predicted torque at zero when no stimulation is delivered to the muscle.

B. State-Space Model Representation

State-space form is basically required for the implementation of Kalman filter. Considering a PHM (l, m, n) as in (4), its state-space form can be compactly written as:

1) process equation

x k = f (x k-1 , u k-1 ) = Ax k-1 + n i=1 B i [u k-1 ] i (7)
2) measurement equation

y k = Cx k (8)
where the subscript k indicates the current time step. In [START_REF] Shimada | Reduction of Muscle Fatigue by Catchlike-Inducing Intermittent Electrical Stimulation in Rat Skeletal Muscle[END_REF], the current state vector

x k = x 1,k , x 2,k , • • • , x q,k T , q = max{l,
A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a 1 1 0 • • • 0 0 a 2 0 1 • • • 0 0 . . . . . . . . . . . . . . . . . . a q-1 0 0 • • • 0 1 a q 0 0 • • • 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , B i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ μ 1,i μ 2,i . . . μ q-1,i μ q,i ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
.

The y k in ( 8) is the measurement of system output. C ∈ R 1×q relates the current state x k to the current measurement y k with the following expression:

C = 1 0 • • • 0 0 .
Note that, in practice matrices A, B i and C might change with each time step or measurement, but here we assume they are constant for simplifying the expression.

IV. IDENTIFICATION USING KALMAN FILTER WITH FORGETTING FACTOR

Recursive estimates are important when parameter identification is needed in real-time, such as when the information is used in adaptive control or real-time diagnosis. The Kalman filter is an efficient recursive filter that estimates the internal states and parameters of a discrete-time system from a series of noisy measurements. Online estimation of the states in x and of the model parameters in θ is investigated simultaneously in this work. Model parameters to be identified are the coefficients relating to both the past measured torque and the past measured MAV of eEMG in [START_REF] Blnder-Macleo | Muscle Fatigue: Clinical Implications for Fatigue Assessment and Neuromuscular Electrical Stimulation[END_REF]. A Kalman filter for state and parameter estimation is performed by regarding the unknown model parameters as elements of the state vector. In this way, the basic Kalman filter algorithm does not need to be modified, except that the state vector x will be augmented with the unknown parameters in θ. That is, the meta-state vector w k has the expression w k = [x k ; θ k ], with subscription k indicating the time step and used in the subsequent explanation. The parameters in θ are assumed to be locally time-invariant or slowly varying compared to the process. Accordingly, the augmented system is described by

w k = F(w k-1 , u k-1 ) (9) y k = Hw k ( 10
)
where

F(w k-1 , u k-1 ) = f (x k-1 , u k-1 ) θ k-1 , H = [C 0 1×[l+m×n] ].
The recursive estimation of the state-space model with Kalman filter consists of two phases, prediction and correction. The main equations are given by: 1) prediction phase:

ŵ- k = F( ŵk-1 , u k-1 ) (11)
P - k = D k P k-1 D T k + Q k-1 (12)
2) correction phase:

K k = P - k H T k (H k P - k H T k + R k ) -1 (13) ŵk = ŵ- k + K k (y k -H k ŵ- k ) (14)
P k = (I -K k H k )P - k (15)
In the prediction phase, at step k, the a priori estimate of the state ŵk , is given by the a posteriori state at the previous step, ŵk-1 , as shown in [START_REF] Hoffer | Closed-Loop, Implanted-Sensor, Functional Electrical Stimulation System for Partial Restoration of Motor Functions[END_REF]. The estimate error covariance P k is propagated according to [START_REF] Riener | Biomechanical Model of the Human Knee Evaluated by Neuromuscular Stimulation[END_REF], wherein Q k is a diagonal matrix containing the process noise covariance, and D k is the process Jacobian with respect to the variables involved, with each element D [i,j] computed by:

D [i,j] = ∂F [i] ∂w [j] ( ŵk-1 , u k-1 )
In the correction phase, K k in ( 13) is called Kalman filter gain, wherein R k is a scalar measurement noise covariance. The updated state is computed in [START_REF] Ding | A Predictive Model of Fatigue in Human Skeletal Muscles[END_REF], and the updated estimate error covariance is given by [START_REF] Chesler | Surface EMG as a Fatigue Indicator During FES-induced Isometric Muscle Contractions[END_REF].

Although the Kalman filter is an effective way of estimating the state and parameters of a discrete-time controlled process, its performance in estimating the time-varying parameters is degraded by the fact that it refers to the entire history of past measurements [START_REF] Mack | Speech Parameter Estimation by Timeweighted-Error Kalman Filter[END_REF]. This is particularly troublesome since the activity of stimulated muscles may vary, with prolonged or repetitive stimulation leading to different muscle states. In order to track the time-varying muscle condition, a forgetting factor λ is deliberately introduced as proposed in [START_REF] Feng | Unscented Kalman Filter for Time Varying Spectral Analysis of Earthquake Ground Motions[END_REF]. Consequently, equation ( 12) and ( 13) can be rewritten as

P - k = D k P k-1 D T k /λ (16) K k = P - k H T k (H k P - k H T k + λ) -1 (17)
Choosing forgetting factor λ ∈ [0, 1] depends on how much we hope the filter to forget the past measurements. The forgetting factor is closer to 1, the filter will forget fewer past measurements. A tradeoff between the smoothness of tracking and lag in detecting the changes in model parameters should be considered when forgetting factor is introduced to a Kalman filter. Usually λ ∈ [0.9, 1] is suitable for most application.

V. RESULTS

A. Muscle Fatigue and Recovery Process

The data from fatigue-recovery tests were used to observe the relationship between torque and MAV of eEMG during intermittent stimulation and recovery. The data during stimulation plateau was used for analysis. The data were firstly treated as the processing step (1) (2) described in section II. C. Then the average torque, MAV and standard deviations were computed simultaneously. The results are depicted in Fig. 6. The torque of ankle plantar-flexion gradually declined after each fatigue-inducing sequence, until around 90% of initial torque in all subjects. After 5 minutes rest, the average torque recovered less than 5% and then remained at the same level in all subjects. Except in subject S1, the torque recovered after 10 minutes rest. As a whole, torque transition in fatigue generally showed a similar tendency in all subjects. The MAV of eEMG represented different transitions among these subjects as depicted in Fig. 6 (b). Although the same tendency can be found in subjects S3 and S5, different tendencies are found among the five subjects. The results in S3, S5 show a simple decline of MAV due to fatigue, S4 shows potentiation phenomenon, while S1, S2 represent somewhat combined characteristics of potentiation and fatigue. Here, we do not focus on the understanding of the different characteristics. However, we were able to confirm that the eEMG-to-torque relationship was not constant as most prior works assumed, but gradually varying, in this intermittent fatigue protocol.

B. Model Identification and Validation Using the Ordinary Least Squares (OLS) Method

The torque was predicted based on past measured torque and past predicted torque as described in section III. A, using random test data. OLS method was performed for model identification. Fig. 7 illustrates the prediction result in subject S5. The data during 97s of random stimulation were separated into two parts at 68s. The model parameters in (4) were identified using the data before 68s. Then torque prediction was calculated using two different methods, as expressed in ( 5) and ( 6) respectively, with the data before and after 68s. It means that the torque prediction after 68s is only driven by eEMG information. The corresponding prediction errors, root-mean-square (RMS) error, are also shown in this figure. From the result of the random stimulation test, we can confirm the feasibility of model identification and its prediction performance. For the next step, we need to investigate this in different muscle fatigue states to verify the model predictability. The torque prediction was conducted using the OLS method with the fatigue-inducing data in 5 subjects. Two different models, the fixed-parameter model and the adapted-parameter model, were tested. The fixed-parameter model can be explained as follows: the model parameters were identified with the data of the first sequence (fatigue1), and were then crossvalidated for all the remaining sequences (fatigue2-fatigue5). The prediction errors are shown in TABLE II. We found that with the fixed-parameter model, the prediction error became higher and higher as the muscle was more and more fatigued. We then supposed that if the model was identified again with the data in the latest sequence, the prediction could be improved. Therefore, the second method, the adaptedparameter model, was proposed. In this approach, the latest data were used to re-identify the model. That is, the torque prediction of fatigue3 was based on the model parameters obtained from fatigue2, and so on. The prediction errors with the adapted-parameter model in all subjects are compared in TABLE II. For example, for fatigue3 of subject S3, the RMS error with the adapted model is 0.0345, as compared to 0.0639 with the fixed model. The average prediction of the adapted model was superior by 16.7-50.8% compared to the fixed model in all subjects.

In Fig. 8, the predictions with the two methods (fixed and adapted model) and in different fatigue conditions are shown in subject S3. Fatigue3 was not plotted in this figure, as there was only a small difference between fatigue3 and fatigue4. Obviously, FES-generated torque declined with the same stimulation as a result of muscle fatigue. The fixed model could be still used for torque prediction. However, in comparison with fatigue1 and fatigue2, the prediction of fatigue4 became less precise based on the fixed model, whereas the adapted model contributed to improve torque prediction accuracy, as the

C. Validation of Time-Varying Model Estimation with Kalman Filter

In this part, the time-varying parameter estimation performed by the Kalman filter was evaluated in simulations and with experimental data. The outlines of the muscle model and Kalman filter were described previously. For a PHM (l, m, n) model, we need to estimate the r = max(l, m) + (l + m × n) dimensional meta-state. The max(l, m) parameters relate to the internal states, the rest relate to the past torque and the past MAV of eEMG. The elements of the meta-state vector were initialized as ŵi (0) = 0, i = 1, 2 • • • , r. The initial output estimate ŷ0 was set at zero. The estimate error covariance was initialized as P 0 = I, where I is an identity matrix.

1) Time-Variant Parameter Tracking in Simulation:

In simulation, invariant parameter tracking was evaluated first with the Kalman filter to investigate the stability of the muscle model. Secondly, in order to investigate the filter's robustness to the time-varying fatigue phenomenon, we slowly changed the model parameters at different instants to imitate changes in muscle condition. The advantage of simulation is that the true parameters are known to be compared with the estimated ones. The simulation model order was chosen as l = 2, m = 2, n = 1 to reduce model complexity, as it is difficult to know how the model output changes when too many parameters change. Thus, four parameters, a 1 , a 2 , μ 11 , μ 21 , were estimated via the Kalman filter algorithm in simulation. At the beginning, all the parameters were kept constant. After 33s or 50s, they were changed linearly (a 1 and μ 21 ) or in steps (a 2 and μ 11 ). The pseudorandom binary sequence (PRBS), which is commonly used in muscle identification [START_REF] Farahat | A Method Identification of Electrically Stimulated Muscle[END_REF], was chosen as model input. Model input, output and the a posteriori estimate of the output are shown in Fig. 9.

Model input

Model input signal Fig. 9. Data set for model identification. Data were generated by simulation using a PHM (2, 2, 1). A Kalman filter with forgetting factor 0.97 was used to identify the model.

The corresponding parameter estimates of the model are depicted in Fig. 10. The solid lines indicate the true parameters, while the dotted lines indicate parameter estimates. All the parameters converge steadily after 5s when the parameters are static. After 33s or 50s, the model parameters gradually vary, and the estimates track the changes well, which implies that the estimation method is suitable for time-variant parameter tracking with a PHM model.

2) Fatigue Tracking Based on Experimental Data: In the model identification described in section V. B, the adaptedparameter model was proved to be able to improve torque prediction, but identification was performed without automatic tracking function. In this section, the identification and validation of time-varying parameters are considered and performed automatically using the Kalman filter. The model order is chosen at (3, 4, 3) as described above. The data from successive series of fatigue-inducing tests were concatenated The time-varying property of the poles may also interpret the resulting time-varying model parameters. Moreover, the locations and movements of the z-plane poles in all subjects present similar characteristics, suggesting that it is possible to assess muscle fatigue dynamics from such information. In general, the damping ratio increases when the muscle is highly fatigued. This matches our intuition concerning the effect of muscle fatigue. prediction performance in time, different prediction horizons, 6s, 18s and 30s, were tested. The idea is illustrated in Fig. 12, which features the estimated model at time instant t0. The torque predictions were computed using (4) and only driven by the eEMG, considering a static system within the prediction horizon. At every torque sample for a given prediction horizon, the RMS error and peak error were evaluated. The prediction errors in all subjects are quantified in Table III. Eighteen seconds was considered to be an appropriate prediction horizon, as it provided a tradeoff between sufficient interval for measurement update in Kalman filter and satisfactory prediction performance. The prediction errors for the 18s prediction horizon in subject S3 are plotted in Fig. 13. The solid blue line indicates RMS error, while the dotted black line indicates peak error. Predictive Performance in Sensing Failure: In order to evaluate tracking robustness in the event of interruption in recording or failure, torque prediction was integrated into the fatigue tracking task. Any event such as a missing torque measurement, was assumed to result in tracking failure. Fig. 14 reveals fatigue tracking performance in subject S3 when muscle torque declined to 76% of the maximum torque during the 100s-stimulation. Assuming measurement was not available from time instant t0 to t1, the estimated model at t0 was used to predict torque until instant t1. When prediction was executed, the model was only driven by the eEMG, while the online model estimation was suspended. To evaluate prediction performance in different muscle fatigue states, this process was repeated until the end of stimulation, where the measurement update in Kalman filter was switched off for 18s.

The prediction tracks the measured torque well, suggesting that if torque measurement is unavailable or unreliable, the predicted torque based on eEMG can be used to bridge such gaps for the prolonged application of FES. Some of the estimated parameters are also depicted in Fig. 14. The middle plot shows the parameters of the past torque, while the parameters of the 3rd power of the past MAV are shown in the lower plot. Although the parameters gradually vary, representing varying myoelectrical mechanical correlation during muscle fatigue, the proposed muscle model and Kalman filter algorithm with forgetting factor provide good prediction performance, as shown in the upper plot.

VI. DISCUSSION

Stimulus-evoked EMG has been previously proposed for the prediction of the mechanical behavior of muscles in FES rehabilitation systems [START_REF] Tepavac | Detection and Prediction of FES-Induced Fatigue[END_REF] [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically Stimulated Muscle[END_REF]. These prediction methods were investigated on the hypothesis of a fixed eEMGto-torque model during ES. Some other researchers found that this relationship was time-variant under different muscle conditions [START_REF] Mizrahi | Myoelectric and Force Characteristics in Transcutanteous Isometric FES[END_REF] [START_REF] Yu | Mechanical and Electromyographic Response to Stimulated Contractions in Paralyzed Tibialis Anterior Post Fatiguing Stimulations[END_REF]. However, they have not proposed a feasible, effective method for torque prediction. The timevarying property of myoelectrical and mechanical muscle behavior is shown in Fig. 6. It implies the limitations of the torque prediction method based on the fixed eEMG-to-torque model and increases the difficulties of torque estimation. The present work proposes the use of the Hammerstein structure to represent muscle contraction dynamics with MAV of eEMG as input, where model estimation is performed by a Kalman filter for fatigue tracking.

The prediction result presented in Fig. 12 reveals good performance of the predictive filter, successively acquiring the mechanical behavior of the muscle under FES. In addition, the error values given in Table III indicate that neither RMS error nor peak error increases significantly when the prediction horizon is extended. It indicates that the proposed eEMGto-torque model properly fits muscle behavior. If prediction quality diminishes when an expanded prediction horizon is selected, this is probably due to the variations of muscle fatigue levels during the horizon. This can occur when prolonged or repetitive stimulation is delivered to the muscle. In this application, 18s prediction horizon is considered as the optimal choice.

The muscle contraction dynamics model has been described in a previous work [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically Stimulated Muscle[END_REF], in continuous stimulation, with the recursive least squares method used for identification. However, there was no significant consideration of different fatigue conditions in intermittent stimulation, as illustrated in Fig. 14. In this study, the proposed estimation method was validated for ankle torque prediction and fatigue tracking using eEMG under isometric condition, it would be promising to verify it in dynamic conditions by introducing a torque-joint angle function to the proposed muscle model as in [START_REF] Farahat | A Method Identification of Electrically Stimulated Muscle[END_REF].

In this work, we did not consider the effect of day-to-day changes, but as the online identification is able to identify parameters for different subjects, which implies that this method itself can work even if there are day-to-day changes if we do not care about how the parameters change depending on the experimental set-up, for example, electrode position. Of course, it is important to investigate the effect of different experimental set-ups. In this case, we may introduce a tuning function to offset the day-to-day variances.

VII. CONCLUSION

The objective of this work is to develop a torque estimation method which can be used to track FES-induced muscle fatigue. In this work, we confirmed that muscle contraction model parameters were time-varying during intermittent stimulation in five SCI subjects. An estimation method of FESinduced torque based on eEMG signals was evaluated in its prediction performance. A time-varying eEMG-to-torque model was employed to represent the myoelectrical and mechanical behavior of stimulated muscles, where the model parameters were estimated by a Kalman filter with forgetting factor. The results of the proposed method for the fatigue tracking task represent a feasible and effective torque prediction performance in isometric condition in all subjects. In terms of improvements in fatigue tracking, when the measurement of torque sensor suffers from recording interruption, the proposed method can bridge these problems and provide sufficiently accurate fatigue tracking only on eEMG measurement. Therefore, we conclude that the proposed estimation method can contribute to precise torque prediction in presence of muscle fatigue. The predicted torque can be further used for adaptive closed-loop FES control, considering the compensation of muscle fatigue. Future work will be extended to adaptive FES control for muscle fatigue compensation in SCI patients. We will also endeavor to further improve the identification by means of multiple eEMG measurements and more selective stimulation. By introducing a joint dynamics model and joint angle sensors along with the proposed method, we will work on FES close-loop control based on eEMG sensing without the use of torque sensors.

Fig. 1 .

 1 Fig. 1.Experimental set-up for electrical stimulation and ankle torque measurement.

Fig. 2 .

 2 Fig. 2. Schematic representation of the experimental sessions.

  Fig.3.An example of processed eEMG and torque. The raw eEMG signal and Mwave were zoomed in (0.5s ∼ 1.1s) to show the details during stimulation increase. The raw eEMG was contaminated by stimulation artifacts. Blanking window was applied to remove artifacts so that Mwave was effectively extracted.

Fig. 4 .

 4 Fig. 4. Muscle model structure. The present work focuses on the contraction dynamics model contained in the rectangle within the dashed line. The memoryless nonlinearity is modeled by an nth-order polynomial function. The linear dynamics is chosen as an ARX model.

  Fig. 5. Criterion of fit for all datasets of fatigue-inducing test in 5 subjects.'•'---S1, '*'---S2, 'o'---S3, ' '---S4, '+'---S5. The vertical line indicates the selected model order[START_REF] Pelletier | Muscle Fatigue Characteristics in Paralyzed Muscle after Spinal Cord Injury[END_REF][START_REF] Blnder-Macleo | Muscle Fatigue: Clinical Implications for Fatigue Assessment and Neuromuscular Electrical Stimulation[END_REF].

Fig. 6 .

 6 Fig. 6. (a) Average torque and (b) MAV of eEMG in all subjects during the fatigue-recovery process. The results were represented as mean±standard deviation.

Fig. 7 .

 7 Fig. 7. The measured and predicted torque obtained by eEMG-to-torque model in random test in subject S5. The vertical dotted line is shown to separate the data into two parts. The model was identified by OLS method using the data in the left part. (1) The torque was predicted based on the past measured MAV of eEMG and past measured torque. (2) The torque was predicted based on the past measured MAV of eEMG and past predicted torque, where the prediction was calculated only by eEMG.

Fig. 8 .

 8 Fig.8. The measured and predicted torque obtained by eEMG-to-torque model in fatigue-inducing protocol in subject S3. The solid blue lines indicate the measured torques. The dotted red lines represent the corresponding predicted torques based on fixed model which was identified using the data of fatigue1. The dashdotted black lines represent the torque prediction of fatigue4 based on adapted model which was identified using the data of fatigue3. All the torques were normalized by the maximum measured torque in fatigue1.

Fig. 10 .

 10 Fig. 10. Convergence and tracking of both static and time-variant parameters. Solid lines indicate true parameters. Dotted lines indicate parameter estimates. All the initial estimates were set at zero.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. The changes in the contraction dynamics of muscle behavior due to fatigue during intermittent stimulation in subjects S1 (left) and S3 (right). The arrows denote the direction of "movement" of the z-plane pole. The plot color was changed every 16s to show the time transition. Torque Prediction Performance: To investigate the torque

Fig. 13 .

 13 Fig. 13. RMS and peak prediction error of 18s prediction horizon in S3.

Fig. 14 .

 14 Fig.14. Fatigue tracking based on eEMG-to-torque model and Kalman filter in subject S3. The model estimated at t0 was used to predict the torque from t0 to t1, while the online estimation was switched off for 18s (with green background color in upper plot). This process was repeated until the end of stimulation. Upper: The solid blue line indicates the measured torque, while the dotted red line indicates the predicted torque. Middle: The estimated model parameters of the past torque. Lower: The estimated model parameters of the 3rd power of the past MAV.

  m}, and u k-1 is the previous model input. A ∈ R q×q relates the previous state x k-1 to the current state x k . B i ∈ R q×1 relates the previous model input [u k-1 ] i to the current state x k . They can be represented as following: