Davide Fonti

Stanislas Lagarde

Francesca Pizzo

Wala Aboubakr

Christian G Bénar

Bernard Giusiano

Prof Fabrice Bartolomei
email: fabrice.bartolomei@ap-hm.fr

Christian Benar

Parieto-premotor functional connectivity changes during parietal lobe seizures

Keywords: parietal lobe seizures, motor semiology, Stereoelectroencephalography, functional connectivity. Highlights

Introduction

Parietal lobe epilepsies (PLE) are defined by seizures originating from the posterior parietal cortex. PLE is a rare medical condition representing up to 6% of focal epilepsies in surgical series [START_REF] Asadollahi | Drug-resistant parietal lobe epilepsy: clinical manifestations and surgery outcome[END_REF][START_REF] Bartolomei | Neural networks underlying parietal lobe seizures: A quantified study from intracerebral recordings[END_REF][START_REF] Francione | Drug-resistant parietal epilepsy: Polymorphic ictal semiology does not preclude good post-surgical outcome[END_REF]V. Salanova et al., 1995a). Furthermore the incidence of PLE is probably underestimated due to misleading clinical and electrical features that can make differential diagnosis difficult from epileptic seizures arising from other lobes or even psychogenic nonepileptic seizures [START_REF] Francione | Drug-resistant parietal epilepsy: Polymorphic ictal semiology does not preclude good post-surgical outcome[END_REF][START_REF] Pilipović-Dragović | Long-term seizure outcome following epilepsy surgery in the parietal lobe: a meta-analysis[END_REF]V Salanova et al., 1995;[START_REF] Sveinbjornsdottir | Parietal and Occipital Lobe Epilepsy: A Review[END_REF]. Functionally, the parietal lobe is a multimodal cortex that receives and integrates somato-sensory, visual and auditory inputs and that participates in motor planning, attention, spatial representation, working and long-term memory, calculation and language processes [START_REF] Orban | Mapping the parietal cortex of human and non-human primates[END_REF]. Posterior parietal cortex is strongly interconnected with sub-cortical and others cortical areas [START_REF] Caspers | Microarchitecture and connectivity of the parietal lobe[END_REF]. Parietal areas are strongly connected with frontal areas for multiple aspects of motor planning, in particular posterior parietal lobe with pre-motor areas (supplementary motor area (SMA), pre-SMA) and pre-frontal areas) [START_REF] Andersen | Intention, Action Planning, and Decision Making in Parietal-Frontal Circuits[END_REF][START_REF] Caspers | Microarchitecture and connectivity of the parietal lobe[END_REF][START_REF] Fogassi | Motor functions of the parietal lobe[END_REF][START_REF] Luppino | Corticocortical connections of area F3 (SMA-proper) and area F6 (pre-SMA) in the macaque monkey[END_REF][START_REF] Wise | PREMOTOR AND PARIETAL CORTEX: Corticocortical Connectivity and Combinatorial Computations[END_REF]. The high connectivity of the parietal lobe explains the involvement during PLS seizures of many extra-parietal regions. Some authors have suggested that the variety of clinical signs observed in PLS could be due to the involvement of distinct systems in the propagation of ictal discharge [START_REF] Francione | Drug-resistant parietal epilepsy: Polymorphic ictal semiology does not preclude good post-surgical outcome[END_REF]V. Salanova et al., 1995a).

One of the most striking features of parietal seizures are motor signs, which are reported in 75 to 88% of PLS [START_REF] Asadollahi | Drug-resistant parietal lobe epilepsy: clinical manifestations and surgery outcome[END_REF][START_REF] Bartolomei | Neural networks underlying parietal lobe seizures: A quantified study from intracerebral recordings[END_REF][START_REF] Francione | Drug-resistant parietal epilepsy: Polymorphic ictal semiology does not preclude good post-surgical outcome[END_REF][START_REF] Kim | Parietal Lobe Epilepsy: The Semiology, Yield of Diagnostic Workup, and Surgical Outcome[END_REF]V. Salanova et al., 1995b;[START_REF] Sveinbjornsdottir | Parietal and Occipital Lobe Epilepsy: A Review[END_REF]. The mechanisms underlying the occurrence of motor signs are poorly understood and no study has specifically addressed this question. SEEG during presurgical evaluation for drug resistant epilepsy allows recording of the electrical activity of cerebral areas (including deep structures) and represents a unique way to study and quantify epileptogenic networks in focal epilepsies [START_REF] Bartolomei | Defining epileptogenic networks: Contribution of SEEG and signal analysis[END_REF]. In a previous study, using the Epileptogenicity Index Method (EI), we observed that neural networks underlying PLS involved different subnetworks in the parietal lobe and that their complexity was higher when maximal epileptogenicity was in the precuneus or in Brodman area 7 [START_REF] Bartolomei | Neural networks underlying parietal lobe seizures: A quantified study from intracerebral recordings[END_REF].

SEEG can also study epileptogenic networks through the estimation of functional connectivity (FC). FC allows an estimate of the functional coupling between brain regions in physiological and pathological brain processes [START_REF] Bartolomei | Defining epileptogenic networks: Contribution of SEEG and signal analysis[END_REF][START_REF] Fornito | Fundamentals of Brain Network Analysis[END_REF][START_REF] Guye | Imaging structural and functional connectivity: towards a unified definition of human brain organization?[END_REF]. For example, in epilepsy FC is altered during both interictal and ictal periods [START_REF] Bartolomei | Defining epileptogenic networks: Contribution of SEEG and signal analysis[END_REF] for review; ). In the latter case, it provides insights into the modifications of brain networks during seizures (desynchronisation or hypersynchronisation depending on the features of the ictal discharge) and how these changes are involved in clinical manifestations [START_REF] Bartolomei | The role of the dorsal anterior insula in ecstatic sensation revealed by direct electrical brain stimulation[END_REF][START_REF] Bartolomei | Defining epileptogenic networks: Contribution of SEEG and signal analysis[END_REF][START_REF] Bartolomei | How does vagal nerve stimulation (VNS) change EEG brain functional connectivity?[END_REF][START_REF] Hagiwara | How can we explain the frontal presentation of insular lobe epilepsy? The impact of non-linear analysis of insular seizures[END_REF].

In this study, we analyzed seizures recorded by SEEG originating from the posterior parietal cortex. Our main objective was to evaluate the changes in FC occurring at seizure onset and during seizure propagation and to correlate these changes to the clinical semiology. In particular we focused our work on motor semiology, which is one of the most frequent presentations of parietal seizures. For this purpose we compared the FC changes during PLS with and without motor semiology. We hypothesized that FC between parietal and extraparietal frontal areas would be altered differently according to the generation or absence of ictal motor symptoms. Moreover, the extension of epileptogenicity in extra-parietal regions can be the cause of some surgical failures [START_REF] Ristić | Parietal lobe epilepsy: The great imitator among focal epilepsies[END_REF]. Therefore, connectivity changes could be a reflection of the extension of epileptogenic networks. We also evaluated FC changes and surgical outcome as a secondary objective of our study.

Methods

Patients

We method as previously defined [START_REF] Bartolomei | Epileptogenicity of brain structures in human temporal lobe epilepsy: A quantified study from intracerebral EEG[END_REF]. The EI is a semi-automatic method to quantify the dynamic of the fast-activity genesis at the beginning of a seizure. The EI is calculated from two important features of the transition from pre-ictal to ictal activity: (i) the redistribution of signal frequencies from lower frequency bands (theta, alpha) toward higher frequency bands (beta, gamma); and (ii) the time delay over which appearance of these highfrequencies occur at the seizure onset. EI values are normalized at the patients level by the maximum value, and then ranges from 0 (no epileptogenicity) to 1 (maximal epileptogenicity).

A threshold > 0.3 was chosen as an indicator of regions included in the SOZ. Indication for video-SEEG and the anatomical targeting of electrodes was established in each patient on the basis of non-invasive information and hypotheses about the SOZ. SEEG explorations were done using intracerebral multiple contacts electrodes (5-18 contacts, length: 2 mm, 1,5 mm apart, diameter 0,8 mm, placed according Tailairach stereotactic method or with frameless method [START_REF] Skoch | Subdural grid and depth electrode monitoring in pediatric patients[END_REF][START_REF] Talairach | New approach to the neurosurgery of epilepsy. Stereotaxic methodology and therapeutic results. 1. Introduction and history[END_REF]. A post implantation computerized tomography was done to exclude complications and for localisation of each electrode contact by using Gardel GUI (https://meg.univamu.fr/wiki/GARDEL:presentation) for fusioning with pre-implantation magnetic resonance imaging [START_REF] Medina Villalon | A software suite for presurgical brain mapping in epilepsy: Intracerebral EEG[END_REF]. SEEG signals were recorded with 256 channel

Deltamed™, sampled at 256 Hz or 512 Hz and recorded on a hard disk ( 16 bit/ samples) without using a digital filter. A high-pass filter was used (cut-off 0.16 Hz at -3 db) to remove very slow variations that sometimes contaminate the baseline.

Clinical data

For each patient included in the study one representative spontaneous seizure recorded during video-SEEG monitoring was selected. The semiology of each seizure was analysed, with particular attention to the presence of limb clonic, tonic/ dystonic phenomena or hypermotor semiology. According to this analysis two groups of seizures were determined : a "motor seizures" group , in which was present at least one of the aforementioned signs, and a "nonmotor seizures" group, in which seizures occurred without motor signs.

Functional Connectivity (FC) analysis

FC analysis was done using the software (https://meg.univ-amu.fr/wiki/AnyWave) [START_REF] Colombet | AnyWave: A cross-platform and modular software for visualizing and processing electrophysiological signals[END_REF]. In each SEEG seizure, we selected three temporal periods for FC analysis:

Background (BG) representing pre-ictal activity, which could terminate up to one second before ictal discharges; Seizure onset (SO) representing seizure onset period and commencing at the beginning of the ictal discharges, Seizure Propagation period (SP) which commenced at the end of SO, when ictal discharges propagated and became more rhythmic and of higher voltages . Each time period was of the same duration in each patient, averaging 16 s across the cohort (Standard deviation (SD) 11), depending on the total duration of each analysed seizure.

We selected eight areas sampled in all patients. 4 parietal areas: 1) precuneus, 2) superior parietal lobule (SPL), 3) inferior parietal lobule (IPL) or posterior cingular gyrus (pCG), 4) parietal operculum (OP) or Post-central Gyrus (PoCg); 2 frontal areas ( supplementary motor area (SMA) and lateral pre-motor area (pmL); and 2 temporal areas (mesial and lateral). For this purpose, we co-registered pre-implantation MRI and post-implantation CT [START_REF] Medina Villalon | A software suite for presurgical brain mapping in epilepsy: Intracerebral EEG[END_REF] and visually defined the anatomical areas. We chose one bipolar channel (from two adjacent contacts) in grey matter for each specific brain area. As implantation schemes vary largely according to SOZ localisation hypothesis, not all regions were explored in all patients selected. To address the changes in functional connectivity occuring between parietal/temporal/frontal regions we studied 16 specific links as indicated in Supplementary Figure 1.

Interdependencies between bipolar SEEG signals were estimated by pairwise nonlinear regression analysis based on the h² coefficient [START_REF] Wendling | Interpretation of interdependencies in epileptic signals using a macroscopic physiological model of the EEG[END_REF], computed using the Anywave software [START_REF] Colombet | AnyWave: A cross-platform and modular software for visualizing and processing electrophysiological signals[END_REF]. Nonlinear regression analysis is aimed at estimating the degree of association between two signals X and Y independently from the linear or nonlinear nature of this association. The basic idea is to describe the amplitude of signal Y as a function of the amplitude of signal X using a nonlinear regression curve and to compute the variance of Y that is explained, or predicted, by X according to this regression curve. This technique has been shown to be particularly suitable for the analysis of EEG/SEEG signals in the context of epilepsy (review in [START_REF] Bartolomei | Defining epileptogenic networks: Contribution of SEEG and signal analysis[END_REF] and achieves high accuracy in the setting of intracranial EEG recordings [START_REF] Ansari-Asl | Quantitative evaluation of linear and nonlinear methods characterizing interdependencies between brain signals[END_REF][START_REF] Wang | A systematic framework for functional connectivity measures[END_REF][START_REF] Wendling | From EEG signals to brain connectivity: A modelbased evaluation of interdependence measures[END_REF]. In a sliding window, a piecewise linear regression is performed between each pair of signals. The h 2 is the coefficient of determination, which measures the goodness of fit of the regression (equivalent to the r 2 used in linear regression). The h 2 is bounded between 0 (no correlation) and 1 (maximal correlation) and is asymmetric. We used a sliding window of 3 seconds with an overlap of 2 seconds and a maximum delay between signals of 100 ms. We used 0.5 Hz high pass filter in all analyses. Because h 2 values are asymmetric, we symmetrised our connectivity matrix using the higher h 2 value between two channels (chan1-> chan2 and chan2->chan1). We thus obtained connectivity graphs, with each bipolar channel representing a node of the graph, and h 2 values between two brain areas the edges. For each seizure, FC graphs of the different periods were compared (SO versus BG, SP versus BG, SP versus SO). For studying the FC changes between parietal, frontal and temporal regions, for each patient, a z-score values of 16 cortico-cortical interaction comparing the different above mentioned periods were obtained . We also regrouped interactions between parietal and temporal or frontal regions.

Statistical analysis

We looked at how the changes in FC may differ in seizures with and without motor manifestations. Firstly, we computed the z-scores of FC changes (measure of the difference between periods of seizure: BG, SO and SP) in parieto-extra-parietal interactions. The zscores, obtained comparing FC graphs, were computed with the difference of h² values derived from two different windows (SO versus BG, SP versus BG, and SP versus SO) and between two specific cortical areas, normalized with SD. In this way z-scores represent a normalised index of FC changes for a given cortico-cortical interaction between different periods.

Finally, we analysed the relation between ictal FC changes and post-surgical seizure outcome.

For this, we used the mean difference and permutation t-test for assessing statistical differences between motor and non-motor seizures and between patients with a positive surgical outcome (Engel's classes I and II) or an unfavourable surgical outcome (Engel's classes III and IV).

We also used χ² test, Fisher's exact test, Student's t test, Spearman'R for the other analyses.

Statistical significance was established for p≤0.05. For Statistical analysis we utilised IBM ®SPSS ®STATISTICS software, 23th edition and Estimation statistics beta (https://www.estimationstats.com/#/) [START_REF] Ho | Moving beyond P values: data analysis with estimation graphics[END_REF].

Results

We studied 22 patients with a SOZ localized in the parietal cortex (see Table 1 In 14 patients who underwent a surgical intervention, outcome (with a mean follow of 7.9 years, range 9 months-15 years) was favourable (Engel class I and II) in 8 patients (57.1%) (Table 1). SEEG was performed with a mean number of 9.64 electrodes per patient (min 6, max 18).

Semiology of motor and non-motor seizures

We analysed 22 SEEG recorded seizures, one for each patient. 13 seizures (59%) were classified as "motor seizure", 9 (41%) as "non-motor seizure". Age at SEEG, was younger (20.9 years +-12.0) in the "motor-seizure" group than in "non-motor seizure" group, (30.8 years +-6.2), (t test for independent samples, p= 0.035). Age at epilepsy onset and epilepsy duration

were not significatively different in the two groups. No significative correlation was found between the two groups of seizures for aetiology and surgical outcome. The clinical features are summarized in Table 1.

Functional connectivity (FC) analysis

The figure 1 illustrates an example of PLS with motor semiology in patient 7.

At the group-level, during seizure onset (SO), most of parietal-extra-parietal links showed a decrease in FC. When looking at SO vs BG comparison, significant lower z-scores values were found in motor seizures (mean difference in z-scores= -0.918 (95% confidence interval (CI) -1.5, -0.368); p= 0.0018; Figure 2A), in particular between OP/PoCg and mesial temporal areas (mean difference in z-scores= -2.52 (95% IC= -3.31,-1.52); p< 0.001; Figure 2B). No other significant differences were found in the studied parietal-extra-parietal interactions.

In contrast, during seizure propagation (SP), most of parietal-extra-parietal links exhibited a tendency to an increase in functional connectivity. Overall, no significant difference was found in parietal-extra-parietal z-scores computed between SP and BG (mean difference in z-scores 0.53 (95% IC -0.179, 1.28); p= 0.149). However, z-score computed between SP and BG windows were significantly higher in motor seizures only for SPL-pmL interaction (mean difference in z-scores 2.5 (95% IC 0.297, 4.46; p=0.037). We also observed that "motor seizures" had a significantly higher FC increase in comparison with "non-motor seizure"

between SP and SO (mean difference in z-scores 1.63 (95% IC 1.01, 2.27); p< 0.001) (see Figure 3A). We observed in seizures with motor signs a significantly higher increase in FC during SP in comparison with seizures without motor signs in some specific links (see Figure 3B):

• Precuneus -Premotor Lateral (mean difference in z-score = 3.69 (95%IC 1.45, 5.95); p=0.0132);

• Superior parietal lobule -Premotor lateral (mean difference in z-score = 3.03 (95%IC 0.754, 5.44); p=0.0344);

• Inferior parietal lobule -SMA (mean difference in z-score = 1.97 (95%IC 0.699, 3.86); p=0.046.

The difference for other links were not statistically significant (Supplementary Table 2). As "motor" and "non-motor" groups differ for age at SEEG, we tested if there was a correlation between age at SEEG and FC z-scores however no significant correlation was found (Spearman R test , p= 0.47 for SO versus BG z-scores, 0.07 for SP versus BG z-scores and 0.98 for SP versus SO z-scores).

Relation with post-surgical outcome

We 4). We did not find significant difference in FC changes at the level of links between each brain area, between patients with favourable versus those with unfavourable post-surgical outcome.

4.Discussion

The main purpose of this work was to study ictal changes of FC in patients with PLE and their relationship with ictal motor signs. Motor signs, including limb tonic, dystonic or clonic contractions or hyperkinetic manifestations, are the most frequent clinical signs of PLS, reported as habitual semiological features in 19 patients (86.4%) in our cohort, this is in agreement with 75-88% reported in previous studies [START_REF] Asadollahi | Drug-resistant parietal lobe epilepsy: clinical manifestations and surgery outcome[END_REF][START_REF] Bartolomei | Neural networks underlying parietal lobe seizures: A quantified study from intracerebral recordings[END_REF][START_REF] Francione | Drug-resistant parietal epilepsy: Polymorphic ictal semiology does not preclude good post-surgical outcome[END_REF][START_REF] Kim | Parietal Lobe Epilepsy: The Semiology, Yield of Diagnostic Workup, and Surgical Outcome[END_REF][START_REF] Rushworth | The left parietal and premotor cortices: Motor attention and selection[END_REF][START_REF] Sveinbjornsdottir | Parietal and Occipital Lobe Epilepsy: A Review[END_REF]. Our main result was that ictal FC is differently modified according to the presence (or the absence) of motor semiology. Overall, FC between parietal, frontal and temporal areas decreased at SO and increased during SP. During the period of SO we observed a more pronounced decrease of FC between parietal and extra-parietal areas in "motor" seizures, particular in the interaction between the most anterior part of parietal lobe (OP/PoCg) and mesial temporal areas. In contrast, during the phase of SP, we observed an increased FC between parietal and premotor areas in "motor" parietal seizures. This pattern of FC changes is in agreement with previous studies, showing in particular that fast ictal onset discharges are associated with decreased synchronization (for review see [START_REF] Bartolomei | Defining epileptogenic networks: Contribution of SEEG and signal analysis[END_REF]. During seizures with motor manifestations, we observed a higher re-increase of FC during SP period between parietal cortex and premotor cortex, the pmL and SMA, in comparison with seizures without motor signs. This result indicates that preferential routes from parietal cortex to frontal cortices are responsible of the prominent motor presentation during PLS.

By definition, our patients disclosed maximal epileptogenicity in the posterior parietal cortex, in particular, half of included patients had maximal epileptogenicity in the precuneus or in the SPL. These two regions are involved in important intrinsic brain networks such as the default mode and motor networks [START_REF] Andersen | Intention, Action Planning, and Decision Making in Parietal-Frontal Circuits[END_REF][START_REF] Bartolomei | Neural networks underlying parietal lobe seizures: A quantified study from intracerebral recordings[END_REF][START_REF] Cauda | Functional connectivity of the posteromedial cortex[END_REF][START_REF] Cavanna | The precuneus: A review of its functional anatomy and behavioural correlates[END_REF][START_REF] Fogassi | Motor functions of the parietal lobe[END_REF][START_REF] Wise | PREMOTOR AND PARIETAL CORTEX: Corticocortical Connectivity and Combinatorial Computations[END_REF] and are therefore highly connected with several brain areas. The posterior parietal cortex is anatomically and functionally well connected with the pmL and with the SMA [START_REF] Andersen | Intention, Action Planning, and Decision Making in Parietal-Frontal Circuits[END_REF][START_REF] Caspers | Microarchitecture and connectivity of the parietal lobe[END_REF][START_REF] Fogassi | Motor functions of the parietal lobe[END_REF][START_REF] Rushworth | Connection patterns distinguish 3 regions of human parietal cortex[END_REF][START_REF] Rushworth | The left parietal and premotor cortices: Motor attention and selection[END_REF][START_REF] Wise | PREMOTOR AND PARIETAL CORTEX: Corticocortical Connectivity and Combinatorial Computations[END_REF]. Evidence from studies in primates, showed that the pmL is highly connected with different parietal regions and most of its inputs arise from the SPL [START_REF] Wise | PREMOTOR AND PARIETAL CORTEX: Corticocortical Connectivity and Combinatorial Computations[END_REF].

Other studies showed a strong connectivity also between the IPL, particularly its rostral part, and the ventral premotor cortex and caudal ventrolateral prefrontal cortex [START_REF] Caspers | Probabilistic fibre tract analysis of cytoarchitectonically defined human inferior parietal lobule areas reveals similarities to macaques[END_REF][START_REF] Rushworth | Connection patterns distinguish 3 regions of human parietal cortex[END_REF]. As well as the IPL, the rostral part of intraparietal sulcus (IPS) is mainly connected to premotor and prefrontal areas [START_REF] Caspers | Microarchitecture and connectivity of the parietal lobe[END_REF][START_REF] Rushworth | Connection patterns distinguish 3 regions of human parietal cortex[END_REF][START_REF] Uddin | Dissociable connectivity within human angular gyrus and intraparietal sulcus: Evidence from functional and structural connectivity[END_REF]. Concerning the underlying structural connectivity, data from human and primates demonstrate the important anatomical role of the superior longitudinal fasciculus with its subcomponents, which can also overlap with the arcuate fascicle. The first bundle creates a dense link between the medial, dorso-lateral SPL and premotor areas, pre-SMA, SMA and prefrontal cortex for body centred action and motor initiative, the second subcomponent links the IPL, anterior and middle IPS with posterior prefrontal and premotor cortices for perception and awareness and finally the third subcomponent links the posterior IPL with ventral premotor areas with posterior temporal areas, probably involved in action imitation [START_REF] Caspers | Microarchitecture and connectivity of the parietal lobe[END_REF][START_REF] Martino | Analysis of the subcomponents and cortical terminations of the perisylvian superior longitudinal fasciculus: A fiber dissection and DTI tractography study[END_REF][START_REF] Schmahmann | Association fibre pathways of the brain: Parallel observations from diffusion spectrum imaging and autoradiography[END_REF]. This strong underlying connectivity between parietal and premotor cortex may explain the preferential ictal coupling we observed in our study.

Our secondary objective was to analyse the association between ictal FC changes and the postsurgical seizure outcome. We found that patients with unfavourable post-surgical seizure outcome showed a higher FC increase between parietal-extra-parietal areas during SP. This result, although derived from ictal FC analysis, seems to be in accordance with interictal evidence that a higher FC between cortical areas involved during seizure propagation and those not involved, is associated with a poorer surgical prognosis [START_REF] Lagarde | Interictal stereotactic-EEG functional connectivity in refractory focal epilepsies[END_REF]. Other interictal studies suggested that more complex networks are associated with poorer surgical outcome in focal epilepsies [START_REF] Antony | Functional connectivity estimated from intracranial EEG predicts surgical outcome in intractable temporal lobe epilepsy[END_REF][START_REF] Bartolomei | Epileptogenicity of brain structures in human temporal lobe epilepsy: A quantified study from intracerebral EEG[END_REF][START_REF] Englot | Global and regional functional connectivity maps of neural oscillations in focal epilepsy[END_REF][START_REF] Nissen | Identifying the epileptogenic zone in interictal resting-state MEG source-space networks[END_REF][START_REF] Schevon | Cortical abnormalities in epilepsy revealed by local EEG synchrony[END_REF][START_REF] Shah | High interictal connectivity within the resection zone is associated with favorable post-surgical outcomes in focal epilepsy patients[END_REF]. Therefore, it seems that ictal FC could be a marker of the extension of epileptic networks and predict seizure outcome following resective surgery. Other studies are needed to better understand if some specific interactions between parietal and extra-parietal areas are particularly important for surgical prognosis and to establish if ictal and interictal FC analysis can represent a possible prognostic biomarker for outcome in epilepsy surgery.

There are some limitations in our study. SEEG has the advantage of directly studying brain activity in several brain regions. However, this technique is limited by sampling and the variability of implantations from one patient to another. We therefore selected the common regions available in our patients and limited our study to the between these regions. We selected our sample by verifying that the maximum EI value was calculated in contacts positioned within the parietal lobe; this method, although semi-automatic, requires an adequate preselection of the SEEG channels and verification by an experienced user, and performs less well in seizure onset patterns without rapid discharge. Another limitation of our work is the variability of the spontaneous seizures selected for each patient, both in terms of the state in which they could occur (sleep, quiet wakefulness, during the execution of tasks)

and in seizure duration. The state where seizures occur can influence the FC of the pre-ictal window (BG) and therefore potentially the z-score values computed.

In conclusion, we found that motor semiology observed during PLS is associated with an increase of synchronisation between parietal and premotor cortex during seizure propagation. It suggests that seizures follow underlying strong connectivity between these two cortices, reflecting underlying structural connections, and leading to a pathological oversynchronisation responsible for the emergence of motor semiology.

Tables and Figures Legends

  retrospectively selected 22 patients affected by drug-resistant focal epilepsy who underwent SEEG for pre-surgical evaluation between April 2005 and September 2019 at the epileptology departement of the Timone Hospital, Marseille, France. The institutional review board of the French Institute of Health (IRB15226) approved this study and written patient consent was obtained. All patients had a seizure onset zone (SOZ) in the posterior parietal cortex . SOZ was established visually and quantitatively by the epileptogenicity index (EI)

Figure 2 :

 2 Figure 2: Differences in z-scores computed between Seizure onset period (SO) and pre-ictal period (BG) in all parietal-extra-parietal interactions studied, between "motor" and "nonmotor" seizures

Figure 3 :

 3 Figure 3: Differences in z-scores computed between seizure propagation period (SP) and seizure onset period (SO) in all parietal-extra-parietal interactions studied, between "motor" and "non-motor" seizures

Figure 4 :

 4 Figure 4: Differences in z-scores computed between seizure propagation period (SP) and seizure onset period (SO) in all parietal-extra-parietal interactions studied, between patients with favourable surgical outcome (Engel's classes I and II) and unfavourable surgical outcome (Engel's classes III and IV).

  was contra-indicated in 5 cases (bilateral/ multifocal epilepsy, high risk of functional damage) (22.7%) and refused in one case (4.5%). The two main aetiologies were focal cortical dysplasia in 8 patients (36.4%) and a stroke in 7 patients (31.8%)

	). The SOZ was
	in the right hemisphere in 8 patients (36.4%), in the left hemisphere in 9 patients (40.9%) and
	bilateral or multifocal in 5 (22.7%). The mean age at epilepsy onset (SD) was 7.1 (±4.3) years
	(range 3 month-14 years). The mean age at SEEG was 24.9 (SD 11.1) years (range 2.5 to 50
	years). Parietal regions disclosing the maximal epileptogenicity were precuneus in 2 patients
	(9.1%), superior parietal lobule (SPL) in 9 patients (40.9%), post-central gyrus (PoCg) in 2 (9.1
	%), inferior parietal lobule (IPL) in 7 (31.8%), parietal operculum (OP) in one (4.5 %) and
	posterior cingulate gyrus (pCG) in one patient (4.5%). A resective surgery (lesionectomy,

cortectomy, or lobectomy) was performed after SEEG exploration in 14 patients (63.6%). One patient was awaiting surgery, one had a long-term improvement after SEEG thermocoagulation, surgery (see Supplementary table

[START_REF] Uddin | Dissociable connectivity within human angular gyrus and intraparietal sulcus: Evidence from functional and structural connectivity[END_REF]

Table 1 . Comparison of clinical features of the motor and non-motor seizures groups.

 1

	N, number of patients; %, percent; SEEG, Stereoelectroencephalography; SOZ, seizure onset
	zone; SD, standard deviation; NS, not significant; ANET, angiocentric neuroepithelial tumour;
	DNET, dysembryoplastic neuroepithelial tumour; FCD, focal cortical dysplasia; PMG,
	polymicrogyria; EI, epileptogenicity index; SPL, superior parietal lobule; IPL, inferior parietal
	lobule; pCG, posterior cingulate gyrus; OP, parietal operculum; PoCg, post-central gyrus; *: one
	patient has no identified manual dominance because of very young age.
	Figure 1.

Ictal Stereoelectroencephalography (SEEG) trace and functional connectivity (FC) changes graphs of patient 7 with hypermotor seizure.

  Ictal SEEG trace, bipolar 13 channels montage, pink parts represent the three periods of

	superior parietal lobule (SPL(7)), GSM9-GSM10 sample supramarginal gyrus (SMG); PoC11-
	PoC12 sample post-central gyrus (PoCg); B1-B2 sample hippocampus (Hip); B11-B12 sample
	anterior medium temporal gyrus (aMTG); LP'11-LP'12 sample contralateral (left) SPL.
	analysis: BG (background), SO (seizure onset) and SP (seizure propagation). In BG period, note
	the interictal spikes in post-central gyrus (PoCg) and inferior parietal lobule (SMG). In SO, this
	activity is replaced by low voltage rapid discharges that involved both parietal and frontal
	areas. During SP, slower rhythmic activities with higher voltage are seen.
	B: Graphic representation of the FC changes (z-score) between the different periods, only
	significative changes are represented, blue links represents significative FC decreasing, yellow
	ones significative FC increasing. Between SO and BG, inferior parietal lobule (IPL), PoCg and
	supplementary motor area (SMA) show decrease of FC with temporal mesial area and central
	areas, FC increase is seen between pre-frontal area and contralateral superior parietal lobule
	(SPL). During SP, in comparison with BG and SO, an important FC increase is seen between
	parietal and frontal areas.
	PM1-PM2 sample medial prefrontal cortex (mPFC); PM5-PM6 sample dorso lateral pre-frontal
	cortex (DLPFC); SA1-SA2 sample supplementary motor area (SMA); SA8-SA9 sample the dorso-
	lateral pre-motor area (pmL) ; R1-R2 sample medial pre-central cortex (mCen); R10-R11
	sample lateral pre-central cortex (lCen); PA1-PA2 sample precuneus (pCu); PA9-PA10 sample

A:

Acknowledgments

We thank Prof Patrick Chauvel, Prof Martine Gavaret, Prof Maxime Guye, Prof Mathieu Milh, Prof. Agnes Trebuchon, Dr Aileen McGonigal, Dr Francesca Bonini, Dr Sandrine Aubert, Dr Geraldine Daquin, Dr Lisa Vaugier, Dr Anne Lépine, Dr Nathalie Villeneuve (Marseille) for the clinical management of some included patients; Prof Didier Scavarda, Dr Romain Carron, Prof Henry Dufour, Prof Jean-Claude Peragut, and Prof Jean Régis (Marseille) for surgical or SEEG procedures. We thank all the paramedical team of the epileptology unit, who allow optimal recordings and clinical testing of the patients explored by SEEG.

We thank Dr Rinki Singh for proofreading and correction of the manuscript.

Authors have no conflict of interest to disclose

Abbreviations


Laurent Mahéo
email: laurent.maheo@st-cyr.terre-net.defense.gouv.fr

Vincent Grolleau
email: vincent.grolleau@univ-ubs.fr

Gérard Rio
email: gerard.rio@univ-ubs.fr

Damping efficiency of the Tchamwa-Wielgosz explicit dissipative scheme under instantaneous loading conditions

Keywords: Dynamical systems, Structural dynamics, Explicit and dissipative time integration algorithm, Tchamwa-Wielgosz scheme Systèmes dynamiques, Algorithme d'intégration temporelle explicite

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

Version française abrégée

Dans le cadre de problèmes de dynamique, l'utilisation des Méthodes des Eléments Finis et des Différences Finies pour les discrétisations spatiale et temporelle fait apparaître des erreurs purement numériques. Parmi les nombreuses méthodes utilisées pour amortir ces oscillations parasites hautes fréquences, on se propose d'étudier et d'analyser l'efficacité amortissante du schéma explicite et dissipatif de Tchamwa-Wielgosz (Eqs. ( 1) et ( 4)) sous des conditions de chargements impulsionnels. La comparaison théorique avec le schéma des Différences Centrées montre le décentrage à droite de ce schéma dissipatif (Eq. ( 8)) pour lequel l'amortissement peut être interprêté comme une perturbation du déplacement (Eq. ( 6)). On montre alors la sensibilité de l'amortissement au ratio du pas de temps de calcul sur le pas de temps critique r t = t t c . Des études 1-D et 3-D montrent ensuite la sensibilité de son amortissement à la discrétisation spatiale utilisée. En fixant les valeurs du ratio r t et du paramètre d'amortissement ϕ, on note tout d'abord l'intérêt d'exploiter les résultats de pertes d'énergie en fonction du nombre de pas de temps pour lequel l'algorithme est appliqué. On remarque alors la quasi insensibilité de l'amortissement au raffinement du maillage pourvu que les dimensions des éléments soient le plus homogène possible. Dans le cas de maillages composés d'éléments de différentes longueurs, l'efficacité amortissante diminue. Même si le chargement est impulsionnel, sa propagation à travers le maillage perd son caractère instantané puisque plusieurs pas de temps sont nécessaires pour « traverser » les plus gros éléments du maillage. En effet, pour chaque élément i, le ratio r i t = t t i c diminue ce qui provoque une perte d'efficacité amortissante. Dans le cas 3-D, la différenciation entre les propagations longitudinales faisant apparaître un ratio r L t = t t L c et transversales (r T t ) qui sont différents confirme ces résultats.

Introduction

In the field of structural dynamics, the space and time discretization procedures used invariably introduce numerical errors. Damping methods are therefore often used to limit the effects of these spurious oscillations. The Bulk Viscosity method [START_REF] Von Neuman | A method for numerical calculation of hydrodynamic shocks[END_REF][START_REF] Landshoff | A numerical method for treating fluid flow in the presence of shocks[END_REF], for instance, consists in adding an artificial viscosity to the material behavior. This method has been widely used in several industrial codes such as Abaqus [START_REF]ABAQUS Theory Manual[END_REF] or LS-Dyna [START_REF] Hallquist | LS-Dyna 3D theoretical manual[END_REF]. Many implicit time integration methods can also be used for this purpose, such as the HHT [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] and WBZ [START_REF] Wood | An alpha modification of Newmark's method[END_REF] methods, which are based on Newmark schemes [START_REF] Newmark | A method of computation for structural dynamics[END_REF] and involve weighting the internal forces and the inertial forces, respectively. These two methods have been generalized by Chung and Hulbert [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method[END_REF]. An interesting review on this subject was published by Fung [START_REF] Fung | Numerical dissipation in time-step integration algorithms for structural dynamic analysis[END_REF] in 2003. Fewer studies have dealt, however, with the use of explicit time integration methods, and the main dissipative schemes of this kind available to our knowledge are those developed by Newmark, Chung and Lee, Zhai, Hulbert and Chung, Tchamwa and Wielgosz [START_REF] Newmark | A method of computation for structural dynamics[END_REF][START_REF] Chung | A new family of explicit time integration methods for linear and non-linear structural dynamics[END_REF][START_REF] Zhai | Two simple fast integration methods for large-scale dynamic problems in engineering[END_REF][START_REF] Hulbert | Explicit time integration algorithms for structural dynamics with optimal numerical dissipation[END_REF][START_REF] Tchamwa | Contributions à l'étude des méthodes d'intégration directe explicites en dynamique non linéaire des structures[END_REF].

Here it is proposed to focus on the family of algorithms involving explicit time integration schemes, especially the Tchamwa-Wielgosz (TW) scheme, because this method makes it possible to account for the high-frequency processes occurring during impact problems. Instantaneous loading conditions were adopted in order to elicit the highest possible frequencies [START_REF] Grolleau | Numerical damping of spurious oscillations: A comparison between bulk-viscosity and dissipative explicit algorithms[END_REF] and thus to study the damping efficiency of this scheme under the most demanding frequency conditions.

Previous studies [START_REF] Soive | Apports à la méthode des éléments finis appliqués aux calculs de structures en dynamique rapide et amortissement numérique[END_REF][START_REF] Grolleau | Une comparaison des schémas d'intégration temporelle explicites de Chung-Lee et Tchamwa-Wielgosz[END_REF][START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF][START_REF] Mahéo | Etude des effets dissipatifs de différents schémas d'intégration temporelle en calcul dynamique par éléments finis[END_REF] have dealt with the damping efficiency at various TW damping parameter values and various calculation time step -to critical time step ratios. It has been established in particular that the greater this ratio is, the more efficient the damping efficiency will be. To complete this research, the influence of the meshing is studied here in the case of homogeneous and inhomogeneous 1-D and 3-D meshes. To prevent the size of calculation time step -to critical time step ratio from affecting the calculations and to obtain the same instantaneous loading conditions throughout, the same ratio was used in all the calculations.

Using viscous or plastic dissipative behavioral laws, as done by Nsiampa with the TW scheme [START_REF] Nsiampa | Comparative study of numerical explicit schemes for impact problem[END_REF], naturally results in a "filtering" process, which causes the highest frequencies in particular to disappear. To prevent the occurrence of this "natural filtering" process, which might make it difficult to assess the efficiency of the TW scheme, a linear elastic conservative law was adopted in the present study.

In Section 2, the conservative Central Differences and dissipative TW schemes are defined in the FEM context. In Section 3, the damping efficiency of the TW method is studied using a novel theoretical approach. It is established that the TW scheme is a decentered scheme and its damping is attributed to a nodal displacement perturbation. In Sections 4 and 5, various homogeneous and inhomogeneous 1-D and 3-D meshes are tested to determine the effects of the refinement and the irregularity of the meshing. The Herezh++ code developed in C++ by Rio [START_REF] Rio | Herezh++: FEM software for large transformations in solids, Dépôt APP[END_REF] at the Laboratoire d'Ingénierie des Matériaux de Bretagne [START_REF]Laboratoire d'Ingénierie des MATériaux de Bretagne (Limatb)[END_REF] was used for the FEM calculations throughout this study. In Section 6, the theoretical results and those obtained on the numerical examples studied are summarized.

Explicit time integration schemes

After the finite element decomposition and time discretization, the weak form of elastodynamics equilibrium leads:

M qn+1 + R int (q n , q n+1 , qn , qn+1 ) = R ext (q n , q n+1 , qn , qn+1 ) (1)
where M is the mass matrix, R int (q n , q n+1 , qn , qn+1 ) is the generalized internal force and R ext (q n , q n+1 , qn , qn+1 ) is the loading force. Taking q to denote the degree of freedom for displacement, q, q and q are the nodal displacement, the velocity and the acceleration vectors, respectively. Time is discretized into time steps t , which are assumed here to be equal for the sake of simplicity, and the subscript n is used at time t = n t.

The conservative central finite difference scheme

The Central finite Difference (CD) method can be written in the form of the following equations (Géradin [START_REF] Géradin | Mechanical Vibrations: Theory and Applications to Structural Dynamics[END_REF]), from which the intermediate times n -1 2 and n + 1 2 have been omitted:

qn+1 =q n + t 2 ( qn +q n+1 ) q n+1 = q n + t qn + ( t) 2 2 qn ( 2
)
where the velocity and the acceleration at time step n can be written:

qn = ( qn+ 1 2 +q n-1 2 ) 2 and qn = ( qn+ 1 2 -q n-1 2 ) t (3)

The dissipative Tchamwa-Wielgosz scheme

The explicit dissipative Tchamwa-Wielgosz (TW) scheme was first presented in 1997 as part of B. Tchamwa's doctoral research work [START_REF] Tchamwa | Contributions à l'étude des méthodes d'intégration directe explicites en dynamique non linéaire des structures[END_REF]. This scheme is based on four parameters. The authors of subsequent consistency studies [START_REF] Tchamwa | Contributions à l'étude des méthodes d'intégration directe explicites en dynamique non linéaire des structures[END_REF][START_REF] Soive | Apports à la méthode des éléments finis appliqués aux calculs de structures en dynamique rapide et amortissement numérique[END_REF][START_REF] Grolleau | Une comparaison des schémas d'intégration temporelle explicites de Chung-Lee et Tchamwa-Wielgosz[END_REF][START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF] have established that the algorithm can be controlled with a single parameter ϕ:

qn+1 =q n + t qn q n+1 = q n + t qn + ϕ( t) 2 qn (4)
Stability has been studied by previous authors. For further details, see Rio et al. [START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF]. However, let us recall that the stability limit of this conditionally stable scheme in the undamped case (ϕ 1) is ω t <

2 ϕ-1 2
where ω is the undamped angular frequency. The evolutions of the spectral radius and the errors in the numerical damping ratio and the relative period have been studied by Nsiampa et al. [START_REF] Nsiampa | Comparative study of numerical explicit schemes for impact problem[END_REF].

Comparison between the CD and the TW schemes

Nodal displacement

In order to compare the displacements obtained upon applying the CD and TW schemes at a given time step, the positions q tw n+1 given by Eq. ( 4) and q cd n+1 given by Eq. ( 2) can be expressed in terms of q n , qn and qn , i.e., in terms of the same position, velocity and acceleration at time step n. Taking the initial position q n and the acceleration qn at time step n, the two schemes differ in the velocity (Eqs. ( 3) and ( 4)). To equalize the two schemes (Eqs. ( 2) and ( 4)) at the time step prior to (n + 1), it suffices to make the two velocities qtw n and qtw n+1 in the TW scheme equal to the velocities qn-1 2 and qn+ 1 2 in the CD scheme, which, after substituting from Eq. ( 3), gives:

qtw n =q n-1 2 =q n - t 2 qn ( 5
)
qtw n+1 =q n+ 1 2 =q n + t 2 qn q tw n+1 = q n + t qtw n + ϕ( t) 2 qn = q n + t qn + (2ϕ -1) ( t) 2 2 qn = q cd n+1 + (ϕ -1)( t) 2 qn (6)
When the parameter controlling the TW scheme is equal to one, ϕ = 1, the displacements calculated at a given time step with the TW and CD schemes will be identical, and if the equation of equilibrium (Eq. ( 1)) does not depend on q and there is no viscosity, the two schemes will be equivalent.

When the parameter controlling the TW scheme is greater than one, ϕ>1, the accuracy of the TW scheme will no longer be of order 2 in terms of a Taylor time development because the expression for the position q tw n+1 will be altered by the extent of (ϕ -1)( t) 2 qn . This expression shows that the damping in the TW scheme depends on the size of the time step used to calculate t [START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF].

The displacement (ϕ -1)( t) 2 q introduces an internal force K(ϕ -1)( t) 2 qn into the equation of equilibrium, Eq. ( 1). In the case of a periodically oscillating displacement q, the displacement (ϕ -1)( t) 2 q will have the opposite sign to that of the displacement q n . The TW scheme therefore yields a generalized expression for the internal force R int , which is K(ϕ -1)( t) 2 qn lower than the internal force obtained in the CD scheme, and this will lead in this case to the damping of the scheme. Note that the higher the frequency of the oscillating displacement, the greater the damping becomes.

Let us now consider an infinite 1-D mesh composed of identical linear bar elements. The natural angular frequency of the nodal displacement is given by ω

0 = 2 K M = 2 c L √ 2 = 2 1 t c
√ 2 where c is the sound speed, K, M, L are the rigidity, the mass and the length of the element respectively and t c = L c is the critical time step of the element. At a given time step value t , let us introduce the calculation time step -to critical time step ratio r t = t t c .T h e expression for the damping ( t) 2 qn of the natural frequency leads to:

( t) 2 qn =-( t) 2 ω 2 0 q n =-( t) 2 ( 2 √ 2 t c ) 2
q n = -8(r t ) 2 q n . Note that at a given r t , this quantity is constant whatever the length of element adopted.

Decentered scheme

A similar expression q dv n+1 for the position at time step n + 1 can also be obtained by performing a second-order Taylor development on the position, taking a decentered expression for the velocity qn+ζ between time steps (n - 1 2 ) and (n + 1 2 ) in terms of the parameter ζ . Upon expressing qn-1 2 and qn+ 1 2 in terms of qn and qn based on Eq. ( 3), we obtain:

qn+ζ = (1 -ζ)q n-1 2 + ζ qn+ 1 2 (7)
Which gives the following equivalent relations: The TW scheme can therefore be said to be a time integration scheme based on a decentered expression for the velocity at time step n with respect to the velocities expressed at half time steps:

q dv n+1 = q n + t qn+ζ + ( t) 2 2 qn = q n + t (1 -ζ)q n-1 2 + ζ qn+ 1 2 + ( t) 2 2 qn = q n + t qn-1 2 + ζ t[q n+ 1 2 -q n-1 2 ]+ ( t) 2 2 qn = q n + t qn + ζ ( t) 2 qn = q cd n+1 + ζ - 1 2 ( t) 2 qn (8)
qn+ζ = (1 -ζ)q n-1 2 + ζ qn+ 1 2 .This gives ϕ = ζ + 1
2 , and since ϕ 1, ζ 1 2 and the scheme is decentered to the right. To obtain the corresponding damping as in previous studies [START_REF] Grolleau | Numerical damping of spurious oscillations: A comparison between bulk-viscosity and dissipative explicit algorithms[END_REF][START_REF] Grolleau | Une comparaison des schémas d'intégration temporelle explicites de Chung-Lee et Tchamwa-Wielgosz[END_REF][START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF], the value of the damping parameter ϕ is taken to be equal to 1.0332 throughout the following numerical study. The scheme is therefore slightly decentered to the right and Eq. ( 7) gives qn+ζ = 0.4668 qn-1 2 + 0.5332 qn+ 1 2 .

One-dimensional damping efficiency under instantaneous loads

Description of the 1-D problem

We take the case of a cylindrical rod with a length of L rod = 200 mm and a cross-sectional area of S = 10 mm 2 , which is discretized into n l = 100 linear bar elements L i in length. One end of the rod is fixed and an instantaneous compression load P 0 = 10 N with a square wave shape is applied to the free end:

-the pressure ramps P 0 at the beginning and -P 0 at the end of the square wave are simulated in 1 t , -and the plateau duration is equal to 48 t .

The bar is assumed to be made of steel with an elastic modulus E = 210 GPa and a density ρ = 8000 kg m -3 .T h e load applied here does not give rise to any plasticity. Linear elastic behavior is therefore adopted throughout the study. The time integration procedure implemented in the code was based on the lumped (diagonal) mass matrix. The global critical time step t c = L min c is about 3.9 × 10 -7 s, where L min is the length of the smallest element and c = E ρ is the sound speed in the material.

To prevent the size of time step from affecting the calculations [START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF], the same global ratio r t = t t c = 90% was used in all the calculations. In this study, the local critical time step t i c = L i c was defined with each element i, which gives the local ratio r i t = t t i c .

Damping efficiency

The theoretical nodal velocity of the node located at the middle of the rod (x node = 100 mm) is plotted in Fig. 1 as a function of time. Two square waves are observed for each round trip, the first is moving from the free to the fixed end (forward, FW) while the second is moving from the fixed to the free end (backward, BW). The sign of the velocity changes each time the wave is reflected from the free end.

For simplicity, the nodal velocities along the 200 mm cylindrical rod are plotted in Fig. 2 for the first, third and ninth forward waves. The 3rd and 9th waves are shifted vertically by 25 mm/s and 50 mm/s, respectively; the red line depicts the theoretical solution of the continuum problem.

In the case of the conservative CD scheme (Fig. 2(a)), the wave seems to oscillate around the theoretical solution during the first RT. For the third and the ninth waves, the oscillations become more and more chaotic (so-called spurious oscillations). Moreover, we observe a monotonic decrease of the slope of the wave fronts. With the TW scheme (Fig. 2(b)), the spurious oscillations are quickly damped and the first forward wave is close to the theoretical solution. The sinusoidal black curves show that the square shape of the wave is progressively altered by the damping of the TW scheme.

We define the overshoot as the difference (absolute value) between the maximum velocity recorded within the rod and the theoretical velocity. The overshoot of the velocity is plotted in Fig. 3. In this figure, the overshoot is only plotted over an interval where the wave is free from any superimposition.

For the TW scheme, the damping efficiency can be observed in Fig. 3: the overshoot of the nodal velocities is completely damped after 60 t. After a propagation time of 2000 t , the damping of the nodal velocity becomes too strong and the maximum velocity does not reach the theoretical value. The value of the damping parameter ϕ = 1.033 appears to be too large in this case. However, ϕ = 1.033 is acceptable within the present framework of our study that focuses on the quick damping of spurious oscillations associated with instantaneous loadings. In Fig. 5(a), the corresponding remaining energy is shown as a function of the time step (green line with diamond marks, "100 el. (vs. t )" referenced). The apparent energy decrease is in line with the previous observations. Moreover, the energy is approximately independent of the reflections of the wave at the bar ends. In the following, the damping of the TW scheme is expressed in terms of the remaining energy.

In the case of the conservative CD scheme, the value of the overshoot is quickly increasing during the first 500 t . After this duration, the overshoot is still oscillating, but the mean value is slowly increasing from 8 mm/sto12mm/s. Furthermore, the energy remains constant during the whole calculation time when using the CD scheme.

Influence of the meshing refinement

The damping obtained with each method was studied on two meshes composed of 100 and 1000 elements with a uniform length. In this case, the local critical time step t i c is equal to the global critical time step t c and r i t = r t . The remaining energy is presented in Fig. 5(a) with respect to time and the number of time steps. When plotted as a function of time, the damping achieved with the TW scheme can be clearly seen here to have depended on the number of elements used in the discretization procedure.

80% of the energy introduced into the system still remained after a propagation time of 7.56 × 10 -5 s and 8.56 × 10 -4 s with meshes consisting of 1000 and 100 elements, respectively. However, when these results were expressed in terms of the number of time steps, the damping turned out to be identical with both meshing configurations. 80% of the energy introduced into the system still remained after a propagation time of about 2120 t with both meshes.

In the case of refined 1-D homogeneous meshes under instantaneous loading conditions, the damping efficiency depends on the number of times the TW dissipative scheme is applied.

Influence of the meshing irregularity

The damping efficiency of the TW scheme was assessed with an inhomogeneous 100-elements mesh. A random shift was applied to each node except for those at both ends. Nodal position was defined as follows:

x new = x old + α 1 Le 2 (1 -α 2 ) (9)
where x new and x old denote the new and previous nodal positions, respectively. α 1 ∈[0, 1] defines the range of nodal positions. α 2 ∈[0, 2] is a random number defining the position of a node in this range. Part of the two 200 mm-meshes is shown in Fig. 4. With the inhomogeneous mesh, the length of the smallest and the largest element is 0.33 mm and 3.81 mm, respectively. Since, the values of the critical time step of the homogeneous and inhomogeneous meshes are 3.91 × 10 -7 s and 0.65 × 10 -7 s, respectively. The remaining energy is presented in Fig. 5(b) with respect to time and the number of time steps. 65% of the energy introduced into the system still remained after a propagation time of 2.24 × 10 -3 s and 4.41 × 10 -4 s with the homogeneous mesh and the inhomogeneous mesh, respectively. However, when these results were expressed in terms of the number of time steps, the damping turned out to be quasi-identical with both meshes. 65% of the energy introduced into the system still remained after a propagation of 6860 t and 7450 t with the homogeneous mesh and the inhomogeneous mesh, respectively. The theoretical analysis in Section 3.1 shows that the damping can be related to the ratio r t by the expression ( t) 2 qn =-8(r t ) 2 q n . In the case of inhomogeneous meshes, the local ratio r i t r t , which means that the damping efficiency decreased.

Three-dimensional damping efficiency under instantaneous loads

Description of the 3-D problem

Here we take the case of a square rod with a length of L rod = 200 mm and a cross-sectional area of S = 64 mm 2 , which is discretized into n e × n e elements across the section and n l elements lengthwise. One end of the rod is fixed and an axial instantaneous compression load P 0 = 10 N with a square wave shape is applied to the free end:

-the pressure ramps P 0 at the beginning and -P 0 at the end of the square wave are simulated in 1 t , -the plateau duration is equal to 48 t .

The elements used were fully integrated 8-node linear hexaedra with 8 integration points.

To prevent the size of time step from affecting the calculations [START_REF] Rio | Comparative study of numerical explicit time integration algorithms[END_REF], the same global ratio r t = t t c was used in all the calculations. In this section, the propagation involves a 3-D strain field. In order to study the influence of the meshing in the longitudinal and the transversal directions, let us take r L t = t t L c and r T t = t t T c to define the calculation time step -to critical time steps ratios in the longitudinal and the transversal directions, respectively.

Influence of the meshing refinement

The damping efficiency of the TW scheme was studied using variously refined 3-D meshes. These meshes can be refined:

-lengthwise and the mesh is named mesh with flattened elements, -across the section of the bar and the mesh is named mesh with elongated elements, -both lengthwise and across the section of the bar and proportionally, in which case the meshes are named meshes with homogeneous elements 1 and 4.

The meshes used are presented in Fig. 6 and their characteristics are listed in Table 1.

For the sake of simplicity, the remaining energy is presented in Fig. 7(a) in terms of the number of time steps alone. Taking the mesh with homogeneous elements 1 as the reference mesh, the differences in the remaining energies are shown in Fig. 7(b).

In the meshes with homogeneous elements 1 and 4, the damping efficiency was found to be practically the same and the difference in remaining energy was about 1%. In this case, r L t = r T t = r t which shows that all these meshes have practically the same damping efficiency. In the mesh with flattened elements, the damping is not so efficient at the beginning. A difference of 2% can be observed in Fig. 7(b) after a propagation of about 2000 t . In this case, r T t <r t and r L t = r t . The axial instantaneous compression load results in a 3-D strain field due to the Poisson's effects, which leads to a dispersive propagation of the wave along the bar. The transversal component of the displacement is perturbed by the change in the meshing refinement, which decreases the damping efficiency, in line with the theoretical study presented in Section 3.1. However, this decrease is smaller than that occurring with the mesh with elongated elements because most of the energy propagates in the longitudinal direction. Other calculations were carried out using a null Poisson's ratio, i.e. without taking the dispersion into account, and applying a homogeneous velocity to the nodes at the free end. These conditions gave 1-D propagation and the damping efficiencies obtained with the meshes with flattened elements and homogeneous elements were identical. It is worth noting that the discretization procedure used in the propagation direction had considerable effects on the damping efficiency.

In the mesh with elongated elements, the damping efficiency decreased sharply at the beginning. A difference of 6.2% can be observed in Fig. 7(b) after a propagation time of about 3000 t . In this case, r L t <r t and r T t = r t . The longitudinal component of the displacement is perturbed by the change in the meshing, which resulted in the largest decrease observed in the damping among all the meshes studied, in line with the theoretical study presented in Section 3.1. Similar results were obtained when a null Poisson's ratio was used in the calculations. This confirms the importance of the discretization procedure used in the propagation direction.

After a propagation time of 10 000 t , the difference between the remaining energies decreased, reaching 1.5% and 3% in the meshes with flattened and elongated elements, respectively. The remaining energies and their differences subsequently tended towards zero.

Influence of the meshing irregularity

Lastly, the efficiency of the TW scheme was studied in the case of an inhomogeneous mesh in the lengthwise direction. Each section is discretized into n e × n e = 4 × 4 elements and the length of the bar was discretized into n l = 100 elements. The same random shift as in Section 4.4 was applied lengthwise.

The remaining energy and the difference of remaining energy are presented with respect to the number of time steps in Figs. 8(a) and 8(b), respectively.

After a propagation time of 3000 t ,a8.3% difference in the remaining energy was observed between the mesh with homogeneous elements 2 and that with inhomogeneous elements 2. In this case, the transversal length was 2 mm and the smallest longitudinal length was 0.33 mm. Thus, r T t <r t and the damping efficiency of the transversal component of the displacement was perturbed by the change of the meshing. In addition, r i t L r t with each element i and the damping efficiency of the longitudinal component of the displacement is also perturbed by the change of the meshing.

Conclusion

The explicit dissipative Tchamwa-Wielgosz scheme is known to be an efficient method of damping spurious oscillations. To complete previous research, a new theoretical study was performed and the influence of the meshing refinement and homogeneity on the damping efficiency was investigated under instantaneous loading conditions.

In the present theoretical study, it was established that the TW scheme was decentered to the right and its damping was attributed to a nodal displacement perturbation. The damping efficiency of the scheme was found to increase with the size of the time step -to critical time step ratio r t = t t c . In the numerical study presented above, various homogeneous and inhomogeneous 1-D and 3-D meshes were tested to determine the influence of the meshing configuration on the damping efficiency under instantaneous loading conditions. In order to be able to compare the various numerical examples, it was necessary to adopt the following two numerical conditions:

-To prevent the size of time step -to critical time step ratio r t = t t c from affecting the calculations and to obtain instantaneous loading conditions throughout, the same ratio was used in all the calculations. -To make the results easier to interpret, the remaining energy was displayed versus the number of time steps.

In the case of the 1-D homogeneous meshes, the results showed that the damping efficiency depends on the number of times the dissipative TW scheme is applied. When plotted as a function of the number of calculation time steps, the damping efficiency was found to be independent of the size of the elements. In the case of the inhomogeneous 1-D meshes, with each element i, the local ratio r i t is lower than the global ratio r t , which confirmed that the damping efficiency decreased. This is in line with the theoretical study, where the damping efficiency was found to be related to this ratio.

In the case of the homogeneous or inhomogeneous 3-D meshes, the results showed the importance of the discretization procedure used in the propagation direction. On similar lines, with each element i, the value of the calculation time step -to critical time step ratio in the longitudinal direction r iL t and that of this in the transversal direction r iT t were always lower than or equal to the value of the global ratio r t , reflecting the decrease in the damping efficiency.

Fig. 1 .

 1 Fig. 1. Square wave propagation versus time.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Wave shape along the 200 mm cylindrical rod. (For interpretation of the colours in this figure the reader is referred to the web version of this article.)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Part of the homogeneous and inhomogeneous 1-D meshes.

Fig. 6 .

 6 Fig. 6. Scale drawing of the homogeneous 3-D meshes.

Fig. 7 .

 7 Fig. 7. Comparison between various homogeneous 3-D meshes. (a) Remaining energy with the 3-D meshes with homogeneous, flattened and elongated elements. (b) Difference in remaining energy, taking the mesh with homogeneous elements 1 as the reference mesh.

Fig. 8 .

 8 Fig. 8. Comparison between homogeneous and inhomogeneous 3-D meshes. (a) Remaining energy with the 3-D mesh with inhomogeneous elements. (b) Difference in the remaining energy taking the mesh with homogeneous elements 2 as the reference mesh.

Table 1

 1 Various 3-D meshes.

	Meshes	Number of elements
		n e × n e × n l
	Homogeneous elements 4 Homogeneous elements 1 Elongated elements Flattened elements	2 × 2 × 50 8 × 8 × 200 8 × 8 × 25 1 × 1 × 200
	Homogeneous elements 2 Inhomogeneous elements 2	4 × 4 × 100 4 × 4 × 100