Abdelghani Chahmi
email: abdelghani.chahmi@gmail.com

Mokhtar Bendjebbar
email: bendjebb_dz@yahoo.fr

Bertrand Raison
email: bertrand.raison@g2elab.grenoble-inp.fr

Détection De Défaillances Sur Un Entraînement Electrique-Approche Signal

Keywords: MOT-CLES-Moteur asynchrone, Modélisation, Défaut rotorique, Défaut statorique, Analyse fréquentielle. I

RESUME-Cet article traite de la détection de défauts dans les entrainements électriques et en particulier ceux à base de machine asynchrone. D'une part, on démontre la caractérisation des processus de défaillance et l'obtention de modèles de ces dysfonctionnements. D'autre part le développement de diagnostic de la machine. Les études de dysfonctionnement dans les entrainements ne pouvant faire appel qu'à peu de données expérimentales, il a été indispensable de nous munir d'outils de simulation qui ont servi à caractériser les comportements défaillants. La détection de défauts fait appel aux techniques de traitement du signal dans les phases de fonctionnement connues (vitesse fixe). Finalement, nous nous intéressons a développé l'algorithme de modélisation de la machine asynchrone saine et défaillant par la signature fréquentielle (apparition des raies, amplitudes de la composante continue,…..).

INTRODUCTION

Les entrainements électriques à base de machine asynchrone sont largement utilisés dans les applications industrielles en raison de leur faible cout, de leurs performances et de leur robustesse. Les systèmes industriels devenant de plus en plus complexes et sophistiqués, il est nécessaire de pouvoir garantir la sûreté de leur fonctionnement. On entend par sûreté de fonctionnement, la capacité qu'aura un système à effectuer les tâches pour lesquelles il a été conçu. Cette capacité peut être entravée par les défauts affectant le fonctionnement du système, entraînant ainsi des conséquences inacceptables sur la sécurité (des hommes ou du matériel), l'environnement et la qualité de production ou de service. Il est essentiel de mettre en place des fonctions de surveillance permettant de garantir le fonctionnement sûr des systèmes. On entend par fonctions de surveillance l'ensemble des actions permettant de détecter, localiser et diagnostiquer tout phénomène anormal survenant pendant le fonctionnement. Ces fonctions de surveillance sont élaborées à partir des techniques dites DLD (Détection et localisation de Défauts) ou FDI (Fault Détection and Isolation).Les références [START_REF] Diao | Sensor Fault Diagnosis for Improving the Availability of Electrical Drives[END_REF]- [START_REF] Faucher | Fault indicator in induction machine using stator current spectral analysis" Power Electronics, machine and drives[END_REF] passent en revue un grand nombre de cas de défaillances, ainsi que des techniques de détection associées. La raison principale des défaillances de machine reste les défauts de roulement à bille [START_REF] Ceban | éthode globale de diagnostic des machines électriques" hèse de doctorat[END_REF]- [START_REF]Surveillance et diagnostic des défauts rotoriques et mécaniques de la machine asynchrone avec alimentation équilibrée ou déséquilibrée" hèse de doctorat[END_REF] qui ont de nombreuses causes telles que la contamination du lubrifiant, une charge excessive ou encore des causes électriques. Dans cet article, un récapitulatif des modèles des défauts, choix de la variable à étudier, de l'influence des boucles de contrôle et les signatures spectrales du courant statorique est donné dans le paragraphe II. Le paragraphe III rappelle quant à lui les effets mécaniques des défauts de roulements et la formulation d'oscillation du couple de charge.

II. LES DÉFAILLANCES ÉLECTRIQUES

A. Les modeles des défauts dans la machine

Il faut noter que cette modélisation est critiquable: les défauts de type court-circuit sont représenté dans le modèle comme des suppressions de spires c'est-à-dire des ouvertures partielles de phase ; cela ne correspond pas tout à fait à la réalité puisqu'il faudrait aussi tenir compte des spires en courtcircuit. Ce modèle constitue un progrès par rapport aux modèles équilibrés ou /et ne considérant qu'une variation de résistance couramment utilisés pour valider les algorithmes. Sous les conditions rappelées au paragraphe précèdent, les équations des circuits électriques font intervenir des inductances propres et mutuelles permettant de définir les flux en fonction des courants [START_REF] Raison | Détection et localisation de défaillance sur un entraiinement électrique[END_REF]. Au niveau du stator, les tensions vérifient :

[ ] [ ] [ ] [ ] ( 1 )
Ou la matrice des résistances s'écrit d'après l'hypothèse de symétrie :

[ ] = [ ] ( 2 )
Au niveau du rotor, sous l'hypothèse que le rotor est assimilable à un rotor avec enroulement triphasés, les tensions s'écrivent :

[ ]=[ ] [ ] d dt [ ] ( 3 )
La matrice R R posséde la méme forme que R S Les flux s'expriment alors grâce aux matrices inductances et aux courants

[ ]=[ ] [ ]+[ ] [ ] ( 4 )
[ ]=[ ] [ ]+[ ] [ ] ( 5 )
Afin de pouvoir rendre compte des interactions entre les phases , nous avons exprimé chaque matrice résistance et inductance en fonction du nombre de spires des enroulements de la machine. Ce modèle permet en particulier de crées des défauts asymétriques que ce soit au stator ou au rotor par la variation du nombre de spires de l'enroulement de la phase en défaut. Le défaut peut être localisé sur une ou plusieurs phases.

[ ] = [ ] = R s [ ] ( 6 )
[ ] = [ ] [ ] ( 7 )
n sa correspond au nombre de spire de la phase a du stator.

[ ] [ ] [ ] ( 8 )
[ ] [ ] [ ] ( 9 )
[ ] = [ ] ( 10 )
B. Choix de la variable à étudier L'objectif de cette approche consiste à déterminer les signaux accessibles les plus significatifs pour la détection de défauts. Pour cela, il est nécessaire de tenir compte d'une part de la bande passante du système bouclé et d'autre part du contenu spectral du défaut à surveiller.

C. De ' f ue e e ou e e o o e

Il ne faut pas perdre de vue que le système comporte des boucles de contrôle des courants (pour assurer le contrôle vectoriel) et de la vitesse [START_REF] Diao | Sensor Fault Diagnosis for Improving the Availability of Electrical Drives[END_REF]. Ces correcteurs ont pour fonction d'assurer la conformité entre la consigne ref et la mesure y en agissent sur la grandeur de commande (u).Intuitivement, il est évident que la commande cherchera à compenser toute variation de la grandeur mesurées disponible pour la détection des défauts (les deux courants mesurés et la vitesse mesurées, le vecteur de tension de commande en sortie du dispositif de contrôle vectoriel et la référence de couple en sortie de l'asservissement de vitesse). Considérons le schéma bloc suivant représentant de manière simplifiée la partie contrôle du couple d'un entrainement électrique similaire à notre application.

D. Les signatures spéctrales du courant statorique

Le but de la caractérisation fréquentielle des défauts consiste, en fonction du défaut considéré, à déterminer la (ou les) fréquences (s) liées (s) au défaut et à discriminer la variable (mesurée ou calculée) la plus pertinente dans laquelle cette fréquence peut être isolée. La figure 2 montre la fréquence électrique f e =100Hz du spectre du courant statorique (Isa) de la phase a de la machine asynchrone sans court-circuit. La figure 3 montre le spectre du courant statorique (Isa) de la machine asynchrone avec court-circuit sur une phase, on note l'apparition de la fréquence du défaut 3f e (300Hz). La figure 4 présente les fréquences qui caractérisent le défaut par l'apparition de deux fréquences au tour de la fréquence électrique (1±2g) ou g représente le glissement.Remarquons que ce défaut sera tré difficile à détecter dans le spectre du courant car l'amplitude de la raie à la fréquence fondamentale est très importante par rapport à celle de la raie recherchée.L'amplitude de défaut considéré consiste en une suppression de 30% du bobinage d'une phase au rotor ou stator. Après avoir souligné le choix de la variable a consideré pour détecteter la signature fréquentielle carectéristique du défaut statorique et rotorique, nous avons présenté les caractéristiques fréquentielles des défauts concernant onduleur , capteur de courant et capteur de vitesse.La figure 5 illustre le spectre du courant statorique (Isa) de la machine saine avec défaut d'onduleur composant ouvert et la machine saine,on notant l'apparition des raies autour de la fréquence électrique f e et autours des fréquences de défauts kf e avec k =1,2,3,……. . Le défaut capteur de courant considéré consiste en un offset sur la grandeur mesuré (50A soit environ 15% du calibre).La figure3 présente le spectre du courant statorique (Isa) de la phase a de la machine lorsque le capteur est défaillant, en retrouve la fréquence de défaut 2f e (200Hz). Fig. 6.Spectre du courant de la phase a de la machine avec défaut sur le capteur de courant.

En tenant compte de la plage en basse fréquence, on remarque que la fréquence de défaut sur le capteur de vitesse est inexistante (les variables concernées sont continues). La figure 7 présente l'allure de la vitesse mesurée pour une erreur du capteur de 5 Rad/s (soit 2,5% de la valeur de référence).On remarque que le contrôle de la vitesse réagit de manière à compenser l'offset qui apparait sur la vitesse mesurée.

III. LES EFFETS MECANIQUES DES DEFAUTS DE ROULEMENTS

Afin de pouvoir utiliser un modèle assez développé de la partie mécanique de l'entrainement électrique, nous avons modélisé une transmission constituée par un arbre flexible avec réducteur. La partie mécanique contient aussi des roulements qui assurent le parfait alignement des axes des pièces tournantes. Nous avons considéré qu'ils n'intervenaient pas dans les équations mécaniques du système lorsqu'ils sont sains. En revanche, lorsqu'un défaut se produit sur ces éléments, nous les prenons en compte par la génération d'un couple de perturbation.

A. Modèle de la transmission

Le couple moteur est transmis à la roue par le biais d'un réducteur à double train épicycloïdal. Nous ne disposions pas de modelés mécanique de cette transmission. Nous avons donc développé un modèle simple permettant de présenter les deux phénomènes qui nous semblaient les plus importants à considérer : l'élasticité de la transmission et l'irrégularité du couple transmis due aux dents de l'engrènement.la figure 8 présente le modèle [START_REF] Faucher | Fault indicator in induction machine using stator current spectral analysis" Power Electronics, machine and drives[END_REF][START_REF] Raison | Détection et localisation de défaillance sur un entraiinement électrique[END_REF]. Les équations de cette transmission sont les suivantes :

1 = 1 d 1 dt f 1 1 1 (11) 1 = 2 d 2 dt f 2 2 2 (12) 1 = ( -2 ) ( 1 -2 ) (13)

B. P e e o p e e 'e g e ge

De nombreuse études ont modélisé l'engrenage par le biais de décomposition en séries de fourier et ont surtout analysé son influence sur les vibrations et les torsions d'arbre. Pour notre étude, nous avons préféré une approche plus intuitive afin d'expliquer les phénomènes physiques entrant en jeu. La validité des hypothèses et de la modélisation est justifiée a posteriori en comparant l'effet du modèle d'engrenage introduit sur la vitesse de rotation de l'axe de la machine sur des relevés issus d'autres travaux [START_REF] Ceban | éthode globale de diagnostic des machines électriques" hèse de doctorat[END_REF]. A partir de l'équation, nous pouvons remarquer que l'engrenage intervient au niveau du couple transmis (Ct).la forme de ce couple en fonction de la position est approchée par les modèles proposés à la figure. Les deux motifs présentés pour le couple transmis permettent de rendre compte des chocs entre les dents. Le premier motif, en pointillés, traduit les chocs entre chaque dent pour la transmission du couple entre les deux parties de l'engrenage. Le deuxième motif, en traits pleins, « représente » le frottement des dents en contact de chaque roue de l'engrenage. Nous avons choisi ce modèle pour la suite des simulations : il traduit une sorte de coefficient de torsion variable en fonction du nombre de dents en contact.

C. Les roulements

Une des raisons les plus fréquentes de défauts dans les machines électriques est liée aux défauts dans les roulements. Nous avons considéré qu'en fonctionnement sain les roulements ne généraient aucun signal. Un défaut sur un roulement entraine une modification de l'alignement de l'arbre de la machine avec l'axe de la machine. Il conduit donc à une modification de la répartition de l'induction au sein de la machine. Cependant, devant la complexité à rendre ces effets, nous avons décidé de ne représenter que le couple perturbateur généré sur l'axe de la machine [START_REF] Boqiang | weak-signal detection and the application in detection of electrical motors faults" Electrical achines and Systems[END_REF][START_REF] Jelassi | Impact d'un court-circuit interne au stator d'une machine asynchrone sur les phénomènes générés par la denture" hèse de doctorat[END_REF].

f = P f [1-( P cos ) 2 ] (15)
Avec BD =diamètre de la bille PD = profondeur du roulement = angle de contact de la bille sur la couronne interne -Pour un défaut sur la couronne interne ou externe du roulement :

{ f = n 2 f [1 - P cos ] f i = n 2 f [1 P cos ] (16)
Où n est le nombre de billes du roulement Nous présentons ici les résultats de simulations effectuées avec les modèles de roulements défaillants introduits. Les simulations ont été effectuées pour une vitesse de 150 rad/s (soit une fréquence mécanique de 23.9 Hz). Les figures suivantes présentent les spectres de la référence de couple (sortie du correcteur de vitesse) pour les différents défauts roulement considérés.

IV. CONCLUSION:

Apres avoir souligné l'influence du contrôle sur le choix de la variable à considérer pour détecter la signature fréquentielle A partir de la vitesse mécanique Les modéles de machine utilisés pour la caractérisation et la validation des algorithmes de surveillance de l'application sont critiquables bien qu'ils permettent une représentation plus réaliste des situations de défaut que les modéles équilibrés couramment adoptés dans la littérature.Ils permettent de représenter des déséquilibres et l'influence des défauts sur l'ensemble des grandeurs. Une carecterisation de degré de validité du modèle développé est à faire.La difficulté rencontrée avec les méthodes de traitement du signal tient à la mise en évidence des raies dont la fréquence et l'amplitude varient avec le point de fonctionnement.

Fig. 1 .

 1 Fig.1. Schéma bloc du système simplifie

Fig. 2 .

 2 Fig.2. Spectre du courant statorique de la phase a.

Fig. 3 .

 3 Fig.3. Spectre du courant statorique avec défaut sur une phase au stator.

Fig. 4 .

 4 Fig.4. Spectre du courant statorique avec défaut sur une phase au rotor.

Fig. 5 .

 5 Fig.5. Spectre du courant statorique de la machine saine avec défaut onduleur et la machine saine.

Fig. 7 .

 7 Fig.7. Vitesse mesurée lors d'un défaut capteur de vitesse.

Fig. 8 .

 8 Fig.8. Modèle d'une transmission flexible

Fig. 9 .

 9 Fig.9. Effet sur le couple transmis par l'engrenage

  Application du défaut capteurOù est la position de l'axe. Tm, amplitude du défaut, est choisie proportionnelle au couple appliqué sur l'arbre. K m est donné par k m = f f où f est la vitesse mécanique du rotor et est f la fréquence correspondant au défaut choisi. f est une fonction des dimensions du roulement et devient : Pour un défaut sur une bille :

Fig. 1 .

 1 Fig.1 .Grandeurs caractéristiques d'un roulement

Fig. 11 .

 11 Fig.11. Spectre de la référence de couple pour un défaut sur la bague externe.

Fig. 12 .

 12 Fig.12. Spectre de la référence de couple pour un défaut sur la bague interne.

Fig. 13 .

 13 Fig.13. Spectre de la référence de couple pour un défaut sur la bille.

Table 1 : Les caractéristiques des défauts électriques

 1 caractéristique du défaut, nous avons présenté les caractéristiques fréquentielles des défauts concernant notre application. Ces fréquences peuvent être variables (fonction de la fréquence électrique ou de la fréquence de rotation mécanique) ou fixe (liées à la MLI). Elles sont résumées dans les tableaux suivants.A partir du courant de phase

		Stator	Rotor	Onduleur	Capteur courant
	fe
	2fe			X	X
	3fe	X		X

Table 2 : Les caractéristiques des défauts mécaniques

 2

		Défaut roulement	Défaut roulement	Défaut roulement
		1 (bague externe)	2 (bague interne)	3 (bille)
	f0	X
	fi		X
	fb			X


Francesca Branzoli
email: f.branzoli-ihu@icm-institute.org

Ece Ercan
email: eceercan@gmail.com

Romain Valabrègue
email: romain.valabregue@upmc.fr

Emily T Wood
email: emilytwood@gmail.com

Mathijs Buijs
email: m.buijs@lumc.nl

Andrew Webb
email: a.webb@lumc.nli.ronen@lumc.nl

Itamar Ronen

Differentiating between

Keywords: diffusion-weighted MRS, intracellular metabolites, DTI, healthy aging, white matter

à 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.

Introduction

Diffusion tensor imaging (DTI) has provided much evidence of age-related changes in human brain white matter (WM) microstructure, see e.g. (Salat, et al., 2005,Sullivan and[START_REF] Sullivan | Diffusion tensor imaging and aging[END_REF].

Lower fractional anisotropy (FA) and higher water mean diffusivity (MD) have been consistently detected in several white matter regions, suggesting alteration of WM microstructure in healthy aging, even when concomitant volume declines are not detectable [START_REF] Coutu | Non-Gaussian water diffusion in aging white matter[END_REF][START_REF] Madden | Diffusion tensor imaging of cerebral white matter integrity in cognitive aging[END_REF][START_REF] Sullivan | Diffusion tensor imaging and aging[END_REF]. DTI analysis has suggested that the primary neurobiological mechanism associated with age-related degeneration in WM is a decline in the structural integrity of myelin [START_REF] Wozniak | Advances in white matter imaging: a review of in vivo magnetic resonance methodologies and their applicability to the study of development and aging[END_REF]. This hypothesis was corroborated by histological analyses performed in the rhesus monkey as well as in human post-mortem samples, showing that age-related alterations in the structure of myelin sheaths precede the loss of myelinated nerve fibers [START_REF] Murray | A quantitative postmortem MRI design sensitive to white matter hyperintensity differences and their relationship with underlying pathology[END_REF][START_REF] Peters | The effects of normal aging on myelin and nerve fibers: a review[END_REF][START_REF] Peters | A review of the structural alterations in the cerebral hemispheres of the aging rhesus monkey[END_REF][START_REF] Sandell | Disrupted myelin and axon loss in the anterior commissure of the aged rhesus monkey[END_REF]. Other studies have suggested that selective axonal loss occurs in healthy aging, providing evidence of "resilience" of thick axons with thick myelin, combined with the loss of thin, poorly-myelinated axons [START_REF] Marner | Marked loss of myelinated nerve fibers in the human brain with age[END_REF][START_REF] Tang | Age-induced white matter changes in the human brain: a stereological investigation[END_REF]. Ex vivo studies also pointed out that, despite alterations in the myelin sheaths, only a few nerve fibers showed axonal degeneration [START_REF] Peters | Aging and the myelinated fibers in prefrontal cortex and corpus callosum of the monkey[END_REF]. While FA is a sensitive marker for WM injury without specificity for axonal injury or demyelination, it has been suggested that increases in radial diffusivity and decreases in axial diffusivity correlate with demyelination processes and axonal damage, respectively [START_REF] Song | Dysmyelination revealed through MRI as increased radial (but unchanged axial) diffusion of water[END_REF][START_REF] Wang | Diffusion basis spectrum imaging detects and distinguishes coexisting subclinical inflammation, demyelination and axonal injury in experimental autoimmune encephalomyelitis mice[END_REF]. Support for this hypothesis was also provided by histological analysis [START_REF] Budde | Axial diffusivity is the primary correlate of axonal injury in the experimental autoimmune encephalomyelitis spinal cord: a quantitative pixelwise analysis[END_REF][START_REF] Hofling | Diffusion tensor imaging detects axonal injury and demyelination in the spinal cord and cranial nerves of a murine model of globoid cell leukodystrophy[END_REF][START_REF] Kim | Detecting axon damage in spinal cord from a mouse model of multiple sclerosis[END_REF]. Although DTI metrics are potentially sensitive measures of alterations in tissue microstructure, they are indirect measures of tissue integrity. Since DTI reflects the molecular diffusion of water, which is present in all tissue compartments, it is influenced by several microstructural factors that affect both the intraas well as the extracellular space, including myelin and axonal density, membrane integrity and axonal orientation dispersion [START_REF] Beaulieu | The basis of anisotropic water diffusion in the nervous system -a technical review[END_REF][START_REF] Douaud | DTI measures in crossing-fibre areas: increased diffusion anisotropy reveals early white matter alteration in MCI and mild Alzheimer's disease[END_REF]. As a consequence, the interpretation of DTI metrics based on the neurobiological mechanisms that determine their values is challenging. The interpretation of DTI data is further complicated by the local architecture of WM tracts. For example, in the case of crossing fibers, the major axis of the tensor model may not be parallel to the myelinated axons, the many separate tracts which propagate within the DTI voxel result in low local FA values (Pierpaoli, et al., 2001,Wheeler-Kingshott and[START_REF] Wheeler-Kingshott | About "axial" and "radial" diffusivities[END_REF].

Diffusion-weighted 1 H magnetic resonance spectroscopy (DW-MRS) provides structural information exclusively related to the intracellular space by probing the diffusion properties of intracellular metabolites, such as total N-acetylaspartate (tNAA = N-acetylaspartate + N-acetyl aspartyl glutamate), total creatine (tCr = creatine + phosphocreatine), and choline compounds (tCho) [START_REF] Nicolay | Diffusion NMR spectroscopy[END_REF][START_REF] Posse | Human brain: proton diffusion MR spectroscopy[END_REF]. The distribution of these metabolites across different cell types in the brain makes DW-MRS not only compartment-specific, but also specific or at least preferential for different cell types [START_REF] Choi | Application of MRS to mouse models of neurodegenerative illness[END_REF][START_REF] Le Belle | A comparison of cell and tissue extraction techniques using high-resolution 1H-NMR spectroscopy[END_REF][START_REF] Urenjak | Proton nuclear magnetic resonance spectroscopy unambiguously identifies different neural cell types[END_REF]. NAA and NAAG are known to be present almost exclusively inside neurons, and in white matter they are exclusively intra-axonal. tNAA diffusion properties in white matter are therefore solely dictated by the structural properties of intra-axonal components and by the properties of the intra-axonal medium, thus providing a specific marker to the axonal contribution to WM degeneration. Soluble Choline compounds are located mainly in glial cells [START_REF] Choi | Application of MRS to mouse models of neurodegenerative illness[END_REF][START_REF] Le Belle | A comparison of cell and tissue extraction techniques using high-resolution 1H-NMR spectroscopy[END_REF][START_REF] Urenjak | Proton nuclear magnetic resonance spectroscopy unambiguously identifies different neural cell types[END_REF] and thus their diffusion properties report on structural changes in glia, resulting from, for example, inflammation [START_REF] Ercan | Glial and axonal changes in systemic lupus erythematosus measured with diffusion of intracellular metabolites[END_REF]. Creatine and phosphocreatine are present in both tissue types and are tightly linked to cell energetics through the creatine kinase cycle (CK), and it has been shown that increase in energy demand modulates the diffusion coefficient of total creatine [START_REF] Branzoli | Functional diffusionweighted magnetic resonance spectroscopy of the human primary visual cortex at 7 T[END_REF]. Previous observations have provided strong evidence for the important physiological role of the CK system for normal brain function, and suggested that energy dysfunction processes are crucially related to cell death not only in neurodegenerative disorders but even in normal aging [START_REF] Andres | Functions and effects of creatine in the central nervous system[END_REF]. So far, DW-MRS has been applied to study normal aging [START_REF] Zheng | The effect of age and cerebral ischemia on diffusion-weighted proton MR spectroscopy of the human brain[END_REF] as well as disease pathophysiology in brain tumors [START_REF] Colvin | New insights into tumor microstructure using temporal diffusion spectroscopy[END_REF][START_REF] Harada | Diffusion-weighted in vivo localized proton MR spectroscopy of human cerebral ischemia and tumor[END_REF], acute cerebral ischemia [START_REF] Harada | Diffusion-weighted in vivo localized proton MR spectroscopy of human cerebral ischemia and tumor[END_REF][START_REF] Zheng | The effect of age and cerebral ischemia on diffusion-weighted proton MR spectroscopy of the human brain[END_REF], multiple sclerosis [START_REF] Wood | Investigating axonal damage in multiple sclerosis by diffusion tensor spectroscopy[END_REF], psychiatric disorders [START_REF] Du | Myelin and axon abnormalities in schizophrenia measured with magnetic resonance imaging techniques[END_REF], and mitochondrial cytopathies [START_REF] Liu | Apparent diffusion coefficients of metabolites in patients with MELAS using diffusion-weighted MR spectroscopy[END_REF]. Here we employed DW-MRS in combination with DTI in order to separately detect and quantify intra-axonal damage and demyelination in healthy aging in the anterior body of the corpus callosum (CC). The integrity of the CC plays a crucial role in cognitive and motor decline, since processes involved in crosshemispheric information exchange depend on the integrity of the callosal fibers (Banich, 1998,Reuter-Lorenz and[START_REF] Reuter-Lorenz | Differential effects of aging on the functions of the corpus callosum[END_REF]. Differences in the DTI metrics measured in the corpus callosum of healthy elderly have been reported in several DTI studies, yet the underlying biological mechanisms for these changes are still unknown [START_REF] Billiet | Age-related microstructural differences quantified using myelin water imaging and advanced diffusion MRI[END_REF][START_REF] Fling | Differential callosal contributions to bimanual control in young and older adults[END_REF][START_REF] Teipel | Longitudinal changes in fiber tract integrity in healthy aging and mild cognitive impairment: a DTI follow-up study[END_REF]. In addition, the CC is a prominent white matter structure that is easy to identify on MRI images, and thus positioning of the volume of interest (VOI) for DW-MRS experiments in a consistent manner across subjects is relatively straightforward. The morphology of the CC allows a relatively easy estimate of the pure intra-axonal tNAA diffusivity, as well as accurate volumetric measurements, of interest as macroscopic correlates of the microstructural indices derived from the DW-MRS and DTI experiments.

In this study, we used DTI and DW-MRS to measure the diffusion properties of water, tNAA, tCr, and tCho in the CC of healthy elderly and young subjects at 7T. We combined the information from these two MR tools to discriminate between age-related myelin alterations and possible concomitant axonal or glial structural changes, and to obtain insight into possible cell energy alterations in normal aging. For an unbiased estimation of the intra-axonal diffusivity of tNAA, we analyzed the DW-MRS data using a model that provides the cytosolic diffusion coefficient of tNAA, free from confounds stemming from the curvature of the CC within the VOI and from microscopic fiber orientation dispersion [START_REF] Ronen | Microstructural organization of axons in the human corpus callosum quantified by diffusionweighted magnetic resonance spectroscopy of N-acetylaspartate and post-mortem histology[END_REF]

MRI/MRS protocol

Anatomical images: A short survey scan and a sensitivity encoding (SENSE) reference scan were followed by a 3D T 1 -weighted gradient echo acquisition to allow for accurate planning of the spectroscopy volume of interest (VOI). Imaging parameters were: field-of-view (AP, FH, LR):

246.4x246.4x174.0 mm 3 , resolution 1x1x1 mm 3 , TR/TE = 4.7 ms/2.1 ms.

DTI protocol:

Multi-slice single-shot 2D spin-echo echo planar imaging (EPI) was performed with the following data acquisition parameters: 60 slices of 2 mm thickness, TR/TE = 9100 ms/58 ms, field of view of 224×224×120 mm 3 with an in-plane data matrix of 112×112, resulting in a data set with resolution of 2×2×2 mm 3 , 15 diffusion weighting directions with b = 1000 s/mm 2 . Fat suppression was performed using spectral inversion recovery (SPIR). Parallel imaging was performed using SENSE with a reduction factor of 3 along the phase encoding direction (AP).

DW-MRS protocol:

A single voxel, point resolved spectroscopy (PRESS) sequence with a bipolar diffusion weighting gradient scheme was used. A VOI of 30 (AP) x 15 (LR) x8 (FH) mm 3 was positioned in the anterior part of the body of the CC (Insets of Figure 1). Sequence parameters were: TR = 3 cardiac cycles (cardiac triggering was performed using the peripheral pulse unit), TE = 121 ms, spectral width SW = 3 kHz, and number of acquired complex sample points NP = 2048.

Diffusion weighting was applied in an interleaved way in two directions with respect to the VOI: left-right (g (0,0,1) , mostly parallel to the callosal fibers), and one that formed a 45° angle between the superior-inferior and anterior-posterior directions of the VOI (g (1,1,0) , mostly perpendicular to the callosal fibers). The VOI position and the gradient directions are shown in Figure 1 (insets). The diffusion weighting was applied using a diffusion time Δ = 60.5 ms, duration of bipolar diffusion group: δ = 34 ms , bipolar gap: 18 μs and 4 increasing gradient strengths g corresponding to b values 330, 960, 1900, 3200 s/mm 2 , and 660, 1920, 3800, 6400 s/mm 2 for the (0,0,1) and (1,1,0) directions, respectively. 32 spectra were acquired for each direction and b value. Water suppression was applied using a frequency selective excitation pulse followed by a dephasing gradient before metabolite excitation. The water signal was only partially suppressed to allow for a sufficient residual water signal for individual phasing of each spectrum. Non-water suppressed spectra were acquired under similar conditions for eddy current corrections. The total duration of the DW-MRS protocol was approximately 15 minutes.

Post-processing and analysis

Image processing: T 1 -weighted images were used for tissue segmentation within the VOI, and further processed in FreeSurfer (http://surfer.nmr.mgh.harvard.edu/) to calculate the intracranial volume, total brain volume, and CC volumes for each subject. To evaluate the callosal atrophy due to aging, the CC volumes were normalized according to the intracranial volume of the same subject.

The maps of the main eigenvector (e 1 ) were derived and confined to the VOI using MIPAV (Medical Image Processing, Analysis and Visualization) software, JIST (Java Image Science Toolbox) and in-house developed MATLAB® codes (Mathworks Inc., Natick MA, USA) as described in [START_REF] Wood | Reproducibility and optimization of in vivo human diffusion-weighted MRS of the corpus callosum at 3 T and 7 T[END_REF]. These maps were used in the modeling procedure described below.

Average DTI volumes were motion corrected with ExploreDTI, and DTI metrics were derived using the DTI toolbox of the FMRIB Software Library (FSL release 5.0, http://www.fmrib.ox.ac.uk/fsl/). The DTI metrics were further analyzed with tract based spatial statistics (TBSS) [START_REF] Smith | Tractbased spatial statistics: voxelwise analysis of multi-subject diffusion data[END_REF]. Statistical differences between older and younger subjects were assessed in FA-MNI152 standard space using 5000 permutations and were corrected for multiple comparisons based on threshold-free cluster enhancement (Winkler, et al.). Three young and one old volunteer were excluded from TBSS analysis due to the poor registration to FA-MNI152. T 1 -volumes were registered to MNI152 space first by affine transformation using FSL FLIRT (Jenkinson, et al., 2002,Jenkinson and[START_REF] Jenkinson | A global optimisation method for robust affine registration of brain images[END_REF], and then by a non-rigid registration using FNIRT [START_REF] Andersson | Non-linear registration, aka spatial normalisation[END_REF]. The same transformation matrices were then used to register the DW-MRS VOI to the MNI152 standard space.

DW-MRS processing:

The diffusion-weighted spectra were quantified using LCModel [START_REF] Provencher | Automatic quantitation of localized in vivo 1H spectra with LCModel[END_REF]) using an appropriate basis set. Diffusion coefficients for tNAA in the two gradient directions, D (0,0,1) (tNAA) and D (1,1,0) (tNAA), were estimated by fitting the tNAA peak estimates from LCModel to monoexponential decay functions as shown for one subject in Figure 2. D (0,0,1) (tCr), D (1,1,0) (tCr), and D (0,0,1) (tCho), D (1,1,0) (tCho) were estimated similarly. The cytosolic diffusion coefficient of tNAA, D(tNAA), was estimated by fitting the tNAA LCModel data to a model described in detail in [START_REF] Ronen | Microstructural organization of axons in the human corpus callosum quantified by diffusionweighted magnetic resonance spectroscopy of N-acetylaspartate and post-mortem histology[END_REF]. Briefly, the model assumes the tNAA to be exclusively intra-axonal, and the axons to be curved cylinders with a known diameter distribution [START_REF] Aboitiz | Fiber composition of the human corpus callosum[END_REF]. The effect of the macroscopic curvature of the CC within the VOI on the diffusionweighted data in both diffusion weighted directions is accounted for by using the distribution of angles between the main eigenvectors e 1 of the DTI data and the diffusion weighting directions. The experimental data is fitted to the model using two fitting variables: D(tNAA), which is independent of the curvature of the tract within the VOI, and the microscopic axonal angular dispersion φ, which describes the average fanning of individual axons about the main propagation direction within a DTI voxel. Student t-test were performed in order to evaluate possible statistical differences between the two groups. The significance level was set to p = 0.05.

Results

Table 1 summarizes the demographic characteristics of young and healthy elderly subjects, as well as the average voxel tissue composition of the two groups, and the CC normalized volumes. Small, not statistically significant differences, were observed between the two groups in the WM fraction (young subjects: WM fraction ~ 64%, elderly subjects: WM fraction ~ 68%) and CSF fractions (young subjects: WM fraction ~ 29%, elderly subjects: WM fraction ~ 24%), while differences in GM fraction were negligible (GM fraction ~ 8% in both groups). However, when the WM and GM fractions were normalized over the brain tissue percentage (100 -CSF(%)), no differences could be observed in the tissue composition (WM ~ 90%, GM ~ 10% in both groups). The difference in normalized callosal volume between the two groups was not statistically significant (p = 0.11).

Figure 3a shows an FA map of a coronal section of the corpus callosum calculated for one of the younger subjects, overlaid with the representation of the main eigenvectors of the diffusion tensor for each voxel, and with the DWS VOI. Figure 3 b andc show histograms of the distributions of the angles  1 and  2 of the main eigenvector e 1 and the gradients g (0,0,1) and g (1,1,0) , respectively, calculated within the masked VOI for the same subject. Histogram bins are of 1°. The two distributions show the highest probability around 13° and around 95°, for  1 and  2 respectively, consistently with the CC morphology.

Figure 4 show maps of statistically significant differences in DTI measures overlaid on the MNI152

T 1 -weighted image, and the cumulative DW-MRS VOI (i.e. the sum of all the VOIs for each patient, following the transformation to MNI152 coordinates). Significantly lower FA, higher MD and higher radial diffusivity (RD) in the older subjects group compared to the younger subjects group were found throughout white matter, including within the DW-MRS VOI. In contrast, no significant differences in AD were measured in the VOI.

Examples of diffusion-weighted spectra acquired in one subject with diffusion gradients applied in the directions parallel and perpendicular to the CC fibres at low and high b values are shown in Figure 1. The signal decay induced by the diffusion weighting in the direction parallel to the callosal fibres was much steeper than when the gradient was applied in the direction perpendicular to the fibres. This is due to the fact that the molecular motion along the axonal direction is essentially unrestricted, as opposed to the motion perpendicular to the axonal membrane. The linewidth of tNAA was ~ 12 Hz, and the Cramer-Rao lower bound (CRLB) was < 15 % for all diffusion weighting conditions for tNAA and tCr, and < 20% for tCho. Figure 2 shows the tNAA normalized signal decay as derived from LCModel analysis, plotted on a logarithmic scale, as a function of the b value, for the two directions (0,0,1) and (1,1,0) (closed and open symbols, respectively). Both the modelled free tNAA diffusion coefficient, D(tNAA), and the tNAA diffusion coefficient measured in the direction of the callosal fibers, D (0,0,1) (tNAA), were found to be slightly smaller in healthy elderly subjects compared to young controls, but these differences did not reach statistical significance. Neither the tNAA diffusion coefficient measured in the direction perpendicular to the fibers, D (1,1,0) (tNAA), nor the the axonal microscopic angular dispersion angle φ differed between the two groups.

The tCr axial diffusion coefficient D (0,0,1) (tCr) was found to be smaller in healthy elderly subjects (p = 0.058), while the other tCr diffusion parameters did not change between the two groups.

None of the tCho diffusion properties was found to differ significantly between the two groups.

Table 2 summarizes the results for metabolite diffusion parameters in young and healthy elderly subjects.

Discussion

In this study we used a combination of tools sensitive to the microstructural features of neural tissue to investigate age-related changes in normal-appearing white matter in healthy aging. The novel feature we introduced here is the incorporation of DW-MRS, a compartment and cell-type measure for specific alterations in the intracellular space of axons and glia. The observed differences in water MD (higher in older adults) and in water FA (lower in older adults) in the CC are consistent with previous studies suggesting the onset of demyelination processes in the seventh decade of age or later in healthy elderly (Salat, et al., 2005,Sullivan and[START_REF] Sullivan | Diffusion tensor imaging and aging[END_REF]. The estimated water RD values were found to be greater in elderly subjects, while no significant differences in the water AD were observed in the DWS voxel, also in accordance with literature reports [START_REF] Burzynska | Age-related differences in white matter microstructure: region-specific patterns of diffusivity[END_REF][START_REF] Coutu | Non-Gaussian water diffusion in aging white matter[END_REF][START_REF] Madden | Diffusion tensor imaging of cerebral white matter integrity in cognitive aging[END_REF]. The observed differences in water diffusion properties of the elderly cannot be solely related to a loss of axonal fibers leading to an increase in the extracellular space, since the older adults did not show evidence of severe atrophy of the CC, as expected in subjects in this age range (60-70 years), and the differences in callosal volumes were below the threshold of statistical significance. In some previous studies the DTI results were interpreted as a consequence of the fact that in healthy aging demyelination processes are dominant with respect to other possible microstructural WM damages such as disruption of intra-axonal components, based on the hypothesis that water AD changes should reflect ongoing intra-axonal damage [START_REF] Madden | Diffusion tensor imaging of cerebral white matter integrity in cognitive aging[END_REF]. However, the assumption that the water AD can be unequivocally associated with degeneration processes occurring in the intra-cellular space as suggested previously [START_REF] Hofling | Diffusion tensor imaging detects axonal injury and demyelination in the spinal cord and cranial nerves of a murine model of globoid cell leukodystrophy[END_REF] seems to be rather simplistic if the complexity of fiber organization is taken into account. Even in a region of the WM such as the CC where areas with no macroscopic intersection of fiber bundles can be isolated, the presence of microscopic fiber bending and fiber crossing makes the interpretation of the AD parameter more complicated, since it may also be susceptible to demyelination effects (Ronen et al., 2013a;[START_REF] Wheeler-Kingshott | About "axial" and "radial" diffusivities[END_REF].

The hypothesis that water AD alone cannot be interpreted as reflecting purely intra-axonal changes was recently suggested in a study by Wood et al., where an increase in water AD and a concomitant decrease in tNAA axial diffusivity were observed in the CC of multiple sclerosis (MS) patients using a methodology similar to the one presented here. In Wood et al., the diffusion tensor of tNAA was calculated to yield, among the other metrics, the axial diffusivity of tNAA. Although this measure could be partially affected by the macroscopic curvature of the CC in the spectroscopic VOI, as well as by the microscopic axonal misalignment in the DTI voxel, an opposite trend was observed in water and tNAA axial diffusivity, showing that these two metrics are not markers of the same biological processes. This point is further emphasized by the confirmation here of the existence of significant axonal misalignment (φ ~ 12 deg), already shown by DW-MRS and histology (Ronen et al., 2013a), and known to affect DTI measurements as well [START_REF] Sotiropoulos | Fusion in diffusion MRI for improved fibre orientation estimation: An application to the 3T and 7T data of the Human Connectome Project[END_REF].

In the current study, we estimated the cytosolic tNAA diffusion coefficient within axons, which is not affected by the interactions with the cellular walls, or by fiber geometry, and it is thus a pure measure of intra-axonal integrity. Our results are the first to show a decrease in D (0,0,1) (tNAA), as well as a decrease in D(tNAA), in an elderly population; however, the results did not reach statistical significance. These differences, if confirmed, may be related to the presence of structural disruptions associated with axonal degeneration, or pruning of thinner axons, and the very small effect is in line with the fact that our group included healthy aging subjects in the seventh decade.

The population size in this study was similar to that in a previous DW-MRS study in MS that investigated axonal differences in normal appearing white matter in the CC [START_REF] Wood | Investigating axonal damage in multiple sclerosis by diffusion tensor spectroscopy[END_REF]. In both studies, significant effects on DTI measures such as FA and RD were detected. In our current study, an optimized acquisition scheme was employed that corrects for inter-subject variability in the CC curvature using a previously demonstrated model. Yet, the differences in D(tNAA) in the CC between the two age populations were negligible. Our results suggest that the intra-axonal damage in the CC of healthy subjects in this age range is small (<10%), in line with the results derived by Peters and Sethares in 2002 from histological measurements in normally-aging macaques, and supporting the idea that cognitive impairments are more strongly associated with axonal integrity rather than with demyelination processes [START_REF] Sandell | Disrupted myelin and axon loss in the anterior commissure of the aged rhesus monkey[END_REF].

The creatine axial diffusivity D (0,0,1) (tCr) was found to be lower in healthy elderly subjects compared to younger controls but this finding did not reach statistical significance (p = 0.058).

This finding may be linked to the processes involving adenosine triphosphate (ATP) production and utilization in the brain, which are closely connected to the conversion of PCr to Cr through the CK cycle [START_REF] Wallimann | Intracellular compartmentation, structure and function of creatine kinase isoenzymes in tissues with high and fluctuating energy demands: the 'phosphocreatine circuit' for cellular energy homeostasis[END_REF]. By enabling saltatory conduction, myelination is also thought to reduce ATP consumption and, in turn, to be related to axonal energy balance [START_REF] Saab | The role of myelin and oligodendrocytes in axonal energy metabolism[END_REF]. A lower ATP production rate in aging brains was further corroborated by a study using phosphorus MRS at 7T [START_REF] Zhu | In vivo NAD assay reveals the intracellular NAD contents and redox state in healthy human brain and their age dependences[END_REF]. Recent studies have convincingly shown that glial cells, especially oligodendrocytes at the axonal level and astrocytes at the synaptic level, supply lactate to neurons as a metabolic substrate for mitochondria to generate ATP, and that this process is altered by demyelination [START_REF] Saab | The role of myelin and oligodendrocytes in axonal energy metabolism[END_REF]. The demyelination processes affecting the white matter of healthy elderly subjects may, therefore, be indirectly reflected in abnormalities in the diffusion of tCr as a consequence of an impaired energy metabolism, possibly yielding an increase in the PCr/Cr ratio, or a dysfunction in the trafficking of these metabolites between glia and neurons [START_REF] Andres | Functions and effects of creatine in the central nervous system[END_REF].

Taken together, our results on tNAA and tCr diffusion support the hypothesis that the two major oligodendrocyte functions of myelination and axonal support are distinct, as suggested by a previous study where no perturbation of axonal integrity was observed in a model of truly dysmyelinated mutant mouse [START_REF] Griffiths | Axonal swellings and degeneration in mice lacking the major proteolipid of myelin[END_REF].

The high specificity of MRS techniques is provided at the cost of an intrinsic low signal-to-noise ratio (SNR), which is further exacerbated by the addition of diffusion weighting and error propagation in the calculation of diffusion measures. The reproducibility of the model employed here for calculating metabolite diffusion metrics was tested previously in a study investigating the inter-and intra-subject variability of both empirical and modelled diffusion properties of tNAA, and evaluating the optimized protocol to be used in clinically relevant acquisition times in the same region of interest [START_REF] Wood | Reproducibility and optimization of in vivo human diffusion-weighted MRS of the corpus callosum at 3 T and 7 T[END_REF]. The acquisition scheme employed in this study was carefully evaluated following the guidelines reported in [START_REF] Wood | Reproducibility and optimization of in vivo human diffusion-weighted MRS of the corpus callosum at 3 T and 7 T[END_REF]. Power calculations in the same study showed that, based on the suggested acquisition scheme, at 7T the minimum number of subjects in each group necessary to detect a 10% difference in D(tNAA), ADC(tNAA) and ADC(tCr) between two groups is 13, 12, and 7, respectively. The acquisition scheme employed in the current study is expected to give less variability than the one in [START_REF] Wood | Reproducibility and optimization of in vivo human diffusion-weighted MRS of the corpus callosum at 3 T and 7 T[END_REF], as it employs similar number of averages with one more b-value. This emphasizes the accuracy and reproducibility of DW-MRS measurements at 7T.

One of the limitations of this study is that, due to the non-negligible chemical shift displacement related to the relatively long refocusing RF pulses in the PRESS sequence, the VOI localization was optimized for tNAA, but not for tCr. Therefore, D (0,0,1) (tCr) may be contaminated by D (1,1,0) (tCr) to a larger extent than in the case of tNAA, due to the imperfect alignment of the VOI with the fibers.

Another limitation is the impossibility of excluding completely grey matter contamination from the VOI, due to the shape of the corpus callosum and the minimum size of the VOI required to achieve a reasonable spectral signal-to-noise ratio. Although an average fraction of GM of about 8 % was present in the VOI of both young and elderly, it is very likely that this may introduce a small bias in the estimation of the metabolite diffusion properties, albeit without affecting the comparison between groups. Additional factors that may have impacted on the statistical significance of our results, and that should be improved in future studies, are the relatively small number of subjects who were recruited, the sub-optimal resolution for the DTI acquisition employed due to time constraints, and the localized spectroscopy method used. Improvements to the latter could be implemented using, for example, a diffusion-sensitized semi-localized by adiabatic selective refocusing (SEMI-LASER) technique that would result in a reduced chemical shift displacement (Scheenen et al.) and B 1 field inhomogeneity compared to the PRESS technique. Although 7T allows higher SNR and higher resolution, due to the limited time allowed for this study the DTI resolution was sacrificed in order to allow enough time for an optimized DW-MRS acquisition and the other scans included in the protocol.

Conclusion

In this work we explored age-related microstructural changes in the corpus callosum in vivo with an approach that combines diffusion tensor imaging (DTI) with diffusion-weighted MRS (DW-MRS) at ultrahigh field. The combination of the overall sensitivity to microstructural changes provided by DTI, and the specificity of DW-MRS to the intracellular space, provides a unique possibility for invivo characterization of subtle compartment-specific changes in tissue microstructure in normal aging. Furthermore, the specificity or preferentiality of certain metabolites to either neurons/axons or to glia endows this approach with an additional level of characterization, unmet by other MRbased methods. Our results show that the intra-axonal microstructural changes as measured by tNAA diffusivity, if present, are below the level of detectability with our DW-MRS protocol in healthy elderly in the age range investigated, suggesting that the intra-axonal integrity is mainly preserved in this population. This result requires confirmation in further studies employing larger population groups for establishing subtle differences in these intracellular measures. Abnormal changes in the tCr diffusivity may be interpreted as a reduction in the energy metabolism related to an ongoing CK dysfunction and coupled to myelin synthesis and transport. Further experiments are nevertheless needed to corroborate this hypothesis. Our findings will be substantiated with preclinical studies of animal models of disease combined with histology for comparison of DW-MRS metrics with a gold standard. The signal decay induced by the diffusion weighting in the direction parallel to the callosal fibers (0,0,1) was much steeper than when the gradient was applied in the direction perpendicular to the fibers (1,1,0). This is due to the fact that the molecular motion along the axonal direction is essentially unrestricted, as opposed to the motion perpendicular to the axonal membrane. Solid lines represent fits to single exponentials.

Error bars represent the Cramer-Rao Lower Bounds (CRLB) derived with LCModel.

FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: tNAA signal decay induced by diffusion weighting. tNAA normalized signal decay,

Figure 3 :

 3 Figure 3: Macroscopic fiber orientation in the corpus callosum. a) Coloured FA map of a coronal

Figure 4 :

 4 Figure 4: Differences in DTI metrics in healthy elderly versus younger controls. TBSS results

  ).

	2 Materials and methods
	2.1 Human subjects
	Thirteen healthy young controls (mean age 27 ± 3 years, age range 23-32 years, 8 males, 5 females)
	and sixteen healthy elderly subjects (mean age 64 ± 3 years, age range 60-71 years, 10 males, 6
	females) gave informed consent according to the procedures approved by the Leiden University
	Medical Center Institutional Review Board. Demographic characteristics are reported in Table 1.
	Subjects had no known history of neurological or psychiatric disorders.
	2.2	MRI scanner/Hardware
	Scans were performed on a 7 Tesla Achieva Philips whole-body MRI scanner (Philips Healthcare,
	Best, The Netherlands) equipped with gradient coils capable of a maximum gradient strength of 40
	mT/m and a slew rate of 200 T/m/s. A head coil consisting of a quadrature birdcage transmit and
	32-channel phased array receive were used for all measurements (Nova Medical Inc., Wilmington,
	MA, USA).

Table 1 :

 1 Demographic characteristics, spectroscopic voxel tissue composition and normalized callosal volumes of healthy elderly subjects and healthy young controls. Values are the mean ± standard deviation.

		Young	Elderly	p-value
	Age (years)	27 ± 3	64 ± 3	1
	Gender (male/female)	8/5	10/6	0.96
	Voxel WM composition (%)	63 ± 8	68 ± 7	0.11
	Voxel GM composition (%)	8 ± 3	8 ± 3	0.95
	Voxel CSF composition (%)	28 ± 8	24 ± 7	0.10
	CC volume (%)	0.24 ± 0.02	0.22 ± 0.03	0.11

Table 2 :

 2 tNAA, tCr, and tCho group average diffusion properties and standard deviation of the microscopic axonal orientation distribution measured in the spectroscopic voxel. Values are the mean ± standard deviation.

		Young	Elderly	p value
	D(tNAA) (µm 2 /ms)	0.53 ± 0.04	0.51 ± 0.04	0.21
	D (0,0,1) (tNAA) (µm 2 /ms)	0.39 ± 0.04	0.37 ± 0.04	0.12
	D (1,1,0) (tNAA) (µm 2 /ms)	0.08 ± 0.02	0.08 ± 0.02	0.79
	D (0,0,1) /D (1,1,0) (tNAA)	4.9 ± 1.1	4.8 ± 1.4	0.78
	ADC(tNAA)	0.24 ± 0.02	0.23 ± 0.02	0.18
	σ φ (deg)	12 ± 5	12 ± 5	0.77
	D (0,0,1) (tCr) (µm 2 /ms)	0.27 ± 0.04	0.24 ± 0.04	0.058
	D (1,1,0) (tCr) (µm 2 /ms)	0.11 ± 0.03	0.11 ± 0.03	0.96
	D (0,0,1) /D (1,1,0) (tCr)	2.5 ± 0.6	2.2 ± 0.8	0.30

Acknowledgements

We thank the support of the IHU and the following grants: Institut des neurosciences translationnelle -ANR-10-IAIHU-06, and Infrastructure d'avenir en Biologie Santé -ANR-11-INBS-0006.


Ekanki Sharma

Antonetaiuliana Bratcu

Claudiu Chiculiţă
email: claudiu.chiculita@ugal.ro

Yvon Bésanger
email: yvon.besanger@g2elab.grenoble-inp.fr

Co-simulation of a Low-Voltage Utility Grid Controlled over IEC 61850 protocol

Keywords: Co-simulation, delay, IEC 61850, Manufacturing Message Specification (MMS), TCP/IP, low-voltage (LV) grid I

This paper presents a co-simulation model using MATLAB ® toolboxes to illustrate an interaction between the communication system and the energy grid, coherent with the concept of smart grid that employs IEC 61850 communication standard. The MMS (Manufacturing Message Specification) protocol supported by IEC 61850, based on TCP/IP, is used for the vertical communication between the Supervisory and Data Acquisition (SCADA) system and Intelligent Electronic Devices (IEDs) embedding the local control of different parts of the smart grid. In this paper an IED supporting the power control of a photovoltaic (PV) plant connected to a low-voltage (LV) utility grid is considered. Communication system consisting of the transport layer and a router placed on the network layer is modeled as an event-driven system using SimEvents ® toolbox and energy grid is modeled as a time-driven system using SimPowerSystems ® toolbox. Co-simulation results are obtained by combining different communication scenarios and time-varying irradiance scenarios for the PV plant when this one is required to provide a certain power in response to a power reference received from SCADA over the communication network. The analysis aims at illustrating the impact that stochastic behavior and delays due to network communication have on the global system behavior.

INTRODUCTION

The concept of smart grid at distribution level characterizes the information flow using latest communication technologies. A smart distribution grid is a modernized electrical grid (fully automated distribution system) that uses information and communication technology to improve reliability, security and efficiency of energy delivery [START_REF] Li | Cosimulation for Smart Grid Communications[END_REF]. So smart grids use complex commu-nication networks in order to allow communication between the DERs (Distribution Energy Resources) to implement functions such as power monitoring, energy management and ADA functions (Advanced Distribution Automation). Smart grids are therefore highly automated cyber-physical systems. The concept of smart distribution grid has led to the development and implementation of new communication standards able to provide interoperability [START_REF] Padilla | Towards Smart Integration of Distributed Energy Resources Using Distributed Network Protocol Over Ethernet[END_REF][START_REF] Manbachi | Real-Time Adaptive VVO/CVR Topology Using Multi-Agent System and IEC 61850 -Based Communication Protocol[END_REF]. This increasing role of communication network in the operation and control of power system results in strong dependence of the global performance on the performance of communication [START_REF] Li | Cosimulation for Smart Grid Communications[END_REF]. As information flows through the communication network, it exhibits variable delays that can cause unstable behavior and some asynchronous events that can lead to loss of data. This highly automated power system is prone to vulnerabilities caused by interdependencies with communication and information systems [START_REF] Rozel | Towards a Common Model for Studying Critical Infrastructure Interdependencies[END_REF].

Recent work has been reported in the literature [START_REF] Rozel | Towards a Common Model for Studying Critical Infrastructure Interdependencies[END_REF][START_REF] Caire | Vulnerability Analysis of Coupled Heterogeneous Critical Infrastructures: a Co -simulation Approach with a Testbed Validation[END_REF] on identifying and modeling these interdependencies. The implementation of widely used advanced protocols like IEC 61850 has been investigated in [START_REF] Manbachi | Real-Time Adaptive VVO/CVR Topology Using Multi-Agent System and IEC 61850 -Based Communication Protocol[END_REF] and [START_REF] Ferrari | Co-simulation of Network Infrastructure for Substation Automation Systems[END_REF]. The detailed simulation of energy systems along with the communication architecture requires hybrid models, i.e., models including both continuous time-based and discrete-event-based simulation aspects. There are different domain-specific simulation tools that have been developed and described in the literature to simulate the dynamic system behavior for power systems in different frameworks, e.g., electrical systems, cyber-physical systems, control or stochastic environmental behavior [START_REF] Stifter | Co-simulation of Components, Controls and Power Systems based on Open Source Software[END_REF]. A commonly tool available for such kind of analysis is Simulink ® of MathWorks that can help in studying the behavior of models with hybrid dynamics. Other than this, developing co-simulation framework is also possible that allows the use of most appropriate expert tool for each domain. But in this case it has a disadvantage of making the different tools to work together efficiently [START_REF] Stifter | Co-simulation of Components, Controls and Power Systems based on Open Source Software[END_REF]; e.g., time synchronization and data exchange of two fundamentally different simulators are major challenges [START_REF] Stevic | A Two-Step Simulation Approach for Joint Analysis of Power Systems and Communication Infrastructures[END_REF]. Concerning the co-simulation framework, there are lot of approaches proposed in the literature that include using free and open-source available tools such as GridLAB-D and PSAT [START_REF] Stifter | Co-simulation of Components, Controls and Power Systems based on Open Source Software[END_REF], or co-simulation tools, such as VPNET [START_REF] Stevic | A Two-Step Simulation Approach for Joint Analysis of Power Systems and Communication Infrastructures[END_REF], combining INSPIRE & HLA [START_REF] Georg | Analyzing Cyber-Physical Energy Systems: The INSPIRE Cosimulation of Power and ICT Systems using HLA[END_REF]. The work proposed in this paper is using a more flexible approach for setting up a kind of co-simulation scenario in MATLAB ® , as compared to the one discussed in [START_REF] Stifter | Co-simulation of Components, Controls and Power Systems based on Open Source Software[END_REF][START_REF] Stevic | A Two-Step Simulation Approach for Joint Analysis of Power Systems and Communication Infrastructures[END_REF][START_REF] Georg | Analyzing Cyber-Physical Energy Systems: The INSPIRE Cosimulation of Power and ICT Systems using HLA[END_REF][START_REF] Jablkowski | A Hardware-in-the-Loop Co-Simulation Architecture for Power System Applications in Virtual Execution Environments[END_REF][START_REF] Tipper | A Microgrid Co-Simulation Framework[END_REF]. Fig. 1 suggests the main idea of the implemented co-simulation architecture.

Use of IEC 61850 protocol over the energy grid results in hybrid global dynamics. This paper deals with the implementation of Manufacturing Message Specification (MMS) protocol supported by IEC 61850 communication protocol standard for controlling a photovoltaic (PV) plant connected to a low-voltage (LV) utility grid. Communication system consisting of the transport layer and a router placed on the network layer is modeled as an event-driven system using SimEvents ® toolbox and energy grid is modeled as a time-driven system using SimPo-werSystems ® toolbox. Mechanisms related to acknowledgments and retransmissions specific to TCP/IP as well as management of router's queue length by means of Random Early Detection (RED) algorithm and congestion avoidance are implemented.

Co-simulation is performed to illustrate interaction between the communication system and the power system in terms of latencies affecting the global behavior of PV system. Fig. 1 depicts the general scenario of the studied interaction. This paper is organized as follows. Section II discusses the IEC 61850-based communication system modelling. It is followed by Section III that describes LV utility grid modelling. In Sections IV and V, the coupled model and co-simulation results are presented, respectively, followed by conclusion of the work in Section VI, the final one.

II. IEC 61850-BASED COMMUNICATION SYSTEM

MODELLING USING SIMEVENTS ® / MATLAB ®

A. IEC 61850-based Communication System

IEC 61850 standard characterizes a comprehensive structure for designing power utility automation system. This standard provides a data flow mechanism to exchange information among the Intelligent Electronic Devices (IEDs) that embed the local control of different parts of smart grid. Fig. 2 shows the communication architecture of the standard, including three protocols: GOOSE (Generic Object Oriented Substation Event) messaging, MMS (Manufacturing Message Specification) and SV (Sampled Values) with two interface buses that support the data exchange.

B. Relation between IEC 61850 Standard and TCP/IP

IEC 61850 standard supports three protocols, namely GOOSE, SV and MMS to allow the data flow. MMS protocol supports MMS traffic that allows SCADA systems placed at supervisory level to access vertically all the IED objects. Its functioning is based on TCP/IP (layer 3). GOOSE traffic allows IEDs to access data horizontally. SV traffic carries the current and voltage samples used in discrete-time control. Out of the three mentioned protocols, MMS is the one that is based on TCP/IP layers and allows traditional supervision operations, i.e., read, write and report.

C. Transport Layer in TCP/IP

Transport layer is an end-to-end protocol that provides communication between two processes (application programs) running over different hosts. This layer sends data in segments. Transport layer supports two protocols: UDP (User Datagram Protocol) and TCP (Transmission Control Protocol). UDP is an unreliable connectionless protocol, whereas TCP is a connection-oriented protocol that provides reliable transport by establishing the connection between client and server before data transfer begins and releasing the connection after data transfer ends up. Hence TCP ensures error and flow control, as well as congestion control [START_REF] Tanenbaum | Computer Networks[END_REF][START_REF] Forouzan | Data Communications and Networking[END_REF].

D. TCP with RED (Random Early Detection) Router

It is difficult to build a router that does not drop packets even when it is overloaded with the incoming flow of packets. In order to allow the router to work properly, RED [START_REF] Floyd | Random Early Detection Gateways for Congestion Avoidance[END_REF] algorithm is implemented, so as to keep a check on the moment of time the router should start discarding packets. Algorithm calculates the average queue length and, if it exceeds a pre-specified threshold, some packets are dropped on random basis. Main idea of the generalized RED algorithm is given below:

for each packet arrival calculate the average queue size (avg) if avg calculate probability with probability : mark the arriving packet else if mark the arriving packet E. Error Control in TCP TCP provides reliability using error control based on mechanism for detecting corrupted, lost, out-of-order and duplicated segments using simple tools like acknowledgment, retransmission and segment reassembly.

• Acknowledgment: TCP uses ACK signals to confirm the receipt of data segments.

• Retransmission: sender retransmits the data segment, in case segment is lost or delayed. This occurs in two conditions: when retransmission timer expires or when sender receives three duplicate ACKs.

Retransmission RTO: one retransmission time-out timer is maintained for all sent, but not acknowledged segments. When timer matures, earliest outstanding segment is retransmitted. The value of RTO is updated based on round trip time (RTT) of segments.

RTT: time needed for a segment to reach the destination and for an acknowledgement to be received.

• Segment reassembly: TCP guarantees data to be delivered to process in ordered sequence.

F. TCP Congestion Control

Congestion control refers to a mechanism where the load is kept below the capacity. Unlike UDP, TCP is able to control the network congestion by using congestion control algorithm that has four parts: Slow Start (SS), Congestion Avoidance (CA), Fast Retransmit and Recovery (FR). The algorithm for slow start and congestion avoidance is described below. threshold, slow start stops and CA starts operating, which causes the additive increase in cwnd. When congestion is detected, FR starts working and sender retransmits the segment. Then immediately fast recovery starts operating, which allows CA algorithm to take over the operation instead of allowing slow start. This operating mode is shown in Fig. 4.

G. TCP/IP-based Communication System Modelling in MATLAB ® /SimEvents ®

SimEvents ® is a model-based simulation MATLAB ® toolbox and is useful to model non-deterministic discrete-event systems, whose states changes upon the occurrence of discrete events. The toolbox is installed as a part of Simulink ® library for MATLAB ® and appears as a block library [START_REF] Xu | Networked Modeling and Simulation Based on SimEvents[END_REF]. Discrete-event systems can be constructed as entity processor units. Various blocks such as generators, queues and servers can be used to process the entities. This part of work includes an efficient implementation of TCP/IP-based communication system in SimEvents ® . The TCP/IP protocol has been analyzed and its above described characteristic strategies have been implemented using SimEvents ® . Basic simulation models reported in [START_REF] Chiculiţă | Contributions on control structures for telepresence systems[END_REF] have been used here.

III. LOW-VOLTAGE (LV) UTILITY GRID MODELLING USING

SIMPOWERSYSTEMS ® /MATLAB ®

A. Photovoltaic (PV) Generation System

In this paper an IED supporting the power control of a PV plant connected to an LV utility grid is considered. A PV system is a power source designed to supply usable solar power by means of photovoltaic cells. Several single PV cells can be connected in series and parallel to form PV modules, and PV modules (or panels) can themselves be connected in series and pa-rallel in order to increase power.

A PV system includes then an arrangement of solar panels that absorb and convert sunlight into electricity. The PV panels produce DC current and are connected to the AC grid by means of inverters; here, a so-called back-to-back structure consisting of a DC/DC boost converter and DC/AC three-phase inverter is used (Fig. 5).

The PV plant consists of an n-by-m array of PV modules, where n is the number of parallel-connected strings, each with m number of series-connected cells. Here n=7 and m=5. The control block diagram of the grid-connected PV plant in Fig. 5 shows the two levels of control structure: an upper-level control containing the active power control and a low-level control in charge with DC/DC converter's voltage control. Notations in Fig. 5

B. Maximum Power Point Tracking (MPPT)

Grid-connected PV systems are usually required to provide a certain level of power whenever possible under the given irradiance conditions. When irradiance level is not sufficient to ensure the required power level, PV plant is controlled to provide the maximum power available. It is thus important to have an effective and appropriate MPPT algorithm for the PV system irrespective of irradiance variation. Lot of MPPT algorithms have been implemented in the literature [START_REF] Yu | A Study on Maximum Power Point Tracking Algorithms for Photovoltaic Systems[END_REF][START_REF] Ishaque | A review of maximum power point tracking techniques of PV system for uniform insolation and partial shading condition[END_REF][START_REF] Bhatnagar | Maximum power point tracking control techniques: State-of-the-art in photovoltaic applications[END_REF][START_REF] Salas | Review of the maximum power point tracking algorithms for stand-alone photovoltaic systems[END_REF][START_REF] Gomes De Brito | Evaluation of the Main MPPT Techniques for Photovoltaic Applications[END_REF], that include, among others, Perturb & Observe (P&O), Incremental Conductance (IC), Hill Climbing (HC).

In this work Hill Climbing algorithm has been used, namely by varying the DC/DC boost converter duty cycle (D) by comparing the power at present and previous time. Thus, if incremental power 0, then D is increased in order to make 0; otherwise, if 0, then D is reduced to make 0. The current operating point is moved towards the maximum power one, characterized by ⁄ 0.

C. Modified Maximum Power Point Tracking

The MPPT algorithm has been modified with the possibility of power control in order to get a certain power level, denoted by . If the required power is larger than the maximum one that PV plant can deliver, then MPPT operation must be ensured. In Fig. 6 the general power-voltage curve of a PV plant for a given irradiance level is shown.

is the power reference coming from an upper level; here the upper level is denoted as a dispatcher or supervisor, where the dispatcher is unaware of the maximum power that PV plant can provide at a given time. Modifications brought to the original MPPT algorithm are explained next. P(n) and P(n-1) denote measured PV power at current time and previous time steps, respectively. .

So in case B sign of voltage variation must be reversed compared to the one for MPPT.

In conclusion, wherever the current operating point is placed, if voltage variation is the one used for MPPT multiplied by , then power regulation at is whenever irradiance level allows it. The modified MPPT algorithm flowchart is given in Fig. 7.

IV. CO-SIMULATION MODEL

It is considered that the PV plant included in the low-voltage distribution grid receives power setpoint, , from a dispatcher which belongs to a more complex and remotely located SCADA system. Therefore, reference signal must be transmitted over a communication network. The protocol used to this end is the MMS from IEC 61850 operating over TCP/IP, used for communication between SCADA and the considered low-voltage distribution grid control system acting as an IED. Fig. 8 shows the coupling between the two models as a single Simulink ® diagram, i.e., the event-driven communication model -detailed in Fig. 9 and the time-driven model of the low-voltage distribution grid together with its control. An event-to-timed based signal conversion block was necessary to couple these two models.

V. CO-SIMULATION RESULTS

Simulation results shown next have been obtained under nominal conditions for the communication system. Thus, the network layer in the communication model is equipped with buffer space (queue) at the router that can store maximum 20 packets at a time. The router buffer has high capacity link at its input (as it can receive large number of packets without delay), but limited channel capacity output of 30000 bytes/second and this can lead to congestion in the network. The communication model has a constant line propagation delay of 0.1 seconds. This part of work is not concerned with the router behavior of distributing packets along the different paths, but instead with the router behavior in case of congestion, as that affects the delay times of the packets travelling through it. The delay experienced on the end-to-end path in the network can be decomposed into a deterministic (e.g., propagation) and a stochastic (e.g., queuing) part. The deterministic delay is easy to estimate depends on factors such as packetization delay, p and service delay. The queuing delay is difficult depends on the congestion in the network. In check on the queuing delay, the network load control. In the upcoming simulation results we ex of increasing load in the network on the delay ex packets travelling through the router. One included in the communication model, i.e., th equal to packet length divided by the channel cap calculated. In real world, power reference is a co could happen that it varies step-wise, so we step-wise power reference signal to the system. This input is first fed to the communication received power reference is fed to the power sy irradiance profile (ramp-up down) with peak W/m 2 , results in Fig. 10 show that power injecte following the profile of input irradiance. The po set up as 8 kW. The first two figures show the e ified MPPT algorithm: as long as the irradian cient for the PV plant to provide the required p at the utility grid is around the required one, i. for the losses in the converters and transform diance level decreases so that has reduce erence , this means that the PV plant is no One more scenario has been examined that of adding a perturbation in the system. A pertu introduced as an extra data flow. After adding the load on the network layer will increase an queue more burdened with the incoming da will drop more packets as compared to the p as to manage the dropping of packets, if space will be increased, then this will cause delay t or evaluate, as it propagation delay t to estimate as it order to keep a d must be under xamine the effect xperienced by the more parameter e service timepacity -has been onstant, at most it have provided a n model, and the stem. Setting the k value of 1000 ed into the grid is ower reference is effect of the modce level is suffipower, the power e., 8 kW (except mers). If the irraed below the refot able to provide the required power any longer. At th starts operating.

Queuing delay depends on the tr hence varies stochastically over the point-to-point delay experienced by the whole system, start and read ti duced as they jointly perform the timers compute the time each entity Emitter Transport level and departi port level. Fig. 10 is reporting the power reference and the received the control system of the grid-conn minal condition (router buffer: 20 channel capacity: 30000 bytes/secon the received reference power sent f linear increase due to the fact tha packets and as they are moving to router, they are adding up the delay this delay, power reference receiv system is delayed. The last non de reference -due to use of a zero-ord superposed plots show. This event i the communication between the S plant control.

ion model coupled with the model of a 10-kW grid connected PV array. t shows the effect urbation has been g this perturbation nd this will make ata. Hence queue revious case. So, e in router buffer to increase more.

So a better solution is to increase th case it has been increased to 60000 cause the service time for the pa previous case. Hence, channel woul more quickly and network congesti can be validated from the result in delay between the sent power refere reference.

hat moment MPPT algorithm raffic load along the path and e time. So as to measure this y the data traversing through imer blocks have been introe timing computation. These y takes between arriving at the ing from the Receiver Transcomparison between the sent power reference which feds nected PV array under the no-0 packets, line delay: 0.1 s, nd). The delay experienced in further to the system shows a at queue has become full of wards the head of the buffer y (queuing delay). Because of ved by the PV plant control elayed sample value is fed as der hold element -as the two is a short-time perturbation in SCADA system and the PV he channel capacity (in this 0 bytes/second), as this will ackets to be less than the ld be able to output the data ion would be avoided. This Fig. 11 that shows shorter ence and the received power

VI. CONCLUSION

This paper investigates how the interaction between the communication system and the energy grid affects the global behavior of a smart grid, as the communication parameters affect the data transmission in the power system. The communication model based on TCP/IP protocol (MMS from IEC 61850) has been incorporated for establishing the vertical communication between the SCADA system and IEDs, where the IED considered here is the local control system of a PV array connected to an LV grid. The MATLAB ® co-simulation results have shown changing the parameters such as input to the communication system, router buffer, and channel capacity changes the data flow in the power system. It can be concluded that if the input data flow is increased, the router buffer can help in reducing the data loss but on the other hand it will lead to increase in delay. A better idea is to increase the channel capacity wherever possible, so as to avoid the network congestion.

Future work will aim at enriching the proposed co-simulation model and using it to derive an as general as possible formal model suitable for a large class of smart grids.

Fig 1 .

 1 Fig 1. Overview of co-simulation approach implemented.

Fig 2 .

 2 Fig 2. IEC 61850 communication architecture [12].

Fig 3 .

 3 Fig 3. OSI mapping of IEC 61850 protocol.

  In each RTT cwnd = 2*cwnd Congestion avoidance: Starts cwnd ssthresh New ack received: cwnd = cwnd + 1/cwnd In each RTT cwnd = cwnd+1 SS starts by setting congestion window (cwnd) to one MSS (maximum segment size), size of cwnd increases by one MSS each time an acknowledgment is received. This results in cwnd increasing exponentially, because at each step it doubles. The congestion window cwnd increases till it reaches ssthresh (slow start threshold) tracked by the sender. When cwnd reaches the

Fig 4 .

 4 Fig 4. Congestion control algorithm [13].

Fig 5 .

 5 Fig 5. Grid-connected PV energy conversion system along with a two-level control structure[START_REF] Bratcu | Maximum Power Point Tracking of Grid-connected Photovoltaic Arrays by Using Extremum Seeking Control[END_REF].

Fig 6 .

 6 Fig 6. Power-voltage curve of a PV plant for a given irradiance level Irr1.

2 : 1 . 1 .CASE B. 2 :

 2112 Voltage must decrease as long as 1 So in case A sign of voltage variation is same as in MPPT. CASE B : When CASE B.1: Voltage must decrease as long as 1 Voltage must increase as long as 1 1

Fig 7 .

 7 Fig 7. Modified MPPT algorithm.

Fig

  Fig 8. Communicati

Fig 9 .

 9 Fig 9.

FFig 11 .

 11 Fig 11. Co-simulation re

Fig 10 .

 10 Fig 10. Co-simulation results in the nominal case.