Ben Ridley
email: ben.ridley@univ-amu.fr

Jonathan Wirsich
email: jonathan.wirsich@univ-amu.fr

Gaelle Bettus

Roman Rodionov
email: r.rodionov@ucl.ac.uk

Teresa Murta
email: teresa.murta.11@ucl.ac.uk

Umair Chaudhary
email: umair.chaudhary@ucl.ac.uk

David Carmichael
email: d.carmichael@ucl.ac.uk

Rachel Thornton
email: rachel.thornton@gosh.nhs.uk

Serge Vulliemoz
email: serge.vulliemoz@hcuge.ch

Mcevoy
email: a.mcevoy@ucl.ac.uk

Fabrice Wendling
email: fabrice.wendling@univ-rennes1.fr

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

Jean-Philippe Ranjeva

Louis Lemieux
email: lemieux@ucl.ac.uk

Maxime Guye
email: maxime.guye@ap-hm.fr

Simultaneous intracranial EEG-fMRI shows inter-modality correlation in time-resolved 1 connectivity within normal areas but not within epileptic regions

Keywords: connectivity, multimodal imaging, resting-state, focal epilepsy, dynamic connectivity 62 63 ECoG -Electrocorticography 65 FC -Functional Connectivity 66 Lobe Epilepsy 67 IED -Interictal Epileptic Discharges 68 icEEG -Intracranial Electroencephalography 69 ICN -Intrinsic Connectivity Network

For the first time in research in humans, we used simultaneous icEEG-fMRI to examine the link 41 between connectivity in haemodynamic signals during the resting-state (rs) and connectivity derived 42 from electrophysiological activity in terms of the inter-modal connectivity correlation (IMCC). We 43 quantified IMCC in 9 patients with drug-resistant epilepsy i) within brain networks in 'healthy' non-44 involved cortical zones (NIZ) and ii) within brain networks involved in generating seizures and 45 interictal spikes (IZ1) or solely spikes (IZ2). Functional connectivity (h 2 ) estimates for 10 minutes of 46 resting-state data were obtained between each pair of electrodes within each clinical zone for both 47 icEEG and fMRI. A sliding window approach allowed us to quantify the variability over time of h 2 48 (vh 2 ) as an indicator of connectivity dynamics. We observe significant positive IMCC for h 2 and vh 2 , 49 for multiple bands in the NIZ only, with the strongest effect in the lower icEEG frequencies.

50

Similarly, intra-modal h 2 and vh 2 were found to be differently modified as a function of different 51 epileptic processes: compared to NIZ, h 2 BOLD was higher in IZ1, but lower in IZ2, while h 2 icEEG 52

showed the inverse pattern. This corroborates previous observations of inter-modal connectivity 53 discrepancies in pathological cortices, while providing the first direct invasive and simultaneous 54 comparison in humans. We also studied time-resolved FC variability multimodally for the first time,

55

finding that IZ1 shows both elevated internal h 2 BOLD and less rich dynamical variability, suggesting 56 that its chronic role in epileptogenesis may be linked to greater homogeneity in self-sustaining 57 pathological oscillatory states.

INTRODUCTION 76

The advent of multimodal approaches for investigating brain activity, and the range of and Carmichael 2014; Holmes et al. 2014;Ridley et al. 2015). While region of interest (ROI) 111 selection in such studies may be motivated by an understanding of the likely candidates for 112 involvement in epileptogenesis [cf. [START_REF] Bettus | Decreased basal fMRI functional connectivity in 611 epileptogenic networks and contralateral compensatory mechanisms[END_REF], they may still fail to account for individual 113 variation in the actual location and extent of pathological activity. A small number of studies have 114 taken advantage of the greater specificity of icEEG relative to a priori regions of interest when 115 parcellating cortices into irritative and physiological/non-involved zones (NIZ) to explore functional 116 connectivity specifically as a function of involvement in different epileptic processes in the interictal 117 resting-state (Bartolomei et al., 2013;Bartolomei et al., 2013;[START_REF] Bettus | Enhanced EEG functional connectivity in mesial 617 temporal lobe epilepsy[END_REF]Schevon et al., 118 2007;Warren et al., 2010). Such studies indicate increased functional connectivity as measured by 119 icEEG (icEEG-FC) in regions involved in interictal epileptiform activity in comparison to spared 120 cortices.

121

To our knowledge, [START_REF] Bettus | Interictal functional connectivity of human epileptic 614 networks assessed by intracerebral EEG and BOLD signal fluctuations[END_REF] is the only study to have used the specificity of the 122 icEEG gold standard in parallel with BOLD-derived FC estimates from the same location, within the 123 same individuals. These authors showed an apparent inter-modality discrepancy: higher FC as 124 revealed by icEEG and lower FC as measured by BOLD in electrophysiologically abnormal regions.

125

However, since both modalities were not acquired simultaneously, confounding inter-session effects 126 cannot be ruled out as an explanation. For example, impacts arising from the time of day of Analyzer software (V1.3). Data was analysed after subtracting the common average reference 230 (minimum number of electrodes included: 30, P01), and filtered according to several frequency ranges 231 of interest: broadband (0.5-100Hz), Delta (0.5-3.4 Hz); Theta (3.4-7.4 Hz); Alpha (7.4-12.4 Hz); Beta 232 ; and Gamma (24-100 Hz).

233

While every recording electrode was used to define an ROI, the icEEG and fMRI data were 234 subject to visual inspection by experts [GB and FB] and overly noisy and likely artifactual timeseries 235 of either modality were excluded from all analyses. Using this approach, 13 leads/ROIs were rejected

236

(6 grid electrodes from P03 and one depth from P06 due to icEEG data, and 6 grid electrodes from

237

P06 due to BOLD data) and 570 were retained.

341

Pathological regions also demonstrated relationships not found in unaffected regions. In 342 particular, significant aberrant negative correlations were found between h 2 icEEG and vh 2 BOLD in all 343 icEEG bands, with the effects being strongest at the highest frequencies.

Finally, IEDs were found to be positively correlated with vh 2 BOLD in IZ2, while correlating

345 with h 2 BOLD in IZ1.
346 347

Impact of zone on unimodal h 2 and vh 2 348

While accounting for confounding covariates, significant main effects (Bonferroni-corrected)

349

of zone were observed on both h 2 and vh 2 for BOLD and across all sub-bands for icEEG-FC (Table 350 4). Post-hoc t-tests indicated a differential impact of epileptic processes on FC as well as disjunction

351

in this impact depending on modality (Tables 56

361

The current work confirms a previously observed discrepancy between icEEG-and BOLD-

362

derived FC findings (h 2 ) [START_REF] Bettus | Enhanced EEG functional connectivity in mesial 617 temporal lobe epilepsy[END_REF][START_REF] Bettus | Interictal functional connectivity of human epileptic 614 networks assessed by intracerebral EEG and BOLD signal fluctuations[END_REF]) and extends it to FC variability over 363 time (vh 2 ); in addition, we shed new light on the differential impact of different epileptic processes.

364

We demonstrate for the first time using simultaneously acquired icEEG-fMRI that epileptic cortices 365 are distinguished by their internal multimodal resting-state FC as a function of their involvement in 366 different, clinically-relevant epileptic processes. Furthermore, we provide the first evidence that the 367 inter-modal connectivity correlation that exists in unaffected regions may be modified in cortices

368

affected by epileptic pathology, when estimating h 2 and vh 2 in unaffected regions and using invasive

436

We note that the variability in findings observed across the literature in 'epileptogenic'

437

regions is reflected in the inter-individual variability of FC estimates (Figure 3, left) and differences

438

between IZ1 and non-involved zones in our sample (Figure 3

515

Our findings indicate that in healthy cortex BOLD h 2 is related to vh 2 icEEG across sub-bands,

516

but that BOLD vh 2 shows a relationship to icEEG vh 2 only in alpha, beta and gamma bands (Table 3).

517

Both findings may reflect the imbalance when trying to predict FC in one modality from the other

532

This could maintain the brain in a state of heightened competition ready to evolve momentarily into a 533 specific state under the influence of internal or external sensory modulation (Deco and Corbetta 2011). Our results show that IZ1 is associated with both an augmentation of h 2 BOLD and a reduction in vh 2 BOLD. This could reflect a dynamic repertoire that is reduced in diseased cortices under the 536 influence of both epileptogenic processes and paroxysmal interictal activity, with a propensity toward

Figure 1

 1 Figure 1 depicts the position of electrodes and their labelling for each participant. For a detailed

,

  Figure 2, top). In particular, IZ1 demonstrated 352 higher h 2 BOLD, and IZ2 lower h 2 BOLD compared to NIZ, while for h 2 icEEG the inverse pattern was found 353 across most bands. Across modalities, regions responsible for generating seizures (IZ1) exhibited the 354 most variance in terms of their respective level of augmentation or disruption relative to the other 355 zones (Figure 3, left & right) and IZ2 exhibited a more consistent relationship with the other zones 356 (Figure 3, middle) of elevated h 2 icEEG and reduced h 2 BOLD. This discrepancy between modalities was 357 also reflected in their variability over the sliding window time series (Figure 2, bottom), with 358 pathological regions exhibiting elevated vh 2 icEEG, and lower vh

  Additionally, networks thought to reflect intrinsic pattern of connectivity (ICNs), are also 406 subject to a variety of modifications in epilepsy (Centeno and Carmichael 2014). An apparent 407 dissociation between interictal electrophysiological and haemodynamic metrics of connectivity in 408 focal epileptic regions was indicated by Bettus and colleagues (2011) when analysing non-409 simultaneous data acquired from the same individuals. They found evidence for a diminution of 410 functional connectivity within the IZ2 as measured by BOLD, and an augmentation in h 2 icEEG within 411 IZ1 and IZ2 which was significant in the beta band. Here, using simultaneously-acquired icEEG-412 fMRI, we confirm this inter-modality discrepancy in connectivity estimates (Figure 2, top), as well as 413 extending it to the dynamic variability of FC estimates over time (Figure 2, bottom). In fMRI 414 studies investigating connectivity in epilepsy, where 'involved' regions are selected on the a priori 415 basis of being likely candidates for involvement in epileptogenesis, a common finding is a reduction 416 in BOLD-FC (Centeno and Carmichael 2014) . To our knowledge, Bettus et al. (2011) is the only 417 study to have used the specificity of the icEEG gold standard in parallel with BOLD-derived FC 418 estimates from the same location within the same individuals, and observe a trend to BOLD-FC419reduction in the IZ1 versus NIZ, though it is non-significant. Interestingly, while an inter-modality 420 discrepancy is also observed in the current results for IZ1, it points in an opposite direction than these 421 limited antecedents might lead one to expect: an augmentation of BOLD-FC and reduction in icEEG-

493

  multiple processes that contribute to neurovascular coupling and hence the BOLD-signal generation 494 cascade are effected in epilepsy, and as such represent sources of unmodeled noise standing between 495 underlying processes and the signals as estimated via such metrics as (de-)activation or connectivity.496The separate contributions of perfusion deficits (Duncan 2009) and altered metabolic demands in the 497 chronic epileptic state and under the specific influence of epileptic transients will likely need to be 498 characterized. Multimodal exploration as traditionally performed and in new combinations is one tool to mutually constrain different sources of non-neuronal noise and mutually support conclusions by 500 filling in missing information (Uludağ and Roebroeck 2014). establishing a relationship between indices of connectivity in simultaneous 504 icEEG and fMRI data, sliding window analyses allowed us to resolve the variability of connectivity 505 over time. Comparable analyses have been applied previously to fMRI data in temporal lobe epilepsy 506 patients, showing modifications in time-resolved connectivity associated with disease burden of 507 memory scores (Douw et al. 2015), and both increases (Laufs et al. 2014) and decreases (Nedic et al. 508 2015) in the variability of FC in regions that might plausibly be considered part of epileptic networks 509 on an a priori basis. Current results, in regions assigned to IZ1 and IZ2 by the electrophysiological 510 gold standard, tend to agree with the latter: vh 2 BOLD is significantly reduced in IZ1 (Figure 2). 511 Furthermore, we provide the first simultaneous icEEG-fMRI evidence in humans of an inter-modality 512 discrepancy in the variability of FC over time in epileptic regions: in opposition to BOLD, vh 2 513 estimates derived from icEEG indicate an augmentation across bands in irritative regions, and 514 specifically in IZ1 in alpha, beta and gamma bands.

518(

  Deligianni et al. 2014), with variability in all icEEG bands imparting information about BOLD 519 connectivity oscillations, but the inverse not being the case. This is in good agreement with previous 520 work indicating covariation in correlation between modalities in comparable bands in rat and macaque 521 studies(Magri et al. 2012; Thompson et al. 2013), and possibly reflects the widely hypothesized role 522 of alpha and gamma bands in particular in communication between neural assemblies at distal and 523 local scales, respectively (Keilholz 2014). As with multimodal relationships between estimates of 524 connectivity, multimodal relationships in the variability in FC across sliding window timeseries 525 appear to be disrupted in pathological cortices in epilepsy.526Whole brain computational simulation approaches comparing empirical data to simulated 527 functional data derived from structural models suggest that the agreement between simulated and 528 empirical data are best when both a low energy spontaneous state and several states of localized high 529 activity are stable states of the system(Cabral et al. 2014). In this context, slow resting-state 530 oscillations represent the dynamic exploration of the different states of the brain's intrinsic functional 531 repertoire instigated by underlying anatomic connectivity and intrinsic noise (Deco and Jirsa 2012).

Fig. 3

 3 Fig.2 Post-hoc t-contrasts between zones in terms of functional connectivity (top) and variability

Table 2

 2 provides a breakdown of the number of ROIs 238 provided for the final analysis by each patient.

	239
	240
	2.4

Functional connectivity analysis 241

  A nonlinear measure of covariance, h 2 , was used to estimate functional connectivity between 242 every pair of electrodes. SeeWendling et al. (2010) for an in depth account of this metric. Briefly, 243 dependency between two signals X and Y derived from the same modality was quantified using:

	263	h 2 values between pairs of electrodes or ROIs were computed over sliding window analyses,
	264	separately for each modality. When computing h 2 icEEG between electrodes a 5s window moving with
	265	0.5s increments was used, while for computing h 2 BOLD between ROIs a 90s window moving with 2s
	266	increments was used. Window size was motivated by the need to be large enough to include
	267	uncorrelated (in time) X-Y values in order to compute a meaningful correlation, and on the assumption
	268	that 50 uncorrelated couples of X-Y values are needed for reliably computing the h 2 coefficient.
	269	Taking the auto-correlation function for EEG to be 100ms (Wendling et al. 2001), 50 samples
	270	provides a 5 second window size. 1.8 seconds was taken to be sufficient to include minimal
	271	autocorrelation given evidence that the autocorrelation coefficient of resting-state BOLD oscillations
	272	(TR=3s) in healthy controls drops of sharply after one second and reaches zero on the order of two
		seconds, and 50 samples yields a 90 second window (Kaneoke et al. 2012).
	276
	244 277	run. An estimate of h 2 and vh 2 was obtained for BOLD, and for broadband and five sub-bands for
	245 278	icEEG signals. A thresholding procedure outlined in Supporting Information 2 was used to verify
	279 246 280 281	ℎ 2 (𝜏) = 1 -that all pairwise relationships constituted evidence for the existence of a correlation between variables 𝑉𝐴𝑅[ 𝑌(𝑡 + 𝜏) ] 𝑋(𝑡) of interest at the 99% level of confidence. 𝑉𝐴𝑅[𝑌(𝑡 + 𝜏)] In order to take advantage of a pre-existing module of the AMADEUS software (Wendling et
	247 282 283	With al. 2010; Wendling 2015) -designed to compute large sets of pairwise h 2 -for our BOLD data and to 𝑉𝐴𝑅 [ 𝑌(𝑡 + 𝜏) 𝑋(𝑡) ] ≙ arg 𝑓 min ensure equivalency of processing it was necessary to resample BOLD data at 250 Hz using a 1-
	284	dimensional 1st order linear interpolation routine available in Matlab (interp1). This is identical to the
	285	procedure employed by Bettus et al. (2011), who also demonstrated the comparability of h 2 estimates
		on highly sampled and normal data.
	257
	258	2010). Note that h 2 is a measure of amplitude/power covariance and indicates the presence and
		strength of relationships (and does not, per se, differentiate between 'correlations' and

ℎ (𝐸[𝑌(𝑡 + 𝜏) -ℎ(𝑋(𝑡))]

2 )

248

249

where h is a nonlinear fitting curve which approximates the relationship between X and Y and t is the 250 time shift that maximizes the value of the h 2 coefficient. 251 Conceptually, h 2 is a normalized non-linear correlation coefficient which quantifies the 252 reduction of variance in Y when X is used to predict Y samples. When h 2 XY=1, Y is fully determined by 253 X, and h 2 XY=0 when no relationship exists between the two signals. While it is known that non-254 linearities can occur in EEG signals especially in epilepsy (Casdagli et al. 1997; van Diessen et al. 255 2015), the relative importance of being able to detect linear and non-linear components of 256 relationships between signals is a subject of ongoing discussion (Netoff et al. 2004; Wendling et al.

273

The h 2 values were averaged across windows over the entire scanning period time in order to 274 get a single estimate of static functional connectivity (h 2 ). Additionally, the variation of h 2 over time

275

(vh 2 ) was computed as the standard deviation of the h 2 estimates for each window across the entire

286 287 2.5 Other metrics: IED count and inter-electrode distance 288

	298 329	Euclidean distances (mm) were calculated between the centres of each ROI as a proxy for 3. RESULTS
	299	inter-electrode distances in situ. 3.1
	300
	301	2.6 Statistical comparisons
	302	The following statistical analyses were performed using JMP version 9 (SAS Institute Inc.,
	303	Cary, NC) on pairwise interactions between electrode/ROIs within the same clinical zone (IZ1, IZ2,
	304	NIZ). Relationships between pairs in different zones were not considered. Table 2 provides
	305	information regarding the number of pairwise interactions provided within each zone, yielding an
	306	overall sample of n=11543.
	307	We were also interested in the extent to which connectivity estimates derived from within
	308	each modality are associated between modalities. We assessed the inter-modal connectivity
	309	correlation (IMCC) via Pearson's r between the sample of pairwise connectivity estimates obtained
	310	from BOLD ROIs (h 2 BOLD) on the one hand, and the sample of pairwise connectivity estimates derived
	311	from icEEG electrodes (h 2	icEEG; broadband plus five sub-bands) on the other. The same procedure was
	312	applied between vh 2 BOLD and vh 2 icEEG. Correlation of h 2 BOLD and vh 2 BOLD in relation to IEDs was also
	313	considered. Thus IMCC was investigated for each of the three clinically-defined zones, and
	314	considered significant at a Bonferroni-corrected level of 0.05/((6+1)*3)=p<0.002.
	315	Differential impacts of clinically-defined zones on mean h 2 and vh 2 was modelled (ANOVA)
	316	using 'zone' as regressor of interest (3 levels) while controlling for several sources of 'spurious'
	289	IEDs counts were calculated automatically with AMADEUS following Bourien et al (2005).
	290 325	The approach involves detecting abrupt increases (high amplitude, short duration spikes) in the mean
	291 326	value of the squared modulus of the signal when passed through a wavelet filter bank, enhancing the
	292 327	fast sharp component of the IED relative to the surrounding background EEG. A single IED count
	293 294 328	was associated with each pairwise interaction between electrodes by taking the average of the individual electrode spike counts. It should be noted that this automated procedure does not
	295	differentiate between paroxysmal activity of different types, and these automatic estimates of IEDs
	296	are distinct from those obtained manually during extended clinical monitoring and used to define
		clinical zones.

variation. Two categorical covariates were included: 'Patient' (9 levels) and 'Electrode Type' (2 317 levels), the latter reflecting a departure with

[START_REF] Bettus | Enhanced EEG functional connectivity in mesial 617 temporal lobe epilepsy[END_REF][START_REF] Bettus | Interictal functional connectivity of human epileptic 614 networks assessed by intracerebral EEG and BOLD signal fluctuations[END_REF]

in which patients were 318 implanted solely with depth electrodes as in P01, P04 and P09 here, but unlike P06 implanted solely 319 with grids and all other patients who had a mixture of depth and grid electrodes. Two continuous 320 covariates, 'IEDs' and 'Euclidean distance' were also included, to account for other possible sources 321 of variability inherent in the spatial sampling and incidence of paroxysmal activity in epileptic 322 networks. Tests were compared to two thresholds for significance: an exploratory level of p<0.05, and 323 on the assumption of multimodal associations a Bonferroni correction of the number of modalities 324 (BOLD plus six icEEG bands) threshold of 0.05/7=p<0.007.

Inter-modal connectivity correlations for h 2 and vh 2 330Table 3

 3 lists the observed correlation coefficients and associated p-values. In non-involved 331 (NIZ) cortices, h 2 BOLD was significantly positively associated with both h 2 and vh 2 of all frequency 332 bands measured by our icEEG data. IMCC was strongest in the slower icEEG bands, and weakest at

	333
		the highest frequencies studied here. h 2 BOLD was negatively correlated with vh 2	BOLD. A negative
	335	association was also exhibited between vh 2 BOLD and vh 2	icEEG in the alpha, beta and gamma bands.
	336	IMCC as observed in the non-involved region was not found in pathological cortices: no
	337	significant relationships were observed between h 2 BOLD and h 2 icEEG outside of the NIZ with the
	338	exception of the delta band in IZ1. Similarly, the only relationship between vh 2 icEEG and vh 2	BOLD in
	339	epileptic cortices was a positive correlation observed in the alpha band in IZ2. Finally, while negative
		correlations between h 2 BOLD and vh 2

BOLD were found in all zones, this negative relationship was 340 stronger in pathological regions, substantially so for IZ1.

Inter-modal connectivity correlation in physiological and pathological cortices 372

  In contrast, we find almost no evidence of the inter-modality correlation observed in NIZ for 396 h 2 and vh 2 estimates in IZ1 and IZ2. One possibility is that each modality captures a different aspect 397 of the functional reorganization of networks in the context of epilepsy[START_REF] Bettus | Role of resting state functional connectivity 608 MRI in presurgical investigation of mesial temporal lobe epilepsy[END_REF] Bettus et al. 398 2011; Duncan et al. 2013; Bartolomei et al. 2013b). In this scheme, slow (<0.1Hz) BOLD signal 399 fluctuations may reflect perturbance of the functional integrity of macro-level intrinsic networks, and 400 the faster and broader frequency range of EEG to abnormal organisation of epileptic networks poised 401 to evolve into hypersynchrony at the onset of seizures(Bartolomei et al. 2013a; Bartolomei et al. Bettus et al. 2008; Bartolomei et al. 2013a; Bartolomei et al. 2013b) and fMRI(Bettus et al.

	402
	403	2013b) or interictal activity (Centeno and Carmichael 2014). Modification of FC within regions
	369 404	characterised by involvement in specifically epileptic processes is supported by icEEG (Schevon et al.
	370	and simultaneous electrophysiological and haemodynamic signals. 2007;
	371
	4.1 390
	391	recording and simultaneous fMRI acquisition, we observe correlations between FC estimated on
	392	haemodynamic and electrophysiological signals across frequency bands in non-involved, 'healthy'
	393	cortex (Table 3). We find the strongest relationships at lower frequencies, possibly reflecting a bias
	394	toward lower bands in maintaining relationships over distant regions due to the relationship between

The non-invasive BOLD signal is an indirect measure of neurological activity, motivating 373 interest in better understanding its electrophysiological correlates, including in terms of icEEG-374 derived connectivity and its relationship to both spontaneous BOLD activity and resting-state 375 functional connectivity. Correlates from electrophysiological activity (as opposed to connectivity) at 376 multiple spatial and temporal scales have been proposed (Keilholz 2014; Tagliazucchi and Laufs 377 2015), with both high-frequency local field potentials (Shmuel and Leopold 2008; Nir et al. 2008; 378 Schölvinck et al. 2010), and lower frequencies (Lu et al. 2007; He et al. 2008; Pan et al. 2013; Lu et 379 al. 2016) being advanced as electrophysiological correlates of spontaneous rsfMRI fluctuations. 380 Others have examined the extent to which connectivity estimates derived from EEG and BOLD are 381 associated with one another, in terms of the intermodal correlation of the connectivity estimates 382 derived from each modality. Invasive electrophysiological recordings in rats indicate relationships 383 between BOLD correlations and EEG band power correlations that are strongest in lower frequency 384 bands, especially delta, but extending to gamma (Lu et al. 2007; Pan et al. 2011). Modulation by 385 exogenous stimulation of rsfMRI functional connectivity in rat whisker barrel cortex has been found 386 to be reflected in the delta range, but not at higher frequencies (Lu et al. 2016). In humans, it has been 387 shown that BOLD covariance (functional connectivity) matrices can be well explained based on EEG 388 covariance matrices -similarly so across bands -while the reverse is only true in the lower frequency 389 bands (Deligianni et al. 2014). Here, in the first study to combine invasive electrophysiological distance and signal delay (Schölvinck et al. 2013).

  While studies that examine BOLD-FC between ICN nodes typically report 423 decreases in patients, modifications within nodes tell a more complicated story. Not only decreases 424 but increases in connectivity are evident in the frontal and temporal nodes of multiple ICNs, both for

	431	connectivity increases reported when it is specifically the vicinity of individually localized
	432	epileptogenic regions under consideration. Stufflebeam et al. (2011) found BOLD-FC increases in
	433	local (<5mm) connectivity in proximity to grid electrodes classified as belonging to the seizure onset
	434	zone in 5 of 6 patients with focal epilepsy. Similarly, Luo et al. (2014) reported an enhancement of
	435	FC in FLE patients in the immediate neighbourhood of cortical epileptogenic zones defined by EEG-
		fMRI data fusion.
	425
	426	TLE (He et al. 2008; Liao et al. 2010) and FLE (Braakman et al. 2013; Widjaja et al. 2013). Though
		increases have been interpreted as potentially compensatory (Bettus et al. 2009; Bettus et al. 2010;

422

FC (Figure2, top).

, bottom and Table 2),

	466	Haylett et al. 2013), but that evoked neurovascular coupling is modified during ictal epileptic activity
	439	, right). Insofar as modifications to
	440	connectivity in epileptic networks which are salient to different modalities might reflect different
	441	disease process, as suggested above, the extent to which these might interact with different aetiologies
	442	is unknown, and is an area that should receive further attention. Indeed, a recent analysis of
	443	magnetoencephalography data in both mesial TLE and focal neocortical patients indicates a mixture
	444	of focal increases and decreases in pre-surgery FC in epileptogenic regions that would later undergo
	445	resection (putative IZ1) versus non-involved homologous regions in the contralateral hemisphere
	454
	455	versus the positive correlation observed between IEDs and vh 2 BOLD and the anti-correlation of vh 2 BOLD
	456	and h 2 icEEG (Table 3).
	457
	458	4.2 Altered neurovascular coupling in pathological cortices
	459	Altered neurovascular coupling has been proposed as an explanation of apparent
	460	discrepancies between icEEG-and BOLD-derived indices of FC (Bettus et al. 2011), and might also
	461	play a role in the absence of correlation seen here. Metabolic anomalies in electrophysiologically
	462	abnormal regions have been demonstrated in both TLE (Guye et al. 2002) and FLE (Guye et al.
		2005), with the addition of both vascularisation defects (Duncan 2010) and blood-brain barrier

(Englot et al. 2015)

. In contrast, the IZ2 may represent a more similar situation across individuals than 446 IZ1 (Figure

3

, IZ2-NIZ contrast), as relatively 'healthy' cortex extending beyond focal lesions that is 447 not so degenerate as to self-generate seizures, and where IEDs are the main form of pathological 448 activity impinging on brain networks. For example, while metabolic abnormalities suggestive of 449 defective neurovascular coupling are observed in IZ2 (Guye et al. 2002; Guye et al. 2005), they are 450 reversible as seen in post-operative patients who have become seizure free after successful resection 451 surgery (Serles et al. 2001). This difference in the relationship to ongoing disease processes may be 452 reflected in the fact that IEDs appear to be disruptive to h 2 BOLD connectivity leading to its reduction in 453 IZ2 (Figure 2, top) and greater variability as measured by vh 2 icEEG (Figure 2

(

  Harris et al. 2013). This is suggestive of a broad comparability between spontaneous and evoked 467 neurovascular coupling which may be subject to modification in the context of epilepsy.

	468
	469	Our results suggest that multimodal differences persist despite an equivalent incidence of
	470	inter-modal between-session factors (as ensured by simultaneous acquisition), though this in itself
	471	does not exclude a difference in impact between modalities. Modelling may be one means of
	472	understanding how such differences contribute to the discrepancies seen here, but it should be noted
	473	that the choice of parameters derived from one signal for the understating of another is highly non-
	474	trivial. For example, an argument for a disruptive effect of interictal epileptic transients such as IEDs
	475	on connectivity has a prima facie plausibility. However, empirical results have been more equivocal
	476	with negative (Nissen et al. 2016), positive (Bartolomei et al. 2013a) and null (Bettus et al. 2008)

associations of IEDs with FC in epileptics networks observed in scalp and invasive electrophysiology, 477 but also overall limited-to-no-effect of removing IED-containing epochs by censoring data

(Bettus et

with resting-state phenomena in the form of IED-and seizure-associated BOLD changes in regions 481 associated with intrinsic connectivity networks, most especially the default mode network (Gotman 482 and Pittau 2011; van Graan et al. 2015). Interictal phenomena may be associated with increases, 483 decreases or unchanged BOLD signal (Bénar et al., 2006; Salek-Haddadi et al., 2006; Stefan and 484 Lopes da Silva, 2013; Thornton et al., 2011), and BOLD changes may be found in regions with no 485 apparent involvement on EEG (Kobayashi et al. 2006). Interpreting the variability of these results is 486 further complicated by recent evidence that coupling is state-and time-dependant in addition to being 487 spatially organized in the resting-state (Feige et al. 2016). While suggestive of a framework for 488 understanding the differential impact between modalities, how altered coupling might pertain to the 489 findings in the current study -the focus of which is FC in epileptic networks rather than resting 490 activity in ICNs per se -is unclear. Indeed, there is evidence of dissociation between different means 491 of characterizing resting state activity in epilepsy, with discordant finding in terms of the amplitude 492 and connectivity of resting state activity in the same patients (Zhang et al. 2015). As noted above,

  . A validation study byCarmichael and colleagues (Carmichael et al. 2012) found that signal reduction was similar for depth and grid electrodes, and was on average less than 570 50% at 5mm from the electrode, the radius of the ROIs used here. Similarly, data from simultaneous 571 EEG-fMRI studies suggests minimal disturbances of temporal signal to noise ratios, with comparable 572 detectability of activation fMRI time-courses in the presence and absence of EEG recording We provide the first evidence derived from intracranial EEG and simultaneous BOLD signals 588 of inter-modality correlation in healthy human cortex in terms of both static functional connectivity 589 and its time-resolved variability. Furthermore, while ruling out differences in the intra-individual Schematic representation of individual implantation schemes 827 projected onto templates in MNI space. Red spheres: ROIs in the primary irritative zone (IZ1); yellow 828 spheres: ROIs in secondary irritative zone (IZ2); black spheres: ROIs in non-involved zone (NIZ).

	587 826	5. Conclusion Lopes R, Moeller F, Besson P, et al (2014) Study on the Relationships between Intrinsic Functional Figure Legends
		Fig.1 Regions of interest (ROIs)
	537
	seizures being one of the states that such influences promote. If a reduced repertoire reflects homogeneity in self-sustaining pathological oscillatory states, this could be reflected in both abnormally high FC and less rich dynamics, consistent with the negative relationship between h 2 and vh 2 of BOLD which we find to be augmented in irritative cortices and at its highest in IZ1. In contrast, IZ2 may represent cortex that is less degenerate in this sense -showing no reduction in vh 2 BOLD compared to NIZ -but is chronically subject to disruption by interictal transients which could be reflected in the correlation of IEDs and vh 2 BOLD (Table 3). 4.4. Limitations and Future Directions The implantation, placement and coverage of intracranial electrodes are motivated first and foremost by clinical necessity, leading to a degree of heterogeneity in our sample. We have attempted to identify and account for several forms of heterogeneity via quantitative proxies including Euclidean distance and automatically-generated IED counts. For example, since the automatic extraction of IEDs is based on locating abrupt and transient signal modulations there is a possibility of misidentifying non-pathological transients as IEDs, as is also the case for manually defined IEDs (Noachtar and Rémi 2009). Supplementary Figure 1 shows the number of 'IEDs' identified by zone at group and individual levels, and though the relative amounts of 'IEDs' detected by the algorithm meet clinical expectations (Bartolomei et al. 2016), a small amount were located in regions classified as NIZ based on experienced clinical analysis. It is likely that these proxies account for only some of the variance involved in estimating connectivity in epilepsy, and more exact measures of propagation distance and a more discerning treatment of different types of interictal electrophysiological phenomena may provide further insights. Additionally, a development of the current work could involve the normalization of BOLD-FC against estimates derived from controls using the ROIs extracted from the electrodes in each patient, to allow the demarcation of differences deriving from spatial sampling and pathology. Other potential influences on FC could be considered as well, including anti-epileptic medications as patients were continued on their normal regimen and were not reduced at the time of scan While simultaneously acquired icEEG-fMRI has the benefit of avoiding the confound of separate acquisitions, it is not without its own technical challenges (Carmichael et al. 2010; Carmichael et al. 2012). In particular, reduction in BOLD signal can occur around electrodes due to the difference in magnetic susceptibility between electrode and tissue in the presence of the scanner's acquisition during fMRI (Luo and Glover 2012). Given the points of consistency our BOLD results share with previous non-simultaneous icEEG-fMRI and unimodal fMRI data (Bettus et al. 2009; Bettus et al. 2011; Mankinen et al. 2011; Stufflebeam et al. 2011; Luo et al. 2014), signal drop-out due to electrodes seems unlikely to be a crucial factor in this work. Our results tend to suggest that despite the loss of MRI signal in the immediate vicinity of the electrodes, the observed effects extend Finally, while we have examined the variability of FC over time for a sliding window analysis, it should be acknowledged that establishing true dynamics requires a demonstration of a distinction with random fluctuations, and that this is particularly controversial in the context of fluctuations in BOLD-derived FC (Keilholz 2014). That said, we demonstrate relationships between haemodynamic and electrophysiological indices of variation and connectivity suggesting that at least some part of the fluctuation in BOLD connectivity reflects a neural origin. incidence of non-stationarities, we observe a lack of inter-modality connectivity correlation in regions subject to epileptic processes, with a confirmation of inter-modality discrepancies in functional connectivity associated with epileptic cortices while establishing for the first time that this discrepancy extends to the dynamical variability of connectivity over time. relationship: rapid changes in spontaneous neural activity in the somatosensory cortex are associated with haemodynamic fluctuations that resemble stimulus-evoked haemodynamics. Eur J Neurosci 38:2902-2916. doi: 10.1111/ejn.12295 Cabral J, Kringelbach ML, Deco G (2014) Exploring the network dynamics underlying brain activity during rest. Prog Neurobiol 114:102-131. doi: 10.1016/j.pneurobio.2013.12.005 Carmichael DW, Thornton JS, Rodionov R, et al (2010) Feasibility of simultaneous intracranial EEG-fMRI in humans: a safety study. NeuroImage 49:379-390. doi: 10.1016/j.neuroimage.2009.07.062 epilepsy might explain seizure semiology. Front Neurol 5:175. doi: 10.3389/fneur.2014.00175 Liao W, Zhang Z, Pan Z, et al (2010) Altered functional connectivity and small-world in mesial Luo C, An D, Yao D, Gotman J (2014) Patient-specific connectivity pattern of epileptic network in strong magnetic field573 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 574 575 576 577 578 580 581 582 583 584 585 586 590 591 592 593 594 595 626 627 628 629 630 631 632 709 710 722 to a large enough brain area beyond the region of drop out around each electrode. 579 frontal lobe epilepsy. NeuroImage Clin 4:668-675. doi: 10.1016/j.nicl.2014.04.006

Table 1 :

 1 Patient Clinical Demography. Abbreviations: yrs=years, Hand.=Handedness, M=male, F=female, R=right, L=left, NL=non-lesional, HS=hippocampal sclerosis, FCD=focal cortical dysplasia, SFG=superior frontal gyrus, MFG=middle frontal gyrus, IFG=inferior frontal gyrus, TL=temporal lobe, Sup.= superior, SMA= supplementary motor area, OrbF=orbitofrontal gyrus. ILAE surgery outcome: I: Completely seizure free, no aura; II: Only auras, no other seizures; III: One or two seizure days per years, ± auras; IV: Four seizure days per year to 50% reduction of baseline seizure days, ± auras; V: Less than 50% reduction of baseline seizure days to 100% increase of baseline seizure days, ± auras; VI: More than 100% increase of baseline seizure days, ± auras. *Hippocampal sclerosis was considered incidental in this patient. Ictal SPECT, PET, MEG and scalp EEG all support a right frontal seizure origin.

				Hand Diagnosis	Onset	Seizure	MRI	Surgical	Outcome
						(yrs)	Frequency		Resection	ILAE Months
	P01	M	26	Right	Bi-	7	Monthly	L. HS	L. TL	I	38
					Temporal
	P02	F	31	Right Temporal	6	Weekly	NL	L. TL	IV	43
	P03	M	28	Right	Frontal	12	Daily	FCD L.	L.	I	22
								SFG\	SFG/MFG
								MFG	Sup. SMA
	P04	M	38	Right	Frontal	8	Daily	L HS*	R. SMA,	I	23
									m. SFG
	P05	F	34	Right	Frontal	7	Daily	FCD L.	L.	I	19
								IFG	SFG/MFG
	P06	F	27	Left	Frontal	3	Daily	NL	L. SFG	I	24
	P07	M	24	Right Temporal	15	Weekly	NL	R.TL	IV	19
	P08	M	34	Right	Frontal	9	Daily	FCD L.	L.MFG/IF	I	11
								MFG	G
	P09	M	32	Right	Frontal	16	Weekly	NL	R.	I	22
									OrbF/IFG

Table 2 -

 2 Regions of interest/electrodes. Numbers of recording electrodes localized to each electrophysiological zone, overall (top) and for each participant (bottom). Pairwise interactions refers to functional connectivity estimates that contribute to the sample analysed at the group level, and include all possible interactions within each individual and each zone. Interactions between zones were not considered.

	Table

Swiss National Science Foundation (141165 and 140332, SPUM Epilepsy).

(Carmichael et al. 2010; Vulliemoz et al. 2011; Carmichael et al. 2012) offers the possibility of

non-stationarities via the only means possible: simultaneous acquisition.

1.1mm in diameter, and 2.3 or 2.4 mm in length. Within a given array, electrodes were spaced at

intervals of 10mm for ECoG/grids/strips and 5mm or 10mm for depth. In situ, and across all arrays

and individual implantation schemes, the mean inter-electrode Euclidean distance was 44±24mm.

Erlangen, Germany) running software version VB15, with a quadrature head transmit-receive RF coil

standard implementation of template subtraction and filtering[START_REF] Allen | A method for removing imaging artifact from continuous EEG 597 recorded during functional MRI[END_REF] in the Brain Vision

'anticorrelations'), with the advantage of not making any assumptions regarding their nature -linear 260 or otherwise(Wendling et al. 2010). Additionally, the use of h 2 allows comparison with our previous

findings in non-simultaneous recordings[START_REF] Bettus | Interictal functional connectivity of human epileptic 614 networks assessed by intracerebral EEG and BOLD signal fluctuations[END_REF]).

Ridley et al. 2015), FC increases have also been associated with worse neuropsychological outcomes

(Holmes et al. 2014). In epileptic networks, increases have been seen in the regional homogeneity of

resting BOLD oscillations(Mankinen et al. 2011), and indeed there appears to be a pattern of

dysfunctions in TLE(Oby and Janigro 2006). Furthermore, evidence from rat models suggests that

neurovascular coupling is comparable between evoked and spontaneous BOLD oscillations (Bruyns-

[START_REF] Bettus | Enhanced EEG functional connectivity in mesial 617 temporal lobe epilepsy[END_REF] Warren et al. 2010b; Bartolomei et al. 2013a). On the other hand, general linear modelling of

IEDs and other epileptic transient provided some of the earliest evidence of an interaction of epilepsy

impairment in pediatric frontal lobe epilepsy.Epilepsia 54:446-454. doi: 10.1111/epi.12044

Bruyns-Haylett M, Harris S, Boorman L, et al (2013) The resting-state neurovascular coupling

23:684-693. doi: 10.1093/cercor/bhs055

Connect 4:769-779. doi: 10.1089/brain.2014.0250

Kobayashi E, Bagshaw AP, Grova C, et al (2006) Negative BOLD responses to epileptic spikes. Hum

Brain Mapp 27:488-497. doi: 10.1002/hbm.20193

Laufs H (2012) Functional imaging of seizures and epilepsy: evolution from zones to networks. Curr

Opin. doi: 10.1097/WCO.0b013e3283515db9

Laufs H, Rodionov R, Thornton R, et al (2014) Altered FMRI connectivity dynamics in temporal lobe

temporal lobe epilepsy. PloS One 5:e8525. doi: 10.1371/journal.pone.0008525

Connectivity of the Default Mode Network and Transient Epileptic Activity. Front Neurol

5:201. doi: 10.3389/fneur.2014.00201

Lu H, Wang L, Rea WW, et al (2016) Low-but Not High-Frequency LFP Correlates with

Spontaneous BOLD Fluctuations in Rat Whisker Barrel Cortex. Cereb Cortex N Y N 1991

26:683-694. doi: 10.1093/cercor/bhu248

Lu H, Zuo Y, Gu H, et al (2007) Synchronized delta oscillations correlate with the resting-state

functional MRI signal. Proc Natl Acad Sci U S A 104:18265-18269. doi:

10.1073/pnas.0705791104

Luo Q, Glover GH (2012) Influence of dense-array EEG cap on fMRI signal. Magn Reson Med

68:807-815. doi: 10.1002/mrm.23299

Magri C, Schridde U,Murayama Y, et al (2012) The amplitude and timing of the BOLD signal

reflects the relationship between local field potential power at different frequencies. J

Neurosci Off J Soc Neurosci 32:1395 -1407 . doi: 10.1523 /JNEUROSCI.3985-11.2012

Mankinen K, Long X-Y, Paakki J-J, et al (2011) Alterations in regional homogeneity of baseline brain

activity in pediatric temporal lobe epilepsy. Brain Res 1373:221-229. doi:

10.1016/j.brainres.2010.12.004

Nedic S, Stufflebeam SM, Rondinoni C, et al (2015) Using network dynamic fMRI for detection of

epileptogenic foci. BMC Neurol 15:262. doi: 10.1186/s12883-015-0514-y

Netoff TI, Pecora LM, Schiff SJ (2004) Analytical coupling detection in the presence of noise and

nonlinearity. Phys Rev E Stat Nonlin Soft Matter Phys 69:017201. doi:

10.1103/PhysRevE.69.017201

Nir Y, Mukamel R, Dinstein I, et al (2008) Interhemispheric correlations of slow spontaneous

neuronal fluctuations revealed in human sensory cortex. Nat Neurosci 11:1100-1108.

Nissen IA, van Klink NEC, Zijlmans M, et al (2016) Brain areas with epileptic high frequency

oscillations are functionally isolated in MEG virtual electrode networks. Clin Neurophysiol

Off J Int Fed Clin Neurophysiol 127:2581-2591. doi: 10.1016/j.clinph.2016.04.013

Noachtar S, Rémi J (2009) The role of EEG in epilepsy: a critical review. Epilepsy Behav EB 15:22-

doi: 10.1016/j.yebeh.2009.02.035

Oby E, Janigro D (2006) The blood-brain barrier and epilepsy. Epilepsia 47:1761-1774. doi:

10.1111/j.1528-1167.2006.00817.x

PalminiA (2006) The concept of the epileptogenic zone: a modern look at Penfield and Jasper's

views on the role of interictal spikes. Epileptic Disord Int Epilepsy J Videotape 8 Suppl

2:S10-15.

Pan W-J, Thompson G, Magnuson M, et al (2011) Broadband local field potentials correlate with

spontaneous fluctuations in functional magnetic resonance imaging signals in the rat

somatosensory cortex under isoflurane anesthesia. Brain Connect 1:119-131. doi:

10.1089/brain.2011.0014

Epilepsia 43:219-227.

Sporns O (2014) Contributions and challenges for network models in cognitive neuroscience. Nat

Neurosci 17:652-660. doi: 10.1038/nn.3690

Stufflebeam SM, Liu H, Sepulcre J, et al (2011) Localization of focal epileptic discharges using

functional connectivity magnetic resonance imaging. J Neurosurg 114:1693-1697. doi:

10.3171/2011.1.JNS10482

Tagliazucchi E, Laufs H (2015) Multimodal imaging of dynamic functional connectivity. Front

Neurol 6:10. doi: 10.3389/fneur.2015.00010

Thompson GJ, Merritt MD, Pan W-J, et al (2013) Neural correlates of time-varying functional

connectivity in the rat. NeuroImage 83:826-836. doi: 10.1016/j.neuroimage.2013.07.036

requirements of the Declaration of Helsinki and the Joint Research Ethics Committee of the NHNN 31 (UCLH NHS Foundation Trust) and UCL Institute of Neurology. 32 33 Data for this was acquired at UCLH/UCL who received a proportion of funding 34 from the Department of Health's NIHR Biomedical Research Centres funding scheme. Processing and 35 analysis were undertaken at CRMBM/CEMEREM. This work was also supported with funds from the part in this study, which was approved by the Joint Research Ethics Committee of the National 164 Hospital for Neurology and Neurosurgery (NHNN, UCLH NHS Foundation Trust


F Radjai

D E Wolf

S Roux Lpmmh

M Jean

J J Moreau

Force networks 1n dense granular media

U sing Contact Dynamics simulations of quasistatically driven assemblies of rigid particles we show that the contact network at every stage of deformation is composed of two complementary subnetworks: a "strong" percolating subnetwork of the contacts carrying a force larger than the mean force, and a "weak" subnetwork of the . contacts carrying a force lower than the mean force. In the strong subnetwork, all contacts are nonsliding and forces have a decreasing exponential distribution. Almost the whole dissipation takes place inside the weak subnetwodc, where the number of normal forces decays as a power-law. The strong subnetwork supports all the deviatoric load, whereas the weak subnetwork contributes only to the mean pressure. The orientation of the induced anisotropy (major principal axis of the fabric tensor) in the strong subnetwork is the same as that of the stress tensor, but the orientation of the induced anisotropy in the weak subnetwork is orthogonal.

INTRODUCTION

The internal mechanical state of a random packing of noncohesive particles is doubly heterogeneous:

1. Contact forces vary by several orders of ma gnitude on a scale far larger than the particle size [START_REF] Dantu | Proceedings of the 4th Interna tional Conference on Soil M echanics and Foun dations Engineering[END_REF][START_REF] Liu | [END_REF];

2. Frictional dissipation ( due to sliding) in qua sistatic deformation occurs only at a small fraction of contacts.

Micromechanical modeling of granular materi als involves a quantitative characterization of both statistical distributions and spatio-temporal cor relations of internal variables on the scale of these heterogeneities. In this paper, on the basis of nu merical results from simulations of a few thousands of particles in two dimensions, we show that the organization of these variables has a basic bimodal char acter, which is also reflected in the structure of the stress tensor. We will first briefly intro duce the simulation method. Then, we will focus successively on force distributions, stress tensor, induced anisotropy, and dissipation.

SIMULATION METHOD

A numerical simulation of a densely-packed system has to tackle a multicontact problem. Every colli sion in such a medium is a multiple collision that cannot be reduced to a set of independent binary collisions since impulsions propagate through the contact network and may even leave the system. In the same way, the frictional resistance to shear is a collective phenomenon involving the mobili zation of friction at kinematically correlated con tacts. Another basic problem is that the contact laws are strongly nonlinear. In the approximation of perfectly rigid particles, the normal force at a contact is not given locally as a fonction of the relative displacement of two particles, but rather as a result of the global geometrical configuration of the whole system and the boundary conditions. Furthermore, the basic Coulomb's law of friction is nonsmooth in the sense that the relation between the relative tangential velocity and the friction force cannot be represented as a mathematical fonction. The Contact Dynamics (CD) method, which was used for the investigations reported in this paper, takes these features into account on a mathematical basis derived from Convex Analysis (Moreau 1994[START_REF] Jean | M echanics of Geometrical Interfaces[END_REF]). On one hand, the con ditions of perfect rigidity and exact Coulombian friction are implemented with no resort to any regularization. On the other hand, all kinematic constraints are simultaneously taken into account, together with the equations of dynamics, in order to determine contact forces and particle velocities in the system. The method is thus able to deal properly with the nonlocal character of the mo mentum transfers in the contact network.

FORCE DISTRIBUTIONS

Let us begin with the probability distribution P N of normal forces N independently of contact orien tations. Numerical results will be presented here for three samples of 4012, 4025, and 4098 partic les, referred to as samples A, B, and C. Particles are contained in a rectangular frame composed of two horizontal and two vertical rigid walls. Sam ple A is biaxially compressed by allowing the in ward motion of a horizontal wall and by applying a confining load on a vertical wall free to move horizontally. Sample B is in relaxation towards static equilibrium under a confining load applied to a free wall. Sample C is in static equilibrium under confining loads. The acceleration of gravity and the particle-wall coefficients of friction are set to zero. The interparticle coefficient of friction is 0.5 in samples A and C, and 0.2 in sample B. In all samples, the particle radii are uniformly distri buted between 3.8 and 7.5 mm.

Figs. 1( a) and l(b) display the plot of P N on semilog and log-log scales in the three samples. In all cases, independently of the confining load, the normal forces lower than the mean normal force (N) have a power-law distribution, whereas the data for forces larger than ( N) are well fitted by an exponential decay:

{ ( N ) -(){ P N oc ( N ) e (3(1-N/ ( N ) ) N < (N), N > (N).
(1)

Within the statistical precision, the values of exponents a and (3 depend only weakly on the me chanical parameters and the preparation conditi ons of each sample. In static equilibrium or in qua sistatic deformation, the value of a is very close to 2 zero. However, it seems that a increases with the degree of dynamics inside the system and decreases with the coordination number. Let us also remark that, since the mean value (N) separates the two parts of the distribution, a and (3 should be related � ogether by the equation /3 2 = (l-a)(2-a), which 1s approximately satisfied in our simulations.

As far as the distribution P r of the absolute values of friction forces T is concerned, we always find a power law decay f or forces lower than the mean friction force (T) and a decreasing exponen tial fonction for friction forces higher than (T), as shown in Figs. 2(a) and 2(b). This distinction between the forces lower than the mean, to which we will refer as "weak" forces, and those larger than the mean, referred to as "strong" forces in the following, has a deep meaning, as we shall see below. Let us only mention here that the expo nential distribution of strong forces does not seem to be a finite-size effect. In other words, a sim ple normalization of the forces with respect to the mean force allows the data to collapse almost on the same curve (Radjai et al. 1996).

The exponential distribution of strong normal forces, has been observed in experiments [START_REF] Liu | [END_REF], Miller et al. 1996). Weak normal forces and friction forces at individual contacts are tech nically difficult to measure, and their distributions have not yet been observed experimentally. The power-law decay or the uniform distribution of weak forces, comprising nearly 60% of contacts, in dicates that the weak forces are generated through a self-similar branching process from strong forces. In contrast, the strong forces are conditionned by the level of the deviatoric load, which sets in ef fect a characteristic force different from the mean force. This point will become more clear by the study of the stress tensor.

STRESS TENSOR

The stress tensor is the relevant variable in tran sition to the continuous limit and, in contrast to PN and Pr, involves both the contact forces and the contact orientations according to the follo wing classical formula (Christoffersen et al. 1981):

O"ij = V(rl) =cErl fr dj, where 0, is the region in sicl e the sample for which the stress tensor is cal culated with a volume equal to V, f { is the i component of the force at the contact c and dC. is the j-component of the intercenter v;ctor d/_ This equation allows to separate the contributions of weak and strong forces to the global stress ten sor. The weak tensor CT w and the strong tensor CT 8 are given by the above equation when the summa tion is restricted to only weak contacts and only strong contacts, respectively.

Among our samples, only sample A, due to its uniform stress field, provides enough statistics for the study of the stress tensor. Figure 3 shows the orientations and the eigenvalues of CT w and CT 8 in sample A at shear peak, i.e. at the maximal va lue of q = �! ��;, ( q = 0.18), where CT 1 and CT 2 are the principal values of the total stress tensor. We see that CT w is almost spherical and contributes no more than 3% to the deviatoric part of the total stress tensor CT = CT w + CT 5 , whereas CT 8 represents al i:n o � t the :v hole deviatoric load, and its major prmc1pal axis has the same orientation as that of the imposed deformation.

The isotropie nature of the weak stress tensor indicates that the weak contacts feel the influence o � the . deviatoric load only as an average over all directwns. The strong contacts directly support the deviatoric load, whose signature appears as a characterstic force in the exponential decay of strong forces.

INDUCED ANISOTROPY

We now consider the distribution P o of contact ori entations. The anisotropie nature of P o has been e � tensively studied in the past (Biarez & Wien d1ck 1963[START_REF] Satake | Theoretical and Applied Me chanics: 257[END_REF], Cambou & Sidoroff 1985, Rothenburg & Bathurst 1989). I-Iere, we define two separate distributions P 8 and Pô for weak and strong contacts, respectively. The polar dia grams of these distributions are shown in Figure 4 for sample A at shear peak. We see that both distributions are anisotropie, although to a lesser extent in the weak subnetwork than in the strong one. The interesting behavior observed here is that the principal direction of Pô coïncides with that of the stress tensor, whereas the principal direc tion of P 8 is orthogonal! This again confirms the bimodal nature of the force network, and clearly shows the complementarity between the two sub networks: The strong chains need "lateral ,, weak forces in order to be stabe. In a quasistatically driven system, the particle ro tations are generally frustrated, with patterns that have been partia. lly studied in the case of regular packings (Radjai & Roux 1995, Radjai et al. 1996).

The possibility of particles to roll over one another provides a very low-dissipative local mechanism of deformation, so that if each particle could roll over its neighbors, then a dense granular system would essentially beha. ve like a liquid. The question we raise here is how the positions of sliding contacts on the contact network are correlated with the con tact forces?

The distinction between weak and strong sub networks again provides a simple key to the pro blem. In Figure 5, we have shown the two subnet works with different gray levels and the positions of the sliding contacts. Almost 8% of contacts are sliding, and a rapid inspection shows that al most all of them are on the weak subnetwork. In other words, all strong contacts are nonsliding and all the dissipation takes within the weak sub network. This remarka. bly simple behavior of the sliding contacts with respect to the force subnet works shows the dynamic nature of the bimodal d;;:,tribution of forces. The strong subnetwork is completely unstable since all i ts contacts are nons-, liding. As a result, it tends to buckle un der the action of the external deviatoric load, giving rise to an orthogonal anisotropy and sliding contacts in the weak subnetwork.

CONCLUSION

The observations briefly presented above provide a simple picture of the interna! state of a granular packing in quasistatic deformation, at least in two dimensions. The contact network is composed of two complementary subnetworks: 1. A "load-bearing" percolating subnetwork composed of strong contacts; 2. A "dissipative" subnetwork composed of weak contacts.

The strong subnetwork supports the whole de viatoric load applied on the system, although it con tains only nonsliding ( nondissipative) contacts. The whole dissipation takes place inside the weak subnetwork. This implies that the scale of hete rogeneity of a gra.nular system due to sliding and nonsliding states of contacts is the same as that of the strong subnetwork. From the point of view of the stress tensor, the weak subnetwork behaves like an intersticial liquid, whereas the strong sub network has a solid-like behavior since it takes over the whole deviatoric load and thus the stability of the system. From the point of view of the fabric tensor ( or contact orientations), the strong sub network bears the primary anisotropy induced by shear, but it gives rise to a secondary orthogonal anisotropy inside the weak subnetwork. This bimodal behavior of the force network opens new perspectives bath for microscopie and macros copie modeling of granular materials in quasistatic deformation. Let us, however, recall that the ana lysis presented here concerns a two-dimensional system. It is not obvious that it applies in the same manner to three-dimensional systems. The refore, the investigation of three-dimensional sy stems, even though their simulation is more invol ved, is crucial to the further establishment of the results presented in this paper.

  Figure 1. ( a) Semilogarithmic, and (b ) log-log plots of the probability distribution P N of normal forces -N.

Figure 3 .

 3 Figure 2. (a) Sernilogarithmic, and (b log-log plots of the probability distribution PT of friction forces T.

Figure 4 .

 4 Figure 4. Polar diagram of the distribution of contact ori entations in weak and strong subnetworks.

Figure 5 .

 5 Figure 5. Network of normal forces in sample A. Forces are encoded as the widths of intercenter segments. Weak and strong subnetworks are shown with two different gray levels. The filled circles show the positions of sliding contacts.