Mariam Al Harrach
email: mariam.al-harrach@inserm.fr

Pascal Benquet

Fabrice Wendling

Long Term Evolution of Fast Ripples during Epileptogenesis

Keywords: High frequency oscillations, fast ripples, epileptogenesis, hippocampus, kainate model, temporal lobe epilepsy, electrode tissue interface, neural network, Accepted

Fast Ripples (FRs) have received considerable attention in the last decade since they might represent an electrophysiological biomarker of the Epileptogenic Zone (EZ). However, the real dynamics underlying the occurrence, amplitude, and time-frequency content of FRs generation during epileptogenesis are still not well understood. This work aims to characterize and explain the evolution of these features. Stereoelectroencephalography (SEEG) signals recorded in a kainate mouse model of Temporal Lobe Epilepsy (TLE) were processed in order to compute specific FR features. Then realistic physiologically based computational modeling was employed to explore the different elements that can explain the mechanisms of epileptogenesis and simulate the recorded FR in the early and late latent period. Results indicated that continuous changes of FR features are mainly portrayed by the epileptic (pathological) tissue size and synaptic properties. Furthermore, the microelectrodes characteristics were found to dramatically affect the observability and spectral/ temporal content of FRs. Consequently, FRs evolution seems to mirror the continuous pathophysiological mechanism changes that occur during epileptogenesis as long as the microelectrode properties are taken into account. Our study suggests that FRs can account for the pathophysiological changes which might explain the EZ generation and evolution and can contribute in the treatment plan of pharmaco-resistant epilepsies.

Introduction

Temporal Lobe Epilepsy (TLE) is one of the most prevailing categories of adult focal epilepsy [START_REF] Salami | Dynamics of interictal spikes and high-frequency oscillations during epileptogenesis in temporal lobe epilepsy[END_REF][START_REF] Téllez-Zenteno | A Review of the Epidemiology of Temporal Lobe Epilepsy[END_REF]. Patients with TLE suffer from recurrent seizures that highly impact their quality of life. These seizures usually start years succeeding an original brain insult such as status epilepticus, brain injury, febrile seizures, or encephalitis [START_REF] Lévesque | High-frequency oscillations and focal seizures in epileptic rodents[END_REF]. Approximately, 30% of TLE patients are resistant to Anti-Epileptic Drugs (AEDs), making a surgical intervention, to remove the epileptic tissue, the therapeutic option of choice [START_REF] Burnos | The morphology of high frequency oscillations (HFO) does not improve delineating the epileptogenic zone[END_REF]. The outcome of the surgery, however, depends on many factors including a reliable identification of Epileptic Zone (EZ) [START_REF] Burnos | The morphology of high frequency oscillations (HFO) does not improve delineating the epileptogenic zone[END_REF][START_REF] Lévesque | High-frequency oscillations and focal seizures in epileptic rodents[END_REF][START_REF] Zijlmans | High-frequency oscillations as a new biomarker in epilepsy[END_REF]. This calls for reliable biomarkers as well as a deep understanding of the mechanisms that underlies epileptogenesis which can help develop improved mechanism-targeted therapeutic interventions.

Fast Ripples (FRs; activity > 250 Hz) which are pathological transient very high-frequency oscillations are considered as reliable biomarkers of the epileptic activity [START_REF] González Otárula | High-Frequency Oscillation Networks and Surgical Outcome in Adult Focal Epilepsy[END_REF][START_REF] Ortiz | Early Appearance and Spread of Fast Ripples in the Hippocampus in a Model of Cortical Traumatic Brain Injury[END_REF]. They were found to reflect the network's underlying pathological modification during epileptogenesis since they can be recorded weeks before the first spontaneous seizure appearance [START_REF] Bragin | High-frequency oscillations after status epilepticus: epileptogenesis and seizure genesis[END_REF][START_REF] Lévesque | High-frequency oscillations and focal seizures in epileptic rodents[END_REF][START_REF] Li | Unit firing and oscillations at seizure onset in epileptic rodents[END_REF].

FRs were also reported to be initially localized and specific to epileptogenic regions [START_REF] Ibarz | Emergent Dynamics of Fast Ripples in the Epileptic Hippocampus[END_REF][START_REF] Wang | Ripple classification helps to localize the seizure-onset zone in neocortical epilepsy[END_REF]. They were discovered in the recordings of epileptic hippocampus and temporal cortex Local Field Potentials (LFPs) in humans as well as rodents [START_REF] Ibarz | Emergent Dynamics of Fast Ripples in the Epileptic Hippocampus[END_REF]. FRs are explained as the result of a cluster of pyramidal neurons that are firing somewhat asynchronously wherein each neuron is firing at a less frequency than the recorded FRs [START_REF] Demont-Guignard | Distinct hyperexcitability mechanisms underlie fast ripples and epileptic spikes[END_REF]. Since their presence was also linked to functional epileptogenic reorganization [START_REF] Zijlmans | High-frequency oscillations as a new biomarker in epilepsy[END_REF], hypexcitability in the initial SOZ might also induce distant pathological plasticity in connected remote networks that could secondarily lead to the generation of FRs [START_REF] Sheybani | Electrophysiological Evidence for the Development of a Self-Sustained Large-Scale Epileptic Network in the Kainate Mouse Model of Temporal Lobe Epilepsy[END_REF].

Two decades after their discovery, the precise pathophysiological mechanisms of FRs, their temporal and spectral features, and their relation to epileptogenecity are still not well understood. Different studies investigated the mechanisms that can account for epileptogenesis.

It was mainly linked to two elements: 1) progressive decreased inhibition and increased excitation of neurons and interneurons that constitutes the pathological epileptic tissue [START_REF] Demont-Guignard | Distinct hyperexcitability mechanisms underlie fast ripples and epileptic spikes[END_REF][START_REF] Dudek | The Time Course and Circuit Mechanisms of Acquired Epileptogenesis[END_REF][START_REF] Li | Unit firing and oscillations at seizure onset in epileptic rodents[END_REF]. 2) the firing delays of these hyperexcitable neurons as a response to their input from other structures [START_REF] Ibarz | Emergent Dynamics of Fast Ripples in the Epileptic Hippocampus[END_REF]. However, no studies investigated the long term evolution of FRs characteristics during epileptogenesis and its relationship to these mechanisms.

Another aspect, which is usually overlooked in studies investigating interictal spikes and HFOs during epileptogenesis, is the characteristics of microelectrode used to record these activities and its impact on the signal of interest [START_REF] Huneau | Shape features of epileptic spikes are a marker of epileptogenesis in mice[END_REF][START_REF] Ibarz | Emergent Dynamics of Fast Ripples in the Epileptic Hippocampus[END_REF][START_REF] Li | Extrahippocampal high-frequency oscillations during epileptogenesis[END_REF][START_REF] Salami | Dynamics of interictal spikes and high-frequency oscillations during epileptogenesis in temporal lobe epilepsy[END_REF]. Even though microelectrodes allow for high spatial selectivity recordings, they suffer from high impedance due to their reduced size which leads to signal distortion [START_REF] Charkhkar | Chronic intracortical neural recordings using microelectrode arrays coated with PEDOT-TFB[END_REF][START_REF] Stacey | Potential for unreliable interpretation of EEG recorded with microelectrodes[END_REF]. Also, they usually have poor electrical coupling with the tissue resulting in scarring and the formation of glial encapsulation tissue around the microelectrode that contributes to even higher impedance in the FR frequency band [START_REF] Johnson | Repeated voltage biasing improves unit recordings by reducing resistive tissue impedances[END_REF][START_REF] Salami | Dynamics of interictal spikes and high-frequency oscillations during epileptogenesis in temporal lobe epilepsy[END_REF]. Accordingly, the microelectrode properties can dramatically affect both the experimental results analysis as well as the simulated ones.

In this work, we proposed to analyze and explain the evolution of FRs in the CA1 subfield of the hippocampus after KA injection. To do so, first, experimental data were analyzed in order to study the evolution of FRs features throughout a month of recordings in the latent period. A realistic computational model was then used to perform a sensitivity analysis of the different parameters that govern the epileptogenesis which are mainly related to the neural network synaptic properties and the external input (afferent volley of action potentials) from the CA3 hippocampus subfield. Finally, using the results of this analysis realistic FRs at the early and late stages of epileptogenesis were simulated to validate our model.

Methods

Experimental Data Acquisition

Three 80 ± 5-day-old C57BL/6J male mice were used for this study. Experiments were implemented based on the Kainate mouse model of epilepsy [START_REF] Suzuki | Morphogenetic effect of kainate on adult hippocampal neurons associated with a prolonged expression of brain-derived neurotrophic factor[END_REF] respecting the European Communities Council Directive of 24 November 1986 (86/609/EEC). This experiment was approved by the ethics committee on animal experimentation of Rennes (agreement N° APAFIS#2327-2015101914507202).

Surgery

The surgery was performed under anesthesia while the mouse was placed in a stereotaxic frame as shown in Figure 1A. Unilateral localized injection of 50 nl of a 20 mM solution of kainic acid (KA; Sigma-Aldrich) in 0.9% NaCl was performed in the mice. The injection was done into the right dorsal hippocampus (Figure 1C). Afterward, all mice were implanted with one bipolar depth electrode in both the LH and the RH hippocampus (CA1 region, Figure 1C). An additional monopolar electrode was implanted in the cerebellum to act as a reference. After implantation, the five electrode wires were welded into a connector to be attached to the acquisition system. This connector is then fixed on the mouse's skull using dental acrylic cement (Figure 1A). Positions of the KA injection (coordinates from bregma: AP=-2.0 mm, ML=-1.5 mm, DV=-2 mm) and electrodes implantation RH, LH (AP=-2.0 mm, ML=±1.5 mm, DV=-1.9 mm) were chosen using mouse brain atlas [START_REF] Paxinos | Paxinos and Franklin's the Mouse Brain in Stereotaxic Coordinates[END_REF]. The bipolar electrodes were created by twisting two stainless steel wires isolated by Parylene C (125 µm diameter) together as presented in Figure 1B. After the surgery, animals were released into their cages to recover for three days in preparation for the recording sessions.

EEG recording

During the recording sessions, each animal was freely moving in a transparent Plexiglas cage, itself placed in a Faraday cage. Animals were first habituated to their environment for one hour and then Stereo-ElectroEncephalography (SEEG) signals were recorded during 4 hours using an EEG monitoring system (Deltamed TM) and were sampled at 2048 Hz. Recordings were performed every three days (3, 6, 9, 12, 15, 18, 21, 24, 27 and, 30 days) between 1and 5 p.m. as indicated in Figure 1D. A video monitoring was also performed during the recordings.

Computational simulations

Neural Network model

In order to investigate the (patho)physiological mechanisms behind fast ripple generation during epileptogenesis, a computational model of the CA1 subfield of the hippocampus was used. This physiologically relevant microscale model was adapted from the work of [START_REF] Demont-Guignard | Analysis of intracerebral EEG recordings of epileptic spikes: insights from a neural network model[END_REF]. It simulates a 3D volume of the CA1 (Figure 2) that The days in bold represent the signals that were chosen for processing and analysis in this study.

includes principal cells (Pyramidal cells, PYR) along with three types of interneurons (O-LM: Oriens Lacunosum-Moleculare, BIS: Bistratified, and BAS: Basket cells) that have been found to play a role in TLE [START_REF] Demont-Guignard | Analysis of intracerebral EEG recordings of epileptic spikes: insights from a neural network model[END_REF][START_REF] Demont-Guignard | Distinct hyperexcitability mechanisms underlie fast ripples and epileptic spikes[END_REF].

Connections between the different cells are based on Gamma Amino-Butyric Acid (GABA)ergic and glutamatergic synapses defined by realistic connectivity patterns of the CA1 wherein the maximum distance between two connected cells were fixed according to the literature [START_REF] Andersen | The Hippocampus Book, 1st ed[END_REF]. For the PYR cells a reduced two compartments model (dendrite + soma) was used wherein various currents that are of interest in epilepsy are taken into account (Figure 2D). One compartment type model was employed for the modeling of interneurons. The Hudgkin-Huxley formalism was used to compute the compartment's membrane potentials. For more in-depth description of the model please refer to [START_REF] Demont-Guignard | Distinct hyperexcitability mechanisms underlie fast ripples and epileptic spikes[END_REF].

The input of the neural network was defined through a group of external PYR cells which mimics the input that CA1 receives from CA3. Accordingly, the neural network was stimulated via these neurons using sets of quasi-synchronous volleys of action potentials as input. Two types of inputs were used in our simulations: (1) low impact input that leads to background signal generation (theta rhythm) via continuous low amplitude action potential volleys that are uniformly distributed throughout the simulation duration. (2) High impact input which describes the stimulation that, combined with network-specific configuration, leads to HFO generation. The occurrence of the latter is defined by a beta law respecting the usual HFO occurrence frequency depicted in the literature. In both input cases, the action potentials of the volleys are gaussianly distributed.

To generate FRs, a cluster of PYR cells was rendered hyperexcitable (𝑃𝑌𝑅 𝐻 ) as depicted in Figure 2A, 2B in pink. This was done by altering the synaptic characteristics of these cells in The Extracellular Local Field Potential 𝑣 𝑝 at a point 𝑝 of the electrode surface is obtained by solving the forward problem between the neural network and point 𝑝. This was done using the dipole theory as follow:

𝑣 𝑝 = → 𝑑 𝑖 . → 𝑢 𝑖 (𝑘) 4𝜋𝜎𝑟 𝑖 (𝑘) 2 𝑁 𝑖=1
Where 𝜎 is the conductivity, 𝑟 𝑖 (𝑘) is the distance between cell 𝑖 and a measurement point 𝑘 The total potential received by the active surface area of the electrode would be obtained by integrating the dipole equation defined above over the total disc area. This is obtained through the following approximation:

𝑉 = ∯ 𝑣 𝑝 ≈ 1 𝑇 𝑣 𝑝 𝑇 𝑝=1
With T the total number of discretization points of the electrode surface area. The discretization step used is equal to 10𝜇𝑚.

Electrode Tissue Interface (ETI)

Microelectrodes present the advantage of being selective enough to record special oscillation such as HFOs due to their smaller size compared to macroelectrodes. However they suffer from high impedance that can distort the signal of interest depending on its respective characteristics [START_REF] Stacey | Potential for unreliable interpretation of EEG recorded with microelectrodes[END_REF][START_REF] Stacey | Signal distortion from microelectrodes in clinical EEG acquisition systems[END_REF]. In order to reproduce the activity recorded by the stainless steel microelectrodes in vivo, an electrode tissue interface model was defined. This model was described by a modified Randle's model as depicted in Figure 3A. It consists of a spreading resistance (𝑅 𝑠 ) in series with the double layer equivalent elements. The spreading resistance, also called the electrolyte resistance, describes the impedance between the working and the counter electrode [START_REF] Franks | Impedance characterization and modeling of electrodes for biomedical applications[END_REF]. It depends only on the geometrical surface of the electrode as well as the tissue resistivity. The double layer is portrayed by a charge transfer resistance (𝑅 𝐶𝑇 ), which accounts for the charge leakage across the double layer, in parallel with a constant phase angle impedance (𝑍 𝐶𝑃𝐴 ) that portrays its capacitive charging [START_REF] Franks | Impedance characterization and modeling of electrodes for biomedical applications[END_REF]. In addition to the electrode compartment of the ETI we added a tissue compartment that accounts for the encapsulation tissue (glial scarring) around the electrode which occurs as a response to the implanted electrodes and causes a furthermore increase in the ETI impedance [START_REF] Charkhkar | Chronic intracortical neural recordings using microelectrode arrays coated with PEDOT-TFB[END_REF][START_REF] Johnson | Repeated voltage biasing improves unit recordings by reducing resistive tissue impedances[END_REF]. The tissue impedance was found to increase during the first month after the metal electrode implantation then stabilizes in value [START_REF] Charkhkar | Chronic intracortical neural recordings using microelectrode arrays coated with PEDOT-TFB[END_REF]. This tissue response is explained by a mechanical mismatch between the electrode and the surrounding tissue [START_REF] Rivnay | Next-generation probes, particles, and proteins for neural interfacing[END_REF]. It can be described by a membrane capacitance 𝐶 𝑡 of the surrounding cells in parallel with tissue resistance 𝑅 𝑡 that describes the extracellular pathways between these cells together in series with an encapsulation resistance for the encapsulation layer𝑅 𝑒𝑛 . The latter element was added to the spreading resistance previously described ( 𝑅 𝑆 = 𝑅 𝑠 + 𝑅 𝑒𝑛 ). Finally, we added the shunt capacitance (𝐶 𝑠 ) to the ground from the active area of the microelectrode to the input amplifier. The impact of the tissue compartment on the impedance variation with frequency for Day 2, 7 and 24 post implantation versus the pre-implanted microelectrode (without encapsulation tissue) is depicted in figure 3B.

The values of the ETI elements were computed in order to mimic the performance of the experimentally used electrodes. These values were obtained by least square minimization fitting of impedance values measured by Electrochemical Impedance Spectroscopy (EIS). The electrode Bode plots are depicted in Figure 3B. All used equations and fitting values for the electrode and tissue compartments are presented in the supplementary document.

Finally, the simulated signal recorded by the microelectrode was computed by calculating its Fast Fourier Transform (FFT), then multiplying it by the ETI transfer function in the frequency domain and finally reconverting the product back to the time domain using the inverse FFT. An example of the simulated signal before and after filtering due to the ETI interface along with a qualitative comparison with a real recorded signal are presented in Figure 3C.

Geometrical features and simulation description

For this study, we simulated a 300 × 300 × 50𝜇𝑚 volume of the CA1 (see Figure 1).

The volume dimensions were chosen based on the detection range of stainless steel wire electrode diameter used in the experimental protocol. This volume consisted of: 1747 PYR cells and 357 interneurons equally divided into the three groups previously mentioned. The number of PYR cells is computed following realistic density values obtained from the literature [START_REF] Andersen | The Hippocampus Book, 1st ed[END_REF][START_REF] Rapp | Preserved neuron number in the hippocampus of aged rats with spatial learning deficits[END_REF]. For the interneurons number it was computed in order to be equal to 20% of the total cells in the CA1 volume [START_REF] Demont-Guignard | Analysis of intracerebral EEG recordings of epileptic spikes: insights from a neural network model[END_REF]. In addition to these cells, 384 PYR cells (22% of total CA1 PYR cell number)

were also simulated to act as the CA3 input for this network. Each simulation was about 0. This model was used to (1) perform a sensitivity analysis of the different parameters that plays a role in the generation of FRs (2) reproduce the fundamental conditions that leads to recreating FR at the early and late stages of epileptogenesis. To these aims, extensive simulations were performed. For the sensitivity analysis of the tissue characteristics and jitter value, we simulated a total of 20 𝑎𝑛𝑎𝑡𝑜𝑚𝑖𝑒𝑠 × 5 𝑝𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟𝑠 × 42 𝑣𝑎𝑙𝑢𝑒𝑠 = 4200 𝑠𝑖𝑚𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠. In order to recreate the early and late LFP with FRs that were observed in the experimental recordings a total of 3 𝑎𝑛𝑎𝑡𝑜𝑚𝑖𝑒𝑠 × 200𝑠 𝑡𝑖𝑚𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑝𝑟𝑜𝑓𝑖𝑙𝑒𝑠 × 2 𝑠𝑡𝑎𝑡𝑒𝑠 = 1200 𝑠𝑖𝑚𝑢𝑙𝑎𝑡𝑖𝑜𝑛𝑠.

Data analysis and processing

Fast Ripples (FRs) detection and segmentation

Feature extraction and statistical analysis

Four features were extracted from segmented signals in order to describe the FRs evolution during epileptogenesis: the Root Mean Square (RMS) value, FR index that describes the FR to total energy ratio and was calculated as

𝐸 𝐹𝑅 𝐸 𝑡𝑜𝑡𝑎𝑙 = 𝐹𝑅(𝑥) 2 𝑡 𝑒 𝑥=𝑡 𝑖 𝐿𝐹𝑃(𝑥) 2 𝑡 𝑒 𝑥=𝑡 𝑖
, the FR Normalized Spectral Entropy (NSE) which was obtained by-

𝑃(𝑓)𝑙𝑜𝑔2[𝑃(𝑓)] 𝑓𝑠/2 0 𝑙𝑜𝑔2(𝑓𝑠/2)
, where 𝑃(𝑓) is the power spectral density of the FR segment (filtered LFP between 250 and 600 Hz using a Butterworth bandpass filter) at frequency 𝑓 , and the median frequency ( 𝑓 𝑚𝑒𝑑 ). We also computed the duration of the FRs as well as their occurrence. For the occurrence we used an automatic detection algorithm to count the number of FRs per recording session. This algorithm, previously developed in our team, is based on the Gabor transform and is used for the detection and classification of HFOs [START_REF] Jrad | Automatic Detection and Classification of High-Frequency Oscillations in Depth-EEG Signals[END_REF].The RMS variation with respect to the four time points studied was normalized for each mice for the statistical analysis. Since values were not normally distributed, all statistical comparison were done using the nonparametric Kruskal-Wallis and Mann-Whitney-Wilcoxon tests. Results with p-value < 0.05 were considered significant.

Results

Evolution of FR during Epileptogenesis

A total of ~ 2800 FRs segments were extracted from the recorded LFPs for the characterization of fast ripple features evolution during epileptogenesis. The changes in the four features (RMS, FR Index, NSE, and 𝑓 𝑚𝑒𝑑 ) during epileptogenesis are presented in the boxplots of Figure 5. The FR amplitude significantly increased between the early and late stages of epileptogenesis as depicted in Figure 5A. However, the FR Index which represents the ratio of the FRs energy to the total energy decreased over time (Figure 5B) which indicates that even though the FR content of the signal increases the slow wave increases more rapidly resulting in this decreasing profile. The NSE revealed an increase in the regularity and a decrease in the randomness and complexity of the FR power spectrum with time. Finally, the median frequency revealed a significant increase only between the very early stage (third day post KA injection)

and very late stage (24 th day post KA injection) of epileptogenesis and no significant variation afterwards. In addition to these four features, we also found a significant increase in the FR duration and occurrences with time (supplementary Figures S2 andS3).

Mechanistic insights from detailed computational model

In order to explain the experimental results, a realistic computational model was used to perform a sensitivity analysis of the FR features with different physiological and biophysical (electrode characteristics) parameters. FR simulation consisted of two main conditions: 1) creating a hyperexcitable (pathological) cluster inside the neural network volume model and 2) stimulating the network with afferent hyper-synchronized action potentials mimicking the input from CA3. In order to understand the key pathophysiological characteristics that directly impact fast ripple evolution during epileptogenesis extensive simulations were performed (4200 simulations).

To attain the first condition, the GABAergic synapses' reversal potentials are increased from -75mV to at least -50mV creating a pathologically depolarizing GABA cluster (hyperexcitable cluster). However, this parameter alone did not explain the FR shape variation with epileptogenesis. Accordingly, other parameters that were found to directly impact the FR content in the recorded signal were also assessed in this analysis. These parameters were the GABAergic conductance as well as the two glutamatergic synapse conductances (AMPA: αmino-3-hydroxy-5-methylisox-azole-4-propionic acid and NMDA: N-methyl-D-aspartate).

Along with these four physiological parameters, we also studied the variation of FRs with the hyperexcitable cluster size (𝑃𝑌𝑅 ℎ 𝑛𝑢𝑚𝑏𝑒𝑟). The variation of the FR features with respect to these five parameter values are depicted in Figure 6. The three parameters that have a significant impact on the FR amplitude are the GABA reversal potential (EGABA), the AMPA receptors conductance (gAMPA) and the cluster radius. Both RMS and FR Index increased with EGABA and gAMPA. However, with cluster size increase the RMS increased while the FR index decreased. For the other two FR descriptors (NSE and 𝑓 𝑚𝑒𝑑 ) apart from the EGABA they did not display a significant sensitivity to the cluster pathophysiological parameters.

The second condition to simulate FRs as previously mentioned was the synchronization level of afferent volley of action potentials used as input to stimulate the modeled neural network.

This volley of action potentials was described by an intra-burst time jitter which is determined by the Gaussian law standard deviation value (Figure 7A). This jitter reflected the synchronization level of the PYR neurons firing of the CA3 subfield which in turn stimulates the CA1 network. The lower the AP's jitter the higher was the synchronization of the CA1 principal cells responses to the CA3 excitatory drive.

The variation of the FR dynamics with the excitatory input jitter is presented in Figure 7B.

Both the amplitude and energy ratio (FR index) of fast ripples were strongly impacted by the jitter level between 8 and 15 ms. Both features decreased with decreased synchronization of the excitatory APs (Figure 7B). The minimum jitter needed to obtain FRs value was found to be around 8 ms, below this value the higher frequencies start to disappear and were replaced by Interictal Epileptic Spikes (IEDs). No impact was depicted on the NSE and median frequency (Figures 7B andC). Note that the same parameters resulted in a high median frequency standard deviation due to the random aspect of the neural network input. Also, we found that the hyperexcitable cluster radius should be at least 80𝜇𝑚 in order to generate FRs which was equivalent to ~300 hyperexcitable pyramidal cells. All parameter variation boundaries used for this analysis were based on realistic values obtained from the literature [START_REF] Demont-Guignard | Analysis of intracerebral EEG recordings of epileptic spikes: insights from a neural network model[END_REF].

Simulated versus Real FRs

In order to simulate realistic FRs at the early and late stages of epileptogenesis, we need to consider both (patho)physiological mechanisms and biophysical aspects related to the recording microelectrode impedance. Based on the sensitivity analysis performed in the previous section, to reproduce early stage FRs we created a hyperexcitable cluster of ~ 300 𝑃𝑌𝑅 ℎ cells (80 𝜇𝑚 radius) wherein we only changed the EGABA to equal -50mV. In order to switch to late stages the epileptic cluster size was increased to ~800 𝑃𝑌𝑅 ℎ cells (120 𝜇𝑚 radius) at the same time gAMPA conductance was also increased from 8 to 9 mS/𝑐𝑚 2 . In addition to these synaptic alterations, the microelectrode impedance increase due to the encapsulation tissue formation was incorporated in the LFP simulations. The impedance increase between the early and late stages (Figure 2B) is mainly portrayed by an increase in tissue resistance 𝑅 𝑡 and a decrease in its capacitance 𝐶 𝑡 (Figure 8A). For the ETI elements values please refer to the supplementary document.

In order to perform a quantitative comparison between the experimental early (Day 3)/ late (Day 18-24) epileptogenesis stage FRs and the in silico ones, we simulated 3 different anatomies × 200 different stimulations of 0.4 ms where we switch from early to late configuration resulting in a total of 1200 events compared to 1100 events in vivo.

The quantitative comparison, between simulated and recorded FRs at the early and late stages of epileptogenesis is displayed in Figure 8B. All statistical comparisons revealed no significant differences between real and simulated signals in terms of FR amplitude increase between early and late stages (normalized RMS), the ratio of FR to total signal energy, median frequency, and spectral content and complexity. Furthermore, visual inspection of the time-series and spectrograms of early/late simulated signals revealed very similar patterns and characteristics to the recorded ones. An example of early and late simulated compared to real signals is provided in Figure 8C.

Discussion

In this work, we presented an analysis of the long term evolution of FRs during epileptogenesis from the KA injection to the epileptic (chronic) state. This was done as follows:

First, the evolution of FR features in vivo during epileptogenesis using quantitative and statistical analysis was performed. Second, a realistic computational model based on physiology was then used to determine the different parameters that can simulate the underlying mechanisms of epileptogenesis and are henceforth portrayed in the FRs evolution. And last, an encapsulation layer equivalent circuit was added to the ETI model in order to simulate realistic FRs through the early and late stages of the latent period. Accordingly, a combination of pathophysiological (synaptic characteristics) and biophysical (microelectrode impedance) parameter alterations allowed to reproduce the evolution of FRs in silico, demonstrating that FR can mirror the evolution of pathological properties related to epilepsy.

FR morphology and frequency content in the Kainate mouse model

The relationship between FR and seizures is of high importance, especially in a surgical context. Different studies investigated the relationship between HFOs patterns or types and the seizure onset zone or EZ. In some, ripple classification was found to help in SOZ localization and surgical procedures [START_REF] Kerber | Differentiation of specific ripple patterns helps to identify epileptogenic areas for surgical procedures[END_REF][START_REF] Wang | Ripple classification helps to localize the seizure-onset zone in neocortical epilepsy[END_REF]. Others found that HFOs morphology does not result in a better delineation of the EZ [START_REF] Burnos | The morphology of high frequency oscillations (HFO) does not improve delineating the epileptogenic zone[END_REF], and FR were found outside the site of kainite injection even in the cortex [START_REF] Li | Extrahippocampal high-frequency oscillations during epileptogenesis[END_REF][START_REF] Sheybani | Electrophysiological Evidence for the Development of a Self-Sustained Large-Scale Epileptic Network in the Kainate Mouse Model of Temporal Lobe Epilepsy[END_REF]. In this study, we chose not to divide FR based on their pattern or morphology. Instead we studied the global profile variation of amplitude and spectral parameters of all FRs.

Nonetheless, we have to point out the different profiles of FRs that were observed in the experimental recordings. These FR profiles vary between regular/irregular amplitude and frequency content as well as with respect to the background content being quiet or oscillating.

They can be divided mainly into two types of FRs depending on whether they occur on a slow wave or independently from it. Regardless of their type FRs exhibit the same evolution during epileptogenesis as previously demonstrated. Therefore, different profiles of FRs can be explained by physiological as well as biophysical parameters that are directly linked to the amplitude and frequency content of FRs. Firstly, the electrode position with respect to the abnormal tissue or cluster should be considered. Based on the simulated signals the position of the recording electrode including its orientation plays a major role in the characteristics of the recorded FR. Secondly, anatomical variability related to the CA3 input characteristics as well the CA1 pathophysiological parameters variance between subject and between multiple hyperexcitable tissues for the same subject could also impact the arising of different FR patterns. In these same lines, based on our results we did not find a significant increase in the median frequency with time. Still, the median frequency at the 24 th day was higher than the 3 rd day post KA injection. This was observed through an increase in FRs with higher frequencies at the end of the experimental recordings. However, between the 3 rd and 18 th day there was no significant differences which indicates that frequency increase cannot be taken as an indicatory index for epileptogenesis evolution.

Circuit mechanisms of FRs generation during the latent period of epileptogenesis

One of the main findings in our study is that FRs evolution during epileptogenesis depends mostly on the local modification of the cluster of hyperexcitable CA1 neurons characteristics and not the input from CA3. It seems that synchronous action potential discharges of CA3 neurons that induce excitatory postsynaptic potentials on both CA1 neurons and interneurons is not a sufficient condition that explains FR generation. Although combined with altered excitability of a small cluster of CA1 neurons, it leads to FR generation, it does not explain their time evolution between early and late stages of epileptogenesis. Based on our results same excitatory input can generate either early or late FR profiles based on the abnormal excitability/inhibitory properties of the target CA1 network. However, the intra-burst jitter 𝜎 𝑗 of the CA3 afferent volleys needs to be between 8 and 15 ms resulting in a somewhat hypersynchronous firing of the CA1 network (Figure 7). Along these lines and based on our sensitivity analysis performed on the hyperexcitable tissue synaptic properties, we can explain the underlying mechanisms during epileptogenesis evolution on a cellular level. The first major physiological parameter is the abnormal GABAergic inhibition that we portrayed by increasing the potential of pyramidal neurons synaptic GABA (-75mV  -50mV). This abnormal depolarization of GABAergic neurons has been referred to in several papers as the origin of pathological HFOs recorded in humans and rodents [START_REF] Demont-Guignard | Distinct hyperexcitability mechanisms underlie fast ripples and epileptic spikes[END_REF][START_REF] Gulyás | Generation of physiological and pathological high frequency oscillations: the role of perisomatic inhibition in sharp-wave ripple and interictal spike generation[END_REF][START_REF] Lévesque | High-frequency oscillations and focal seizures in epileptic rodents[END_REF][START_REF] Shamas | On the origin of epileptic High Frequency Oscillations observed on clinical electrodes[END_REF]. However, none of these studies explain the evolution of this pathological condition during time. Our simulation study as well as the experimental evolution of FRs seem to support the progressive epilepsy hypothesis mentioned in previous works [START_REF] Dudek | The Time Course and Circuit Mechanisms of Acquired Epileptogenesis[END_REF].

The progressive increase of EGABA over time during epileptogenesis could be explained by an enhanced down-regulation of chloride transporter, such as KCC2, induced by the epileptic activity itself [START_REF] Wang | Double-edged GABAergic synaptic transmission in seizures: The importance of chloride plasticity[END_REF]. Indeed, the large amount of glutamate released by synchronous hyperexitable neurons can potentially activate extrasynaptic NMDAR, and induce calcium influx, activation of mCalpain that trigger KCC2 [START_REF] Wan | M-Calpain Activation Facilitates Seizure Induced KCC2 Down Regulation[END_REF]. KCC2 disruption diminish the capacity and dynamic of chloride [START_REF] Doyon | Chloride Regulation: A Dynamic Equilibrium Crucial for Synaptic Inhibition[END_REF], progressively shift EGABA and favor the emergence or aggravate epileptic [START_REF] Burman | Excitatory GABAergic signalling is associated with benzodiazepine resistance in status epilepticus[END_REF][START_REF] Kurbatova | Dynamic changes of depolarizing GABA in a computational model of epileptogenic brain: Insight for Dravet syndrome[END_REF]. These mechanisms, induced after status epilepticus might be aggravated by a progressive loss of inhibitory interneurons [START_REF] Cossart | Dendritic but not somatic GABAergic inhibition is decreased in experimental epilepsy[END_REF] presented in Figure 9A, as well as increased excitability of pyramidal cells caused by axon sprouting [START_REF] Dudek | The Time Course and Circuit Mechanisms of Acquired Epileptogenesis[END_REF] and abnormal glutamatergic LTP [START_REF] Abegg | Epileptiform activity in rat hippocampus strengthens excitatory synapses[END_REF][START_REF] Franchini | Synaptic GluN2A-Containing NMDA Receptors: From Physiology to Pathological Synaptic Plasticity[END_REF][START_REF] Scimemi | Epileptogenesis is associated with enhanced glutamatergic transmission in the perforant path[END_REF]. This hypothesis is in line with our results wherein a combination of decreased inhibition (EGABA) and increased excitatory conductance (gAMPA) along with an increasing pathological tissue size can reproduce very realistic signals that mirror the recorded ones (Figure 8). Furthermore, this hypothesis can potentially explain the increase in the FR occurrence during the latent period. Considering, the same input for CA3, we found that the number of FRs in a time segment increases when we switched from early stage to late stage configuration as shown in Figure 9.B. Another explanation would be the increase in the excitatory input frequency.

Tissue reaction to implanted microelectrode and FR observability

An important aspect of this study is the biophysical modeling of the recording microelectrode and its direct impact on the FR observability and evolution analysis. This was visible in the case of the RMS and FR index evolution with time where the increase in the impedance between the early and late stages hides the real increase in the FR amplitude and the decrease in its energy ratio in the latent period. Depending on the type and characteristics of the microelectrode this can have critical consequences on the post processing steps as well as interpretation of physiological mechanisms [START_REF] Johnson | Repeated voltage biasing improves unit recordings by reducing resistive tissue impedances[END_REF][START_REF] Mcconnell | Bioimpedance modeling to monitor astrocytic response to chronically implanted electrodes[END_REF]. In our case, this can be observed in the case of the RMS evolution for S2 (Figure 3A) where we can see a decrease in the RMS amplitude between the 3 rd and the 12 th day which is inconsistent with the global pattern. This can be due to an important increase in the glial response between these dates resulting in a decreased observed FR energy which does not portray the real physiological changes. The ETI equivalent tissue parameter values for the encapsulation tissue were adapted from [START_REF] Charkhkar | Chronic intracortical neural recordings using microelectrode arrays coated with PEDOT-TFB[END_REF]. However, the glial response can be different for different microelectrode implantations which was not incorporated in this work and could explain this further. The ETI transfer function magnitude variation for the different periods (pre-implantation and post-implantation days 2, 7, and 28), depicted in supplementary Figure S4, and clearly portrays the impact of the encapsulation tissue on the FR band attenuation. This begs the need for improved types of microelectrodes that have smaller impedance and are more tissue compatible to reduce the mechanical mismatch and scarring that additionally alters the recorded signal.

Conclusions

For the first time, in this work, physiological and biophysical modeling of neural network combined with electrode tissue interface was performed in order to track FR evolution during epileptogenesis. Our findings suggest that FRs evolution can be an indicator of the epileptogenesis stage. The computational modeling indicates that changes in both physiological (neural network) and biophysical characteristics (electrode tissue interface model) can allow for a realistic simulation of the FR at the early and late stages of epileptogenesis. The GABAergic, as well as glutamatergic AMPA synaptic properties of pyramidal cells, govern the underlying circuit mechanism changes in the latent period and not the afferent input from CA3.

Nonetheless, a combination of the two conditions (hyperexcitable cluster and hypersynchronized firing from CA3) are needed to generate FRs. Finally, the encapsulation tissue impact of the microelectrode impedance increase is critical in the analysis of the amplitude and energy of FRs during epileptogenesis.

Figure 1 :

 1 Figure 1: Experimental protocol of intracortical EEG recordings of awake epileptic mices. (A) A mouse during and after electrode implantation. (B) The stainless steel wire microelectrodes used in this protocol. (C) Scheme of the multisite intracortical montage. (D) Timeline of experimental design indicating recording days.

Figure 2 :

 2 Figure 2: Neural network model structure. (A) The simulated 3D volume is designed with the essential elements needed to reproduce the behavior of CA1 subfield (hippocampus). It includes in addition to and pyramidal cells (PYR), oriens-lacunosum moleculare (OLM), Basket (BAS), and Bistratified (BIS) interneurons. The PYR cells in pink are the hyperexcitable PYR (𝑃𝑌𝑅 𝐻 ) cluster. (B) The extracellular field potential (𝑉) observed by the electrode is obtained by solving the forward problem using the dipole theory and then integrating the result over the active electrode surface. (C) The synaptic interactions between the different cells in the network. (D) The reduced CA1 pyramidal neuron two compartments model.

  on the electrode surface, → 𝑢 𝑖 (𝑘) is the corresponding unit vector and → 𝑑 𝑖 is the dipole moment of cell 𝑖 (Please refer to Fig.1 (B)).

Figure 3 :

 3 Figure 3: Microelectrode tissue interface Model. (A) The equivalent circuit comprising two compartments: 1) the electrode compartment that includes the charge transfer resistance 𝑅 𝐶𝑇 in parallel with the constant phase element impedance both in series with the spreading resistance 𝑅 𝑆 ; 2) The glial tissue compartment that includes a tissue capacitance 𝐶 𝑡 in parallel with a tissue resistance 𝑅 𝑡 . (B) Bode plot of the stainless steel microelectrode used in experimental recordings and simulated by the ETI model. The purple band highlights the FR band. (C) An example of the LFP with FRs simulated by the neural network model (left) before (orange) and after (blue) adding the ETI compared with an experimental signal sample (right).

  4 s long with sampling frequency of 25000 Hz. The final signals after ETI filtering were later down-sampled to match the experimental signals' sampling frequency of 2048 Hz.

For

  the experimental signals, LFPs segments with FRs were manually selected by inspecting the bandpass filtered signals (Butterworth; 250-600 Hz). Artifact-free isolated interictal FRs of about 0.3~0.5 ms were detected and exported for segmentation. The same FRs automatic segmentation pipeline was applied on both experimental and simulated signals. It consists of first computing the spectrograms of the preselected (experimental signals) or simulated FRs segments in the FR band (250-600 Hz) as depicted in Figure 4B. This spectrogram is then considered as an image and segmented using the watershed algorithm. This allowed us to automatically segment the FR part of the signals for feature extraction. In brief, this algorithm consists of multiple steps that processes the grey level image obtained from the spectrogram. It includes binarization, noise reduction, sure background and foreground as well as contour area detection, labelling and finally watershed segmentation. The final outcome of this algorithm is presented in Figure 4B. This final boundary marking is used to obtain the start and end of the FR segment (Figure4A). Supplementary FigureS1shows an example of steps of this algorithm applied on the spectrogram image in the FR band. Figure4Cdepicts an example of the time evolution of experimentally recorded FRs for one of the mice after segmentation using the watershed algorithm and realignment.

Figure 4 :

 4 Figure 4: FR segmentation and evolution. A: An example FR segmentation with the raw data segment (top) and the FR band-filtered segment (bottom) of experimental signal before and after FR marking. B: the corresponding spectrogram (left) with the watershed pipeline result (right). C: Evolution of manually detected FRs in one KA treated mouse throughout the recording month.

Figure 5 :

 5 Figure 5: FR evolution during epileptogenesis for the three mice (S1, S2, and S3). Normalized RMS (A), FR index (B), normalized spectral entropy (C), and median frequency (D) variation during the third, 12 th , 18 th , and 24 th day post kainic acid injection. The results of statistical comparisons between the different recording days using the Kruskal-Wallis test are represented as follow: 𝑛𝑠: 𝑛𝑜𝑛𝑠𝑖𝑔𝑛𝑖𝑓𝑖𝑐𝑎𝑛𝑡, * : 1.0 × 10 -2 < 𝑝 ≤ 5.0 × 10 -2 , * * : 1.0 × 10 -3 < 𝑝 ≤ 1.0 × 10 -2 , * * * : 1.0 × 10 -4 < 𝑝 ≤ 1.0 × 10 -3 , * * * * : 𝑝 ≤ 1.0 × 10 -4 .

Figure 6 :

 6 Figure 6: Impact of hyperexcitable cluster tissue parameters: the GABAergic synaptic potential (EGABA) and conductance (gGABA), the glutamatergic conductances gNMDA and gAMPA of the PYR cells in the cluster and the Cluster Radius/ 𝑃𝑌𝑅 ℎ total number on FR features: Normalized RMS, FR index, normal spectral density (NSE) and the median frequency.

Figure 7 :

 7 Figure 7: The impact of the afferent input of the CA3 PYR cells on the fast ripples features. (A) The input of the CA3 PYR cells is modeled as a volley of action potentials (APs) that are distributed respecting a Gaussian law with a specific standard deviation 𝜎 𝑗 (that we call jitter) which describes the level of synchronization of these APs as well as the PYR cell responses. (B) The variation of the FR features (normalized RMS, FR index, NSE and𝑓 𝑚𝑒𝑑 ) with the jitter value between 8 and 15 ms. (C) An example of the FR spectral content and LFP shape with varying Jitter.

Figure 8 :

 8 Figure 8: Simulation of Fast Ripples (FRs) at the early and late stages of epileptogenesis. (A) Physiological and biophysical conditions required to model the transition between early and late stages. (B) Violin plots representing the normalized RMS, FR index, NSE and median frequency at early and late stages for simulated compared to real FRs. All statistical comparisons (Kruskal-Wallis tests) between simulated and experimental signal features were ns (nonsignificant). (C) Time-amplitude and time-frequency qualitative comparison between real and simulated signals at early and late stages of epileptogenesis. 𝑛𝑠: 𝑛𝑜𝑛𝑠𝑖𝑔𝑛𝑖𝑓𝑖𝑐𝑎𝑛𝑡; 𝑝 > 0.05.

Figure 9 :

 9 Figure 9: Pathophysiological changes of CA1 synaptic connections. (A) Hypothetical scenarios of normal versus epileptic synaptic tissue circuits in the early and late stages of the latent period. The early stage is described by a loss of GABAergic interneurons inhibition of principal cells. The Late stage is a combination of more loss of GABAergic interneurons and new excitatory connections between PYR cells as a result to axonal sprouting leading to furthermore increased excitability in a small cluster of the neural network. (B) An example of increasing glutamatergic conductivity and the hyperexcitable cluster size leads to the increase in the simulated FRs occurrence.

Funding Sources

M. Al Harrach was funded by the French research agency (01/10/2019, 30/09/2020, N° ANR-18-CE19-0013-01). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 855109, ERC-SyG2019)

Conflict of Interests

None of the authors have potential conflicts of interest to be disclosed.