Seizure epicenter depth and translaminar field potential synchrony underlie complex variations in tissue oxygenation during ictal initiation

&NA; Whether functional hyperemia during epileptic activity is adequate to meet the heightened metabolic demand of such events is controversial. Whereas some studies have demonstrated hyperoxia during ictal onsets, other work has reported transient hypoxic episodes that are spatially dependent on local surface microvasculature. Crucially, how laminar differences in ictal evolution can affect subsequent cerebrovascular responses has not been thus far investigated, and is likely significant in view of possible laminar‐dependent neurovascular mechanisms and angioarchitecture. We addressed this open question using a novel multi‐modal methodology enabling concurrent measurement of cortical tissue oxygenation, blood flow and hemoglobin concentration, alongside laminar recordings of neural activity, in a urethane anesthetized rat model of recurrent seizures induced by 4‐aminopyridine. We reveal there to be a close relationship between seizure epicenter depth, translaminar local field potential (LFP) synchrony and tissue oxygenation during the early stages of recurrent seizures, whereby deep layer seizures are associated with decreased cross laminar synchrony and prolonged periods of hypoxia, and middle layer seizures are accompanied by increased cross‐laminar synchrony and hyperoxia. Through comparison with functional activation by somatosensory stimulation and graded hypercapnia, we show that these seizure‐related cerebrovascular responses occur in the presence of conserved neural‐hemodynamic and blood flow‐volume coupling. Our data provide new insights into the laminar dependency of seizure‐related neurovascular responses, which may reconcile inconsistent observations of seizure‐related hypoxia in the literature, and highlight a potential layer‐dependent vulnerability that may contribute to the harmful effects of clinical recurrent seizures. The relevance of our findings to perfusion‐related functional neuroimaging techniques in epilepsy are also discussed. HighlightsLaminar LFP synchrony during seizures influences tissue oxygenation responses.Seizures with deep epicenters were associated with reduced cross‐laminar synchrony.Seizures with deep epicenters were associated with tissue hypoxia.Effects observed despite increases in CBF and conserved neural‐hemodynamic coupling.Suggests possible laminar susceptibility to tissue hypoxia during ictal events.


Introduction
Understanding the physiological cascade of blood flow and metabolism during epileptic seizures can provide insights into the pathophysiology and pathological consequences of such events, and enhance the interpretation of diagnostic perfusion-related neuroimaging techniques in epilepsy. There is ongoing debate about whether the supranormal metabolic demands made by epileptic seizures are adequately met by local delivery of oxygenated blood through functional hyperemia. The notion that cerebral perfusion changes might not sustain these intense energy demands and lead to ischemia-hypoxia and neurodegeneration was advocated by a number of early studies (Spielmeyer, 1930;Davis et al., 1944;Cooper et al., 1966;Plum et al., 1968;Dymond and Crandall, 1976). However, this hypothesis fell largely out of favor with subsequent reports of hyperoxia during repetitive seizures with blood flow exceeding metabolic demand (Franck et al., 1985;Ingvar, 1986), and excitotoxic-related neuronal injury in the absence of hypoxia (for review see Meldrum, 2002). Notwithstanding, recent work has begun to come full-circle and provided renewed support for cerebrovascular dysfunction as a parallel phenomenon in epilepsy. Transient decreases in brain tissue oxygenation have recently been observed in the seizure onset zone during onset of ictal-like discharges (Bahar et al., 2006;Zhao et al., 2009), that are intensified in brain regions distal to arteries and arterioles, and can even serve to predict seizure duration (Zhang et al., 2015(Zhang et al., , 2017. Other work has also provided support for ischemia-hypoxia during ictal events by demonstrating that ictal neurodegeneration is spatially coupled to abnormal capillary blood flow (Leal-Campanario et al., 2017), in keeping with computational models suggesting that altered capillary flow patterns can adversely affect oxygen extraction efficacy (Jespersen and Østergaard, 2012). These findings have therefore provided important insights into neurovascular ictal dynamics and the spatial role of microvasculature on these responses in the superficial plane of the cortex, although the question why tissue oxygenation deficits during ictal events have not been routinely observed remains unanswered. One relevant, but neglected, line of inquiry concerns the influence of the differential laminar evolution of ictal activity on accompanying hemodynamic and tissue oxygenation dynamics, given depth-specific variations in angioarchitecture (Masamoto et al., 2004;Blinder et al., 2013), tissue oxygen tension (Feng et al., 1988), and neurovascular coupling (Poplawsky et al., 2015). Laminar differences in seizure initiation and propagation could therefore significantly modulate ensuing cortical tissue oxygenation, thereby providing a mechanistic basis for variability in seizure outcomes and seizure-related oxemic observations in the literature. Unfortunately, previous methodological approaches have been unsuitable and unable to test this prospect directly.
Elucidating the extent to which the relationship between hemodynamic changes and underlying neural activity, known as neurovascular coupling, is conserved during epileptic seizures, compared to normal conditions, is also essential to the correct application of diagnostic perfusion-related neuroimaging techniques, such as blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI). BOLD-fMRI signals are typically interpreted in terms of underlying neural activation within the framework of a general linear model (Logothetis et al., 2001) that is invariant across health and disease, and are being increasingly employed in combination with electroencephalogram (EEG) recordings to localize hemodynamic correlates of electrographic ictal activity (Salek-Haddadi et al., 2002;Tyvaert et al., 2008). Furthermore, ascertaining the stability of the relationship between blood flow and volume is central to refining assumptive estimates of cerebral metabolic rate of oxygen consumption (CMRO 2 ) transients from human 'calibrated' fMRI (Davis et al., 1998;Kida et al., 2007), which can provide more faithful measures of neural activity compared to BOLD signals, and are therefore of potential value to ictal-fMRI studies. Nevertheless, whether normal neurovascular and blood-volume coupling is conserved during epileptic seizures remains to be fully clarified.
Accordingly, we have sought to address these unresolved questions through development of a novel multi-modal methodology that enables concurrent measurement of tissue oxygenation, blood flow and total hemoglobin concentration, alongside laminar neural activity, in the urethane anesthetized rat. We aimed to examine the role of laminardependent initiation and propagation of recurrent acute neocortical seizures, induced by 4-aminpyridine, on ensuing cerebrovascular responses. Further, through comparison to functional activation by whisker stimulation and hypercapnia, we endeavored to describe the response dynamics of all variables along a continuum of brain activation, and interrogate the stability of neural-hemodynamic and blood flow-volume coupling during ictal-like discharges.

Animal preparation and surgery
All procedures were conducted with approval from the UK Home Office under the Animals (Scientific Procedures) Act of 1986. Female hooded Lister rats (total N ¼ 12 weighing 250-350 g) were kept in a 12hr dark/light cycle environment at a temperature of 22 C, with food and water ad libitum. Animals were anesthetized with intraperitoneal injection of urethane (1.25 g/kg) and anesthetic depth was monitored throughout surgery and experiments. Urethane anesthesia preserves excitatory/inhibitory synaptic transmission (Sceniak and MacIver, 2006) and neurovascular coupling (Berwick et al., 2008), as well as providing persistent anesthesia reminiscent of natural sleep (Pagliardini et al., 2013). Core/body temperature of animals was maintained at 37 C through use of a homoeothermic blanket and rectal probe (Harvard Apparatus). Animals were tracheotomized and artificially ventilated with medical air. End-tidal CO 2 was monitored throughout experiments. Blood-gas measurements were acquired during experiments to inform and regulate ventilator settings (e.g. air flow, breathing rate, tidal volume), so as to maintain blood-gas parameters within normal physiological ranges (mean AE SEM: PO 2 : 91.1 AE 4.94 mmHg, PCO 2 : 34.7 AE 3.69 mmHg, SO 2 : 96.8 AE 0.73%, pH: 7.4 AE 0.03). The left femoral artery and vein were cannulated to allow the measurement of arterial blood pressure and phenylephrine infusion (0.13-0.26 mg/h to maintain normotension between~100 and 110 mmHg), respectively . The animal was secured in a stereotaxic frame throughout experimentation and the skull overlying the right parietal cortex thinned to translucency so as to expose the somatosensory cortex. A layer of cyanoacrylate glue was subsequently applied over this region to reduce optical specularities from the brain surface during imaging and preserve skull integrity.

Two-dimensional optical imaging spectroscopy
Two-dimensional optical imaging spectroscopy (2D-OIS) was used to produce 2D images over time of total hemoglobin (Hbt) concentration, using a heterogeneous tissue model as described previously (Berwick et al., , 2008Kennerley et al., 2009). Briefly, the cortex was illuminated at four different wavelengths (495 AE 31 nm, 559 AE 16 nm, 575 AE 14 nm and 587 AE 9 nm FWHM) using a Lambda DG-4 high speed filter changer (Sutter Instrument Company, Novata, CA, USA) and image data recorded at 8 Hz using a Dalsa 1M30P camera (Billerica, MA, USA, each pixel representing~75 μm 2 ) synchronized to the filter switching (2 Hz/wavelength). Data were subsequently subjected to spectral analysis consisting of a path length scaling algorithm (PLSA) employing a modified Beer-Lambert law in conjunction with a path-length correction factor for each wavelength used, based on Monte Carlo simulations of light transport through tissue (Wang et al., 1995;Berwick et al., 2005Berwick et al., , 2008. Baseline concentration of Hbt was estimated at 104 μmol/L (Kennerley et al., 2009). The somatosensory barrel cortex was localized prior to experiments to guide implantation of a multi-channel electrode and multi-sensor probe. Here, the left (contralateral) whisker pad was briefly stimulated electrically using subcutaneous electrodes (30 trials, 2s, 5 Hz, 1.2 mA, 0.3 ms pulse width) and resultant 2D-OIS data averaged and subjected to the above spectral analysis. Spatiotemporal changes in Hbt were analyzed using statistical parametric mapping (SPM) which produced a Z-score activation map that was used to localize the cortical region activated by whisker stimulation (as described previously, Harris et al., 2013).

Multi-sensor probe recordings
A multi-modal sensor comprising of laser Doppler flowmetry (LDF) and luminescence-based oxygen probe (450 μm diameter, Oxford Optronix, UK) was used to obtain concurrent measures of cerebral blood flow (CBF) and tissue oxygenation (tpO 2 ), respectively (Trübel et al., 2006). The laser Doppler probe provided a measure of relative changes in CBF, while the luminescence probe provided an absolute measure of tpO 2 (AE0.1 mmHg, 0.5 mm diameter sphere sampling volume), respectively. A 750 nm low-pass filter was attached to the 2D-OIS camera lens to prevent cross-talk from the LDF probe operating at 830 AE 10 nm, and the transmission curve of the filter accounted for during spectral analysis. The multi-sensor probe was attached to a stereotaxic holder and inserted into the somatosensory 'whisker barrel' cortex under microscope guidance to a depth of 500 μm (layer 4) (Fig. 1A). Signals were amplified and digitized using a 2-channel OxyFlo Pro system (Oxford Optronix, UK) and recorded using a CED Power 1401 and Spike2 software (Cambridge Electronic Design, Cambridge, UK). Sensor data were recorded continuously during each experiment and sampled at 1 Hz.

Laminar electrophysiology
A multi-channel microelectrode (16 channels with 100 μm spacing, site area 177 μm 2 , 1.5-2.7 MΩ impedance; Neuronexus Technologies, Ann Arbor, MI, USA) was placed in a stereotaxic holder and carefully implanted into the right barrel cortex, using a micro-manipulator and microscope, to a depth of 1500 μm (Harris et al., 2013;Boorman et al., 2015), as close as possible to the multi-sensor probe (Fig. 1A), and such that the upper-most electrode was located at the cortical surface. During epilepsy experiments the multi-channel microelectrode was coupled to a fluidic probe (Neuronexus Technologies, Ann Arbor, MI, USA) to enable drug delivery at a cortical depth of 1500 μm. Electrodes, and multi-sensor probes, were allowed to stabilize for at least 1 h prior to recording experiments. The microelectrode was connected to a preamplifier and data acquisition device (Medusa BioAmp/RZ5, TDT, Alachua, FL, USA). Neural data were recorded continuously throughout each experiment and sampled at 6 kHz.

Experimental Paradigm
Multi-modal responses to whisker stimulation (N ¼ 5) were investigated using a 15 trial protocol, each delivering a 16-s train of electrical pulses (5 Hz, 1.2 mA intensity, 0.3 ms pulse width) with an inter-trial 1500 μm depth at 4AP infusion site. C) Percentage CBF changes in all animals (N ¼ 7) during first 2500s of recurrent seizure activity. D) Marked variability in percentage tpO 2 changes in all animals (N ¼ 7) during first 2500s of recurrent seizure activity (color coded as in C). E) Normalized and averaged LFP layer 4 sink responses to 16s whisker stimulation (1.2 mA, 5 Hz, N ¼ 5, top) and accompanying averaged CBF (red) and tpO 2 (blue) percentage changes during stimulation (bottom, N ¼ 5). Inset shows a representative laminar LFP response profile to first pulse in stimulation train (averaged across trials). F) Example end tidal CO 2 trace (top), normalized instantaneous LFP power as a function of depth (middle, averaged across animals, N ¼ 5), and averaged CBF (red) and tpO 2 (blue) percentage changes (bottom, N ¼ 5), during 5% hypercapnia. G) Same as in F but for 10% hypercapnia.
interval of 70 s, to the right whisker pad following a 10-s baseline period (Harris et al., 2013). Hypercapnia experiments (N ¼ 5) were conducted after the stimulation protocol and consisted of seven interleaved trials of 190s duration with normocapnia on odd trials, and hypercapnia on trials 2, 4 and 6. No stimulation was delivered during hypercapnia experiments. Hypercapnia was induced in two sequential experimental runs where the fraction of inspired CO 2 (FiCO 2 ) was increased to either 5% (N ¼ 5) or 10% (N ¼ 5) . Hypercapnia trials for each condition were extracted with a 20s normocapnic baseline and 190s post-hypercapnic period. As hypercapnia has been demonstrated to have a potent anti-convulsant action (Tolner et al., 2011), we conducted acute seizure experiments in a separate group of animals (N ¼ 7). Fluidic microelectrodes were loaded with the potassium channel blocker 4-aminopyridine (4-AP, 15 mM, 1 μl, Sigma, UK) to induce recurrent focal seizure-like discharges in right barrel cortex (Harris et al., 2014a(Harris et al., , 2014b. 4-AP was intra-cortically infused at a depth of 1500 μm over a 5-min period (0.2 μl/min), using a 10 μl Hamilton syringe and syringe pump (World Precision Instruments Inc., FL, USA) (Harris et al., 2014a(Harris et al., , 2014b. Multi-sensor, electrophysiological and 2D-OIS measures were recorded concurrently for 7000s with 4-AP injection initiated following a 280s baseline period.

Data and statistical analysis
All data analysis was conducted using custom written Matlab™ scripts (Mathworks, Natick, MA, USA). Laminar neural recordings were bandpass filtered between 1 and 300 Hz using a 300 th order finite impulse response filter and de-meaned to yield local field potential (LFP) data. Cross-laminar synchronization of neural activity was investigated through use of a mean phase coherence (MPC) algorithm, described in detail previously (Mormann et al., 2000;Schevon et al., 2007). The MPC algorithm was sequentially applied to pairs of LFP signals acquired from two different electrode sites along the multi-channel array (16 channels yielding 120 unique combinations of paired electrode data). Briefly, the instantaneous phase of each LFP signal was first calculated using the Hilbert transform, with instantaneous phase values of both signals projected onto the unit circle in the complex plane, and the mean phase coherence computed according to the expression (Mormann et al., 2000): where ϕ x and ϕ y are the instantaneous phases of two input signals, N is the number of samples of the data set, and t is the time index. Euler's formula was then used to transform the above equation, such that MPC is restricted to the interval [0,1], where MPC ¼ 0 denotes complete asynchrony and MPC ¼ 1 represents strict phase locking between both signals (Mormann et al., 2000). Importantly, MPC is not dependent on the amplitude of the input signals, and is insensitive to constant phase delays between two signals such as that arising due to neural transmission latency between two recording locations (Mormann et al., 2000;Schevon et al., 2007). Synchrony between paired LFP signals was obtained by computing the MPC in 200 ms time-bins, applying a 5-point moving average filter, and downsampling the resultant timeseries to a sampling frequency of 1 Hz (with use of an 8th order lowpass Chebyshev Type I infinite impulse-response filter to minimize aliasing). Each MPC timeseries was normalized to its mean MPC value during the 280s baseline period prior to 4-AP infusion. LFP synchrony between different cortical layers was examined by grouping and averaging normalized MPC data relevant to the source-target cortical laminae of interest, defined according to cortical depth; supragranular (SG < 400 μm), granular (GR, 400-700 μm), infragranular (IG, 800-1100 μm) and deep infragranular (dIG, 1200-1500 μm).
Cerebral metabolic rate of oxygen consumption (CMRO 2 ) was calculated from CBF and tissue oxygen tension (tpO 2 ) measurements according to the following relationship (Gjedde, 2005): where tpO 2 is the tissue oxygen tension and P 50 and h are the halfsaturation tension and Hill coefficient of the oxygen-hemoglobin dissociation curve, 36 mmHg and 2.7, respectively. C a is the arterial oxygen concentration (8 μM/ml). The effective diffusion coefficient of oxygen in brain tissue, L, was calculated for each experiment using trial-averaged baseline tpO 2 values, and baseline values of CBF and CMRO 2 previously reported in rats by our laboratory and others, 47 ml/100 g/min (He et al., 2007) and 219 μM/100 g/min (Zhu et al., 2002), respectively. L was found to be 7.27 AE 0.89 μM/100 g/min/mmHg (N ¼ 22 separate stimulation, hypercapnia and seizure experiments).
Time-series of Hbt responses during the different experimental conditions were extracted by selecting a region of interest (ROI) over acquired Hbt image data. As propagatory patterns during seizures can be highly variable and dynamic (Ma et al., 2012), we elected to manually position a ROI such that it was immediately adjacent to the multi-sensor probe and depth electrode, and captured a well-mixed vascular compartment (i.e. arterial, venous and parenchymal). The mean Hbt time-series of all pixels in the ROI were subsequently averaged for each animal and normalized to the baseline period prior to 4-AP infusion.
Since evoked cortical responses to stimulation exhibit a great deal of stability and reproducibility due to the fine control of stimulation parameters, Hbt image data acquired during each animal's sensory stimulation condition were averaged over the period 0-20s after stimulus onset, and pixels possessing >80% of the maximum change in evoked Hbt automatically selected. The mean Hbt time-series of all pixels within the ROI were then extracted for whisker stimulation and hypercapnia conditions in turn. Hbt data from each animal during stimulation and graded hypercapnia were normalized to the mean baseline value across trials in each respective condition.
Recurrent seizure onset and offset in electrophysiological recordings was estimated by conducting power spectrum analysis on broadband field potential data from the deepest electrode channel (i.e. corresponding to the 4-AP infusion site at 1500 μm, representative example shown in Fig. 1A) using a Gabor transform. This consisted of short-term Fourier transform and 1s Gaussian window function applied sequentially to the entire timeseries (Harris et al., 2014b), and identifying the time-points where total broadband power first exceeded, and latterly fell below, 5 standard deviations above the mean baseline (i.e. 280s period prior to 4-AP infusion). Recurrent seizure duration was 5325.9 AE 825.1s (N ¼ 7) and the seizure initiation period was defined as a 300s time-window following ictal onset. The locus of maximum activation during recurrent seizures (i.e. seizure epicenter) was found by calculating the total instantaneous power (integral of squared magnitude) of broadband laminar field potential data during the entire ictal period, and identifying the electrode channel exhibiting the greatest change in magnitude. Seizure epicenters were located on average at a depth of 1143 AE 113 μm (N ¼ 7), consistent with cortical layer 5b and our previous work (Harris et al., 2014a). Seizure-related LFP response profiles were obtained by identifying depth-negative peak responses during ictal activity whose magnitude was less than 5 standard deviations below the mean baseline activity. LFP profiles were then inverted and normalized to the largest peak response in each animal.
The layer 4 current sink in laminar neural recordings was localized by averaging the whisker evoked response to the first stimulus pulse across trials and identifying the electrode channel associated with the earliest robust negative field potential deflection; in this study the L4 current sink was located 580 AE 111 μm (N ¼ 5) below the cortical surface. Trial-bytrial Layer 4 LFP data were then extracted for each animal and normalized such that the peak (depth-negative) response evoked by the first pulse of the stimulation train was unity across trials (Harris et al., 2013) S.S. Harris et al. NeuroImage 171 (2018) 165-175 (representative animal shown in Fig. 1E). Stimulation-related LFP response profiles were subsequently extracted by identifying the peak response to each electrical pulse in the stimulation train. Electrophysiological changes during hypercapnia were assessed by computing the instantaneous power of laminar LFP data and averaging over 1s time-bins to match the sampling frequency of the multi-sensor probe. Instantaneous power data from each electrode channel was normalized to the mean baseline value during the three trials in each hypercapnia condition.
In order to account for differing data ranges across experimental conditions, trial-by-trial relationships between variables during stimulation, hypercapnia and seizures were analyzed following standardization of the data using z-scores (i.e. with mean 0 and standard deviation 1). Timeseries are presented as a percentage change from baseline, and errors are given as the standard error from the mean (SEM), unless otherwise stated. Linear regression models were fitted using ordinary least squares, and differences between regression slopes examined using analysis of covariance with Tukey-Kramer correction for multiple comparisons. A 1-tailed Wilcoxon signed rank test was used to test significant changes in individual variables, and significant differences between conditions examined using a Kruskal-Wallis test with Tukey-Kramer correction for multiple comparisons. Exploratory correlation analyses involving micro-array data were corrected for using the Benjamini and Yekutieli (2001) procedure for controlling the false discovery rate (FDR) that is suitable for all forms of test-statistic dependency. Multiplicity-adjusted p values < 0.05 were defined as statically significant, with p < .01 deemed to be highly significant.

Variable tpO 2 during recurrent seizures is related to laminar-dependent changes in LFP activation and synchrony
We next investigated candidate processes that might underpin the observed variability in tpO 2 during recurrent seizures. CBF onset dynamics, which typically consisted of robust and rapid increases, was first ruled out as a possible driver of the observed tpO 2 variability, since no significant correlation between both variables was observed during ictal onset (p ¼ .24, N ¼ 7) (example from two different animals with tpO 2 changes of opposite polarity are shown in Fig. 3A-B). Close inspection of laminar LFP waveforms, however, suggested that the vertical propagation of ictal activity differed between animals where tpO 2 increased or decreased, with the latter associated with a comparatively decreased activation in upper cortical layers (example shown in Fig. 3C-D). Mean phase coherence (MPC) analysis of LFP data from paired microelectrode channels was therefore conducted to more reliably assess seizure propagation though variations in laminar synchrony during seizure initiation. This revealed characteristic differences between cases where tpO 2 increased or decreased during seizure onset (example shown in Fig.  3E-F). Notably, the MPC between upper and lower cortical layers exhibited an increase during seizures accompanied by an increase in tpO 2 (red areas in Fig. 3E), whereas seizure onsets associated with a decrease in tpO2 showed the reverse (blue areas in Fig. 3F).
The relationship between these coherence patterns and associated changes in tpO 2 was confirmed using correlation analysis, which revealed MPC between supragranular and underlying laminae, and tpO 2 , to be significantly linearly correlated (SG-GR, Pearson's r ¼ 0.91, p ¼ .02; SG-IG, Pearson's r ¼ 0.91, p ¼ .02; SG-dIG, Pearson's r ¼ 0.93, p ¼ .02; N ¼ 7, FDR corrected for multiple comparisons), at variance with other cortical laminae combinations (Fig. 4A). Cross-laminar changes in MPC during seizure onset, across all cortical depths studied, were also tightly coupled to changes in tpO 2 (Pearson's r ¼ 0.91, p ¼ .004, N ¼ 7, Fig. 4B). We also found that the seizure epicenter depth was inversely monotonically related to mean tpO 2 changes (Spearman's rho ¼ À 0.9, p ¼ .01, N ¼ 7, Fig. 4C left panel), such that seizure epicenter depths >1200 μm (consistent with L5b/6 border and beneath) were more likely to be associated with decreased tpO 2 at seizure onset. In addition, seizure epicenter depth and cross-laminar changes in MPC during seizure onset were also negatively correlated (Spearman's rho ¼ À 0.81, p ¼ .035, N ¼ 7, Fig. 4C right panel). Importantly, these relationships ( Fig. 4B and C) remained significant at multiple time-windows post-ictal onset (Supplemental Fig. 1). Taken together, these results suggest that the epicenter depth of seizure activity and subsequent propagatory dynamics are closely interlinked, and play an important role in shaping tissue oxygenation changes during the initiation of recurrent seizure activity.

CBF-Hbt coupling is preserved during sensory stimulation, hypercapnia and recurrent seizures
Total hemoglobin concentration (Hbt), which is proportional to cerebral blood volume (CBV) under the assumption of a constant hematocrit, also increased markedly, albeit more slowly than CBF, during recurrent seizure activity (peak Hbt: 28.8 AE 2.2%, N ¼ 7, Fig. 5A top panel, 5B). Hbt responses to stimulation and graded hypercapnia also displayed archetypal features, with slower onset and return to baseline dynamics compared to CBF (maximal Hbt: 8.3 AE 0.9%, 8.3 AE 1.2%, 20.4 AE 4%, for stimulation, 5%, and 10% FICO 2 , respectively, N ¼ 5 in all cases, Fig. 5A bottom panels, 5B). Of note, no significant difference in peak Hbt was observed between 10% FICO 2 and recurrent seizure conditions (p ¼ .38, Kruskal-Wallis test with post hoc Tukey Kramer analysis, Fig. 5B), mirroring that seen with CBF (p ¼ .22, Kruskal-Wallis test with post hoc Tukey Kramer analysis, Fig. 2A), and suggesting that higher levels of hypercapnia can induce perfusion-related responses comparable to recurrent seizures.
Trial-by-trial analysis was subsequently performed to assess the relationship between CBF and Hbt within and across experimental conditions. CBF-Hbt linear correlations were observed across experimental conditions (Fig. 5C) that were highly significant in the case of stimulation and pooled hypercapnia (Pearson's r ¼ 0.66, p < .01, N ¼ 75, and r ¼ 0.72, p < .01, N ¼ 30, respectively), and approached significance with respect to recurrent seizures (Pearson's r ¼ 0.59, p ¼ .16, N ¼ 7). Importantly, linear regression slopes were not significantly different between conditions (slope difference between whisker stimulation and seizures, p ¼ .97; slope difference between whisker stimulation and hypercapnia, p ¼ .93; and slope difference between hypercapnia and seizures, p ¼ .92), suggesting that a similar relationship between CBF-Hbt is preserved across a broad continuum of brain activation.

Neural-hemodynamic coupling is preserved during sensory stimulation and recurrent seizures
We then asked whether the relationship between neural activity and perfusion related measures was maintained across normal (whisker stimulation) and aberrant (recurrent seizures) functional activation. Trial-by-trial analysis revealed a significant linear correlation between summated L4 LFP profile activity and CBF during whisker stimulation (Pearson's r ¼ 0.7, p < .01, N ¼ 75, Fig. 6A, blue), and summated LFP responses in the ictal epicenter (site of maximal ictal activation) and CBF during recurrent seizures (Pearson's r ¼ 0.81, p ¼ .01, N ¼ 7, Fig. 6A, red). Linear regression slopes did not differ significantly between whisker stimulation and recurrent seizure conditions (p ¼ .63), although the latter did exhibit a slightly steeper gradient (β ¼ 0.7 and 0.86, respectively). Correspondingly, we also observed a significant linear correlation between L4 LFP activity and Hbt during whisker stimulation (Pearson's r ¼ 0.39, p < .01, N ¼ 75, Fig. 6B, blue), and LFP responses at the site of maximal ictal activation and Hbt during recurrent seizures (Pearson's r ¼ 0.96, p < .01, N ¼ 7, Fig. 6B, red). Yet again, no significant difference was found between regression slopes during whisker stimulation and recurrent seizures (p ¼ .22), with the latter displaying once more a steeper gradient relative to the former (β ¼ 0.39 and 0.9, respectively). Taken together, these results suggest that neural-hemodynamic coupling remains remarkably robust across normal and ictal forms of functional activation.

Discussion
The principal finding of this study is that of a previously unreported link between laminar LFP synchrony and depth-related ictal hypoxia. Notably, repetitive seizures with deep cortical epicenters were associated with prolonged periods of tissue hypoxia and reduced cross-laminar LFP synchrony, despite considerable increases in perfusion and preserved neural-hemodynamic, and CBF-Hbt, coupling. Our data suggest that functional hyperemic mechanisms at the pial-arteriole network level (to which our CBF and Hbt measures are biased) continue to operate Fig. 4. Inter-relationship between tpO 2 , laminar LFP synchrony and seizure epicenter depth. A) Laminar synchrony (MPC) spanning all laminar permutations (SGsupragranular, GRgranular, IGinfragranular, dIGdeep infragranular) against mean percentage change in tpO 2 during seizure onset. Note significant linear correlations between MPC in supragranular and underlying laminae and tpO 2 (top). B) Significant linear correlation between cross laminar synchrony (all laminae) and mean percentage change in tpO 2 during seizure initiation. C) Significant inverse monotonic correlation between seizure epicenter depth and mean percentage change in tpO 2 (left) and cross laminar synchrony (right) during seizure initiation. A-C) Linear (grey) and non-linear functions (red) fitted using least squares regression. Correlation coefficients and significance (FDR corrected for multiple comparisons), and r 2 goodness-of-fit coefficient (red) for non-linear functions, given in each case.
S.S. Harris et al. NeuroImage 171 (2018) 165-175 normatively during repetitive seizures, but are insufficient to maintain normoxia during the onset of deep lying focal seizures, most likely as a result of local depth-dependent microvascular susceptibilities. These findings provide new insights into seizure-related cerebrovascular dysfunction and may be relevant to the interpretation of perfusionrelated neuroimaging signals in epilepsy.
Relationship between cortical tissue oxygenation, depth of the ictal epicenter, and laminar LFP synchrony, during the onset of recurrent seizures Layer 6 is believed to play a key role in cortical gain control and thus sensory processing, with optogenetic activation of L6 corticothalamic neurons in visual cortex found to suppress activity across all laminae Fig. 5. Hbt changes during recurrent seizures, sensory stimulation and graded hypercapnia, and Hbt-CBF coupling. A) Percentage changes in Hbt during first 2500s of seizure activity in individual animals (top) and averaged across animals during sensory stimulation and graded hypercapnia (bottom, N ¼ 5). B) Averaged peak percentage change in Hbt during whisker stimulation (N ¼ 5), 5% and 10% fraction of inspired CO 2 (FICO 2 , N ¼ 5 in each case) and recurrent seizures (N ¼ 7). Statistical comparisons made using Kruskal-Wallis tests with Tukey-Kramer correction. Individual comparisons made using 1-tailed Wilcoxon signed-ranks tests. *P < .05 and **P < .01. Error bars are SEM. C) Trial-by-trial analysis showing strong linear correlation between standardized Hbt and CBF during stimulation ('w', N ¼ 75), hypercapnia ('h', pooled across graded FICO 2 , N ¼ 30) and recurrent seizures ('s', N ¼ 7). Grey linear functions fitted using least squares regression. Correlation coefficients and significance given as insets in each case. Significant differences between all combinations of regression slopes were tested using analysis of covariance with Tukey-Kramer correction for multiple comparisons and given as insets in each case (top).  , N ¼ 7). B) Trial-by-trial analysis of relationship between standardized LFP and Hbt during sensory stimulation (blue crosses, N ¼ 75) and recurrent seizure activity (black open circles, N ¼ 7). A-B) Linear functions fitted using least squares regression with Pearson correlation coefficients and significance given as insets in each case. Significant differences between pairs of regression slopes were tested using analysis of covariance and given as insets in each case (bottom). through recruitment of translaminar fast-spiking inhibitory neurons (Olsen et al., 2012;Bortone et al., 2014;Kim et al., 2014). Recurrent seizures whose epicenter reside in this cortical layer may therefore be predisposed to being more spatially restrained by an intense overlying inhibitory surround (Prince and Wilder, 1967;Schevon et al., 2012), leading to an overall decrease in cross laminar neural synchrony and 'functional disconnection' from surrounding cortex (Warren et al., 2010). Functional isolation of seizure-generating cortex may prevent pathologically firing neurons in the ictal focus from being mollified by the surrounding network and allow these to develop hyper-intense excitatory states (Uhlhaas and Singer, 2006). It may also underpin the accompanying robust decrease in tpO 2 during the initiation of these seizures. We have previously reported the presence of smaller, more transient, 'dips' in tissue oxygenation during onset of singular 4AP discharges in the ictal focus (Bahar et al., 2006;Zhao et al., 2009), mirroring that seen during human seizures (Cooper et al., 1966;Dymond and Crandall, 1976). Importantly, it has since been demonstrated that the magnitude of this dip, in the same seizure model, is positively correlated to the horizontal distance between the seizure focus and surface arteries and arterioles (Zhang et al., 2015(Zhang et al., , 2017. This is presumably due to steep tpO 2 gradients decreasing as a function of distance from such vessels, and rendering distal tissue comparatively more susceptible to hypoxia during metabolically demanding events (Sakad zi c et al., 2010;Kasischke et al., 2011).
It is thus possible that recurrent seizures in deep laminae may also be vulnerable to similar heterogeneity in tissue oxygenation landscape in the depth dimension, given that the absence of a significant correlation between tissue oxygenation and CBF changes during seizure onset, as well as a preserved neural-hemodynamic relationship, argues against a perfusion-related cause for the observed fall in tpO 2 . Indeed, previous work in rat using open and closed skull preparations has shown that cortical tpO 2 decreases as a function of cortical depth (Feng et al., 1988). More recent work in rat has demonstrated deep cortical layers to possess longer and more heterogeneous microvascular transit times, suggestive of reduced oxygen extraction efficacy, in contrast to layer 4 transit times which were comparatively shorter and homogeneous, in keeping with optimal oxygen extraction (Jespersen and Østergaard, 2012;Merkle and Srinivasan, 2016). As a result, recurrent seizures with deep cortical epicenters may be particularly vulnerable to microregional tissue hypoxia given a comparatively low oxygen environment that is compounded by sub-optimal local oxygen extraction efficacy as well as a reduced microvasculature density (Masamoto et al., 2004). This recipe for hypoxia may therefore play a key role in shaping the deleterious effects of status epilepticus and epilepsia partialis continua, particularly in epilepsy syndromes associated with cortical laminar disruption, such as focal cortical dysplasia. Furthermore, it may play an important role in driving hypoxia-induced angiogenesis (Krock et al., 2011) and cerebrovascular remodeling that is frequently observed in epilepsy (Marchi and Lerner-Natoli, 2013).
In contrast, recurrent seizures with middle layer (L4-L5) epicenters were associated with increased cross laminar LFP synchrony, most likely as a consequence of travelling waves of ictal excitatory synaptic activity swelling into adjacent cortical laminae (Smith et al., 2016). The increased laminar propagation of such seizures is presumably supported by strong excitatory connections between layer 4 and other cortical laminae, particularly layers 2/3 and 5, and reciprocal connectivity between layer 5 and layers 2/3 (for review see Feldmeyer, 2012). Since axonal arbors of L2/3 pyramidal neurons project to neighboring cortical columns, supragranular layers are thought to play a key role in the horizontal spread of excitation during sensory processing (Petersen et al., 2003), and have been recently shown to be the first laminae to be recruited during lateral propagation of seizure-like activity (Wenzel et al., 2017). An interesting possibility, therefore, is that, relative to deep epicenter events, seizures with middle-layer epicenters may also be associated with increased lateral seizure propagation. Furthermore, given the close correspondence between synaptic activity and hemodynamic signals (Logothetis et al., 2001), and recent reports of laminar dependent hemodynamic responses (Poplawsky et al., 2015), such seizures may therefore be able to maintain a tissue oxygen surplus by propagatively activating and availing themselves of local laminar neurovascular mechanisms, and exploiting a cortical angioarchitecture that is already optimized and accustomed to the elevated metabolic needs of the cortical input layers; such as increased microvascular (Masamoto et al., 2004) and arteriole primary branch density (Blinder et al., 2013), and red blood cell content and flow changes (Srinivasan and Radhakrishnan, 2014).
Our finding of a relationship between laminar field potential synchrony and tissue oxygenation is in keeping with a burgeoning body of work showing distributed patterns of cerebral hemodynamics to be spatiotemporally coupled to synchronized neural activity in the resting state (Ma et al., 2016) and during inter-ictal events (Khoo et al., 2017), and that neural synchrony, not activity, shapes the relationship between LFPs and the BOLD signal (Hermes et al., 2017). Furthermore, our observation of both decreases and increases in cross-laminar LFP synchrony complements emerging work that has begun to challenge the traditional view that seizures are solely 'hypersynchronous' events (Penfield and Jasper, 1954), and suggests that decreases in neural synchrony during ictal evolution may be, paradoxically, crucial for initiation and maintenance of seizure activity (for review see Uhlhaas and Singer, 2006). While it is as yet unclear why LFP synchrony between supragranular and underlying laminae is preferentially correlated to changes in tpO 2 during the onset of recurrent seizures, it is interesting to note that GABAergic inhibitory neurons densely populate the superficial cortex (Meyer et al., 2011) and could play a key role in local neurovascular communication in the upper layers were it to be present. Our future work will address this directly through multi-depth electrophysiology and tpO 2 recordings, and multi-laminar optogenetic activation of neuronal populations.

Preserved neural-hemodynamic and CBF-Hbt coupling
Our results have also shown that the relationship between neural activity and perfusion-based measures (i.e. CBF and Hbt) in whisker barrel cortex is linear, preserved, and remarkably robust during functional activation by sensory stimulation and intense recurrent seizures. This finding is significant as it provides support for the assumption of linearity in neural-hemodynamic coupling that is routinely made when interpreting perfusion-based neuroimaging signals, such as BOLD fMRI, in terms of underlying neural activity (Logothetis et al., 2001). This is additionally relevant given the rise of simultaneous clinical EEG-fMRI studies to localize hemodynamic correlates of electrophysiological epileptic activity during pre-surgical work-ups, and which are being increasingly conducted during ictal, rather than interictal periods, since this may provide a more accurate identification of the epileptogenic zone (Salek-Haddadi et al., 2002;Tyvaert et al., 2008). However, it is important to note that the observed linear regression slopes, while providing an excellent approximation, were somewhat steeper in the recurrent seizure condition compared to stimulation. This suggests that neural-CBF/Hbt coupling may be intrinsically non-linear at pathological extremes, a trait that may be more discernable during dynamic changes at seizure onset and offset (Ma et al., 2012;Harris et al., 2014a).
Similarly, we also found that that the relationship between CBF and Hbt (i.e. CBV) across all experimental conditions was well approximated by a linear function that was comparable along a continuum of focal to global hemodynamic changes. These similarities indicate that the flowvolume relationship is consistent irrespective of the mechanisms by which changes in flow are induced, and are of relevance to the calibration of the BOLD contrast in quantitative fMRI (Davis et al., 1998;Kida et al., 2007). This finding confirms and extends our previous work showing equivalent Grubb flow-volume coefficients during sensory stimulation and hypercapnia (Jones et al., 2001), and that of linearity between flow and volume during initiation and evolution of singular ictal-like discharges using non-concurrent techniques (Zhao et al., 2009, ictal discharge duration: 96.9 AE 3.1s, peak increase in CBF:55.6 AE 7.8%, peak increase in Hbt: 7.2 AE 1.5%). Notwithstanding, magnitude-related models do not capture temporal dynamics of flow and volume changes, with previous work confirming that Hbt/CBV onsets and returns to baseline more slowly than CBF (Figs. 1 and 5; Jones et al., 2001;Kida et al., 2007), due to factors such as the viscoelastic nature of blood vessels (Mandeville et al., 1999) and/or capillary hyperemia (Hillman et al., 2007). Further studies using dynamic biophysical models are therefore needed to confirm that temporal flow-volume differences are not exacerbated under pathological conditions.

Conclusion
In summary, we have demonstrated, using a novel multi-modal approach, that laminar patterns of seizure initiation and evolution, evinced by translaminar synchrony, play an important role in shaping ictal cerebrovascular responses, and that these can be associated with prolonged periods of tissue hypoxia despite profound increases in perfusion and conserved neural-hemodynamic coupling. These findings provide new insights into the neurophysiological correlates of seizurerelated cerebrovascular responses, and reveal a hitherto unknown differential laminar susceptibility to tissue hypoxia during ictal events. It is thus tempting to speculate that laminar-specific predisposition to cerebrovascular dysfunction may play an important role in shaping the deleterious effects of recurrent seizures in epilepsy syndromes associated with laminar dysfunction, such as focal cortical dysplasia. However, future research, perhaps using simultaneous intracranial depth EEG and neuroimaging, is needed to confirm these findings in the clinical setting, and determine whether the observed mechanisms can contribute to neurocognitive impairment in chronic epilepsy.

Funding
This work was supported by the Medical Research Council (grant number MR/M013553/1); and Epilepsy Research United Kingdom (grant number P1501).