A Ferrari

C Richard

Non-parametric Community Change-points Detection in Streaming Graph Signals

Keywords: Graph signal processing, streaming graph signals, non-parametric change-point detection, graph filtering

Detecting changes in network-structured time series data is of utmost importance in critical applications as diverse as detecting denial of service attacks against online service providers or monitoring energy and water supplies. The aim of this paper is to address this challenge when anomalies activate unknown groups of nodes in a network. We devise an online change-point detection algorithm that fully benefits from the recent advances in graph signal processing to exploit the characteristics of the data that lie on irregular supports. Built upon the kernel machinery, it performs density ratio estimation in an online way. The algorithm is scalable in the sense that it is spatially distributed over the nodes to monitor large-scale dynamic networks. The detection and localization performances of the algorithm are illustrated with simulated data.

I. INTRODUCTION

The problem of change-point detection is defined as finding a switching time instant from which the statistical properties of a signal change, namely, the observations are drawn from a distribution up to the point before it, and are drawn from another distribution after this point. Detecting anomalous events in network-structured time series data has become relevant in a variety of applications. Examples can be found in network security such as detection of denial of service attacks or cyberintrusions in computer networks. Examples can also be found in medical monitoring, such as automatic detection of epileptic seizures in electroencephalogram signals, or detection of mental workload or drowsiness in brain computer interfaces. This problem is of great interest in a wide range of disciplines, such as physics, biology, and social sciences, where graphs provide a powerful machinery for effectively capturing correlations among inter-dependent time series data.

The anomalies in network-structured time series data often tend to activate communities of nodes in the graphs. The problem then consists of deciding, based on noisy measurements at each node of a graph, whether the underlying unknown graph signal is in a nominal state over the whole graph, or if there exist one or more clusters of nodes with anomalous activation. Over the last years, graph signal processing (GSP) has emerged as a new research field aiming at extending the well-developed tools for processing time-domain signals to signals defined over the nodes of a graph while exploiting the underlying connectivity information. Beyond its use for sampling [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF], filtering [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] or modeling [START_REF] Nassif | Distributed diffusion adaptation over graph signals[END_REF] for change-points detection over graph signals [START_REF] Sharpnack | Changepoint detection over graphs with the spectral scan statistic[END_REF], [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF]. These works address critical issues raised in [START_REF] Addario-Berry | On combinatorial testing problems[END_REF], [START_REF] Arias-Castro | Detection of an anomalous cluster in a network[END_REF]. Indeed, if the generalized likelihood ratio (GLR) test is a natural solution for detecting anomalous clusters of activity in graphs, it has high computational complexity because it requires scanning all well-connected clusters of nodes and testing them individually. In [START_REF] Sharpnack | Changepoint detection over graphs with the spectral scan statistic[END_REF], the authors first address this issue by incorporating, into the detection problem, the properties of the graph topology through its spectrum. They analyze the test statistics, which is based on spectral patterns of the combinatorial Laplacian, and show how it is indeed related to the problem of finding graph cuts. Following the same line of reasoning, they derive another detector in [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF], called Graph Fourier Scan Statistic (GFSS). Based on a relaxation of its initial form in [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF], it consists of a low-pass filter built upon the graph Fourier transform. Because of its particular form, the performance of this detector can be precisely characterized. None of these works address the online non-parametric detection problem from streaming graph signals, nor do they implement the detectors in a distributed manner.

The aim of this paper is to propose a flexible strategy built upon the (reproducing) kernel machinery to perform density ratio estimation in an online and distributed manner across nodes, while using GSP framework to enhance the contrast produced by possible anomalous activities over unknown groups of nodes and then make their detection easier.

II. PROBLEM FORMULATION

We consider an undirected graph G = {V, E, W } with p vertices V = {1, . . . , p}, m edges (i, j) ∈ E ⊂ V × V , endowed with a p × p weighted adjacency matrix W where W i,j = W j,i ≥ 0 is the connection strength of (i, j) ∈ E.

We denote by y t (i) the p i -dimensional signal sample observed at time instant t and vertex i. We shall assume that the observations y t (i) are independent w.r.t. t and i. In this paper, we consider that a change in the distribution of y t = y t (1) . . . y t (p) may occur in an unknown cluster C * at an unknown time instant t r , that is,

t < t r : y t ∼ p i=1 f i y t (i) , t ≥ t r : y t ∼ p i=1 f + i y t (i) ,
where ∀i ∈ C * :

f + i = f * i , and ∀i ∈ C * : f + i = f i . (1)
The objective is to derive an online algorithm to process the streaming data y t in order to detect a change-point as defined in [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF], where C * , t r and all f i and f * i are unknown.

A solution to this problem was recently proposed in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF] for scalar measurements (p i = 1) and model y t (i) = m t (i) + e t (i) where e t (i) is a scalar i.i.d zero-mean Gaussian noise. The changes under consideration occur on the mean vector m t as:

t < t r : m t = m, t ≥ t r : m t = m + δ, (2)
where δ is a non-zero vector only on C * which does not vary within this cluster. Built upon GSP to extract information about the clusters, the core detection functionality implemented at each node in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF] is ensured by the adaptive strategy proposed in [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF]. Two linear models are learned simultaneously using two learning rates, on the one hand to estimate the baseline trend, and on the other hand to adapt to abrupt changes. The resulting test statistics, which consists of the distance between these two models, only applies to specific scenarios such as ( 2).

Here we address problem (1) with a flexible non-parametric strategy in the case where a) the measurements at each node are multivariate and b) all distributions f i and f * i are unknown.

III. GFSS APPROACH TO CLUSTERED DETECTION

We first recall the principles of GFSS detection introduced in [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF], with a focus on the filtering process to perform spectral clustering. Next, we introduce the likelihood ratio used to address problem [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF]. None of these ingredients address the online detection problem from streaming graph signals, nor do they implement the detectors in a distributed manner.

A. GFSS spectral clustering

Work [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF] focuses on a static graph signal consisting of a scalar measurement at each node as in [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] with m = 0. Let L be the graph Laplacian of G and let us define by {u i } i=1,...,p the set of orthonormal eigenvectors of L with {µ i } i=1,...,p the corresponding eigenvalues. Given y, the GFSS is defined as:

t GFSS (y) = g y 2 , (3)
where g y is the graph-filtered signal:

g y = p i=2 h * (µ i )(u i y)u i , (4)
with h * (µ) the frequency response of the filter defined in [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF]:

h * (µ) = min 1, γ µ , µ > 0, (5)
and γ > 0 a tuning parameter; see [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF], [START_REF] Tremblay | Design of graph filters and filterbanks[END_REF], [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF] for details on graph filtering.

To get more insight into statistics (3), it is important to recall the role played by the eigenvectors u i of the graph Laplacian matrix L in spectral clustering. Consider first the ideal case of a graph with k > 1 connected components or communities, and denote by C i with i = 1, . . . , k the set of vertices in component i. Each eigenvector u i of the Laplacian matrix is proportional (with the factor 1/ |C i |) to the indicator function of the ith connected component C i , and u i y is proportional to the sum of y on C i . Consequently, (u i y)u i in (4) assigns to each vertex n of C i , the mean value of the signal samples y(n) at all vertices n within C i . Finally, as the number of components C i is unknown, the frequency response in (5) penalizes large numbers of components in [START_REF] Sharpnack | Changepoint detection over graphs with the spectral scan statistic[END_REF]. When the k components C i are connected by few edges, we can assume that this analysis which relies on the properties of the eigenvectors u i of the Laplacian matrix is still approximately valid. Note that this kind of approximation is the starting point of spectral clustering methods; see, e.g., [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], [START_REF] Shi | Normalized cuts and image segmentation[END_REF].

B. Likelihood ratio test

Consider problem [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF], and assume that all probability densities are known. We denote by t the p-dimensional vector containing the log-likelihood ratios of the measurement y t (k):

t (k) = log f + k (y t (k)) f k (y t (k)) . (6)
The likelihood ratio test is a natural solution for addressing a change-point detection problem such as [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF]. Nevertheless it suffers from its computational complexity since it consists of scanning all possible well-connected clusters of nodes C k and testing them individually. Substituting y by t in (4), and considering the resulting graph signal g t , it can be expected that the entries of g t will be approximately equal to:

g t (k) ≈ 1 |C i | log n∈Ci f + n (y t (n)) n∈Ci f n (y t (n)) , (7)
for all node k in C i . This means that, by exploiting the graph topology, g t assigns to each node of G the log-likelihood ratio of the cluster to which it belongs. The clustered structure of the graph is promoted in the log-likelihood ratio g t by the graph filter [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF], which makes its use attractive to address the change-point detection problem [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF].

IV. DETECTION ALGORITHM

The approach outlined in Section III relies on the knowledge of the likelihood ratio at each node n and time instant t. We shall now devise a non-parametric detection scheme based on an estimate ˆ t of t . Computing g t from t using (4) needs to be carried out in a centralized manner. In order to make the detection algorithm scalable, we shall investigate a distributed implementation across the nodes of the network.

A. Online likelihood ratio estimation

Non-parametric change-point detection algorithms were introduced to handle scenarios where no prior information on the data distribution and the nature of the change is available. They are usually classified into two categories, supervised and unsupervised methods, depending if training data are available. This paper focuses on unsupervised algorithms as their flexibility allows them to handle many practical situations. Some of them are based on subspace or manifold tracking methods [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF], [START_REF] Chen | Sequential change-point detection based on nearest neighbors[END_REF], [START_REF] Xie | Change-point detection for highdimensional time series with missing data[END_REF], and aim to detect abrupt changes in the sequence of residuals. Their parametric form makes them pre-designed for tracking specific statistics. Others are model-free as they do not rely on particular parametric assumption [START_REF] Kifer | Detecting change in data streams[END_REF]. They however tend to be less accurate in high-dimensional problems because of the curse of dimensionality. To overcome this difficulty, recent contributions have focused on estimating the ratio of probability densities directly without going through density estimation [START_REF] Sugiyama | Density Ratio Estimation in Machine Learning[END_REF], [START_REF] Aminikhanghahi | Real-time change point detection with application to smart home time series data[END_REF], [START_REF] Liu | Change-point detection in time-series data by relative density-ratio estimation[END_REF]. In this paper, we are interested in an algorithm, denoted as NOUGAT (Nonparametric Online chanGepoint detection AlgoriThm) which implements this strategy [START_REF] Bouchikhi | Nonparamametric online change-point detection with kernel lms by relative density ratio estimation[END_REF], [START_REF]Kernel based online change point detection[END_REF]. It estimates the density ratio on two consecutive intervals of the time series: a set of recent samples, say, the test set, and the set of samples that came before, say, the reference set. Unlike the algorithm considered in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF], it naturally manipulates multivariate observations at each node.

NOUGAT aims to detect changes in the distribution of y t (k) at each vertex k by estimating a model for r t (k) -1, where

r t (k) = f + k (y t (k))/f k (y t (k))
is the density ratio (with some abuse of notation) between the probability densities f + k (y) and f k (y) of the data on a test and reference interval, respectively:

y t-(Ntest-1) (k), . . . , y t (k) ∈ R p k ×Ntest , ( 8
)
y t-(Nref+Ntest-1) (k), . . . , y t-Ntest (k) ∈ R p k ×Nref , (9)
with N test and N ref the number of samples in the test and reference intervals, respectively. It is shown in [START_REF]Kernel based online change point detection[END_REF] how r t (k) -1 can be estimated, in a reproducing kernel Hilbert space K endowed with kernel κ(• , •), by a kernel model. Note that r t (k) -1 is preferred to r t (k) because it leads to an unbiased estimator under the no change-point hypothesis.

Using this model, we arrive at the estimation of t (k) in ( 6):

ˆ t (k) = log θ t-1 (k) κ k ω (y t (k)) + 1 , (10)
κ k ω (•) = κ k ω1 (•), . . . , κ k ω L (•) , (11)
where κ k ωq (•) = κ(• , y ωq (k)), with q ∈ {1, . . . , L}, are the elements of a dictionary of size L which ensures the executability of the algorithm. This dictionary can be designed in an online or offline way using a standard strategy as in [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. For the sake of clarity, we shall assume in the sequel that the dictionary has been designed beforehand.

Similarly to the derivation of the KLMS in [START_REF] Richard | Online prediction of time series data with kernels[END_REF], the parameter vector θ t (k) in (10) can be estimated recursively at each vertex k by a stochastic gradient descent algorithm:

θ t (k) = θ t-1 (k) -µ (H t (k) + νI)θ t-1 (k) + e t (k) , ( 12
)
where µ is the adaptation step and ν is a ridge regularization parameter. Quantities H t (k) and e t (k) are defined as:

e t (k) = h ref t (k) -h test t (k), (13)
h test t (k) = 1 N test t i=t-(Ntest-1) κ k ω (y i (k)), (14)
h ref t (k) = 1 N ref t-Ntest i=t-(Ntest+Nref-1) κ k ω (y i (k)), (15)
H t (k) = 1 N ref t-Ntest i=t-(Ntest+Nref-1) κ k ω (y i (k))κ k ω (y i (k)). (16)

B. Distributed implementation

The test statistics t defined in [START_REF] Addario-Berry | On combinatorial testing problems[END_REF], and estimated in [START_REF] Tremblay | Design of graph filters and filterbanks[END_REF], does not take into account the graph topology. The aim of Section III was to show how graph-filtering t to g t improves the localisation of the cluster containing the abrupt change. The filtering operation defined in (4) requires centralized computational resources and hence cannot be scaled to large networks. Distributed implementations of graph filters have received much attention in the literature [START_REF] Shuman | Chebyshev polynomial approximation for distributed signal processing[END_REF], [START_REF] Sandryhaila | Finite-time distributed consensus through graph filters[END_REF], [START_REF] Segarra | Distributed implementation of linear network operators using graph filters[END_REF], [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF], [START_REF] Coutino Minguez | Advances in distributed graph filtering[END_REF]. They are based on the approximation of the filtering operator by a polynomial in L for which, similarly to a FIR filter, the computation of the n-th entry of L t depends on t (n) and on t (m) where m ∼ n. Recently, [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF], [START_REF] Coutino Minguez | Advances in distributed graph filtering[END_REF] proposed novel filter approximation strategies where, as in an IIR filter, each node recursively aggregates the previous outputs of its neighbors, allowing for the computation of a larger family of frequency responses with a limited memory footprint.

These filters iterate on a static graph signal. In the context of online change point detection, we propose to apply the ARMA K filter defined in [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF] to the streaming signal ˆ t , t ≥ 0, defined in (10):

x q,t = ψ q Lx q,t-1 + ϕ ˆ t , x q,-1 = 0, ∀q = 1 . . . K, (17)

ĝˆ t = K q=1 x q,t + c ˆ t . (18)
If we assume that: i). there is no change point, i.e. t = , ii) ˆ t is an unbiased estimate of t , see [START_REF]Kernel based online change point detection[END_REF] and iii) max {|ψ |}ρ(L) < 1 with ρ(•) the spectral radius of its matrix argument, it can be shown that lim t→∞ E{ĝ ˆ t } = ĝ where ĝ is graph filtered by

h(µ) = c + K q=1 ϕ q 1 -ψ q µ . ( 19
)
The proof is omitted in this paper due to lack of space. This suggests to estimate g t by ĝˆ t where c and {(φ q , ψ q )} q=1,...,L are compute in order that h(µ) approximates h * (µ) in [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier Scan Statistic[END_REF]. The corresponding distributed change point detection algorithm is described in Alg. 1. Update θ t (n) using ( 12)-( 16)

V. EXPERIMENTS

The graph G used for simulations is provided in Fig. 1a. It is unweighted and contains p = 250 vertices and m = 2508 edges. In order to generate y t as in (1) with a change on C * , the communities C i of G were unfolded using [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. Cluster C * is colored in green khaki color at the center right of the graph.

The change-point (1) was set to t r = 1400 and defined as:

∀i ∈ C * : f + i ≡ N (0, C + i ), and ∀i ∈ C * : f i ≡ N (0, I),
with:

C + i = 1 0.9 0.9 1 . (20)
The parameters of NOUGAT were set to N ref = N test = 128, µ = 0.01 and λ = 0.01. Note that given these windows lengths, see ( 8)-( 9), NOUGAT detection is expected at t * = t r +128. A Gaussian kernel κ with bandwidth √ 0.3 was used. The kernel dictionary size was set to L = 80 with 40 elements sampled from f i and 40 elements sampled from f + i . Parameter γ of the GFSS filter (5) was set to γ = 0.3, and K in ARMA K approximation to K = 4 [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF].

Figure 2 evaluates the detection performance of Algorithm 1 by comparing ˆ t clearly benefits from the graph topology through GFSS filter. In order to show that this improvement comes from the ability of ĝˆ t to better localize the nodes where a change occurs, Figs. 1b and1c represents the two test statistics on G at time t = t * . To confirm this result Monte Carlo simulation were performed to estimate ROC curves in both cases. Probabilities of detection and false alarm are defined (e.g. for | ˆ t * (n)|) as:

P d = Prob | ˆ t * (k)| > ξ | k ∈ C * , (21)
P fa = Prob | ˆ t * (k)| > ξ | k ∈ C * . ( 22
)
Figure 3 clearly shows the benefit of using graph filtering in the proposed algorithm.

VI. CONCLUSION

This paper introduces a novel algorithm to detect a change point on a streaming graph signal which is located on an unknown cluster of the graph. The algorithm is fully nonparametric: it does not require any model of the change point and applies to multivariate measurements on the nodes of the graph. It also achieves scalability to large-scale dynamic networks through: i) a fully distributed implementation across the nodes, ii) an online processing of the measurements. The detection and localization performances of the algorithm are illustrated using simulated data.

  signals on graphs, to cite a few key research directions, GSP was recently used This work was funded in part by ANR under grant ANR-19-CE48-0002.

Algorithm 1 :

 1 Distributed graph change point detection 1: for t = 1, 2, . . . do 2:// Each node updates separately its likelihood ratio 3:for n = 1, . . . , k do 4:

Fig. 1 :

 1 Fig. 1: Graph topology with coloured clusters C i and test statistics over the graph.

2 2 and ĝˆ t 2 2

 22 , for illustration purposes only since these two statistics are global. As expected, ĝˆ t 2 2

g t 2 2 Fig. 2 :Fig. 3 :

 223 Fig. 2: Global change point detection statistics over G.


Muhammad Arshad

Abdelhamid Maali

Cyrille Claudet

Laurent Lobry

Francois Peters

Elisabeth Lemaire

An experimental study on the role of inter-particle friction in the shear-thinning behavior of non-Brownian suspensions †

µm in diameter), immersed in a Newtonian liquid, for different normal loads ranging from 10 to 1000 nN. It is shown that the inter-particle friction coefficient decreases with the load, contrarily to what is expected for macroscopic contacting bodies. The experimental friction law is then introduced into the viscosity model proposed by Lobry et al. and the results are compared to the viscosity of suspensions made of the same particles dispersed in the same liquid as those used for AFM measurements. The very good agreement between the measured viscosity values and those predicted by the model of Lobry et al. with the friction coefficient measured by AFM as input data shows the relevance of the scenario proposed by Lobry et al. and highlights the close links between the microscopic friction properties of the particles and the macroscopic rheological behavior of suspensions.

Introduction

The understanding and physical description of the rheology of concentrated suspensions has undergone important developments over the last 10 years, in particular by taking into account the frictional contacts between particles 1,2 . Several numerical [3][4][5][6][7] , theoretical 8 and experimental [9][10][11] works have evidenced the marked influence of the inter-particle friction on the rheology of suspensions. In particular, friction is manifested by a strong increase of the suspension viscosity.

These studies resulted in a considerable progress in describing and interpreting the non-Newtonian behaviors manifested by concentrated non-Brownian suspensions. In particular, the scenario proposed by Seto et al. 12 , Mari et al. 4 and Wyart & Cates 8 to explain the discontinuous shear-thickening (DST) observed in some highly concentrated non-Brownian suspensions turned out to be very valuable. This scenario is based on the existence of a transition from a flow regime where the particles do not come into direct contact (frictionless rheology) because the shear stress is not sufficient to allow the particles to counteract the repulsive inter- actions to a flow regime where the particles experience frictional contacts. This model has given rise to several experimental validations, including those of Guy et al. 13 , Clavaud et al. 14 , Comtet et al. 9 and Hsu et al. 15 . In the first two cases, the authors compare the rheological behavior of several suspensions whose DST characteristic stress is tuned either by varying the particle size 13 or the range of the repulsive force between particles 14 . Comtet et al. 9 present a direct measurement of the inter-particle friction by atomic force microscopy which reveals two regimes: for small loads, the contact is lubricated whereas, beyond a critical normal force, the contact appears to become frictional. The authors establish a quantitative relationship between this critical load and the stress associated with DST. At last, Hsu et al. 15 combines the tuning of the DST transition by varying the roughness of the particles and measurements of inter-particle friction coefficient by AFM, using an accurate method they had previously developed 16 . Their results illustrate in an extremely clear way the strong links between the rheology of suspensions and tribological properties of the particles. Another recent paper 17 that combines AFM and rheometric measurements shows how adhesion modifies the DST transition.

The shear-thinning behavior exhibited by most non-Brownian suspensions has been less studied but it has been often reported in the literature [18][19][20][21][22][23] . In colloidal suspensions, when adhesive forces are present between particles, shear-thinning can be explained by a competition between the adhesive forces under the effect of which the particles tend to form flocs and the hydrodynamic forces which break these aggregates. These attractive colloidal forces are known to be active in Brownian suspensions but their effects in non-Brownian suspensions have been little studied. It is only recently that Richards et al. 23 used the constraintbased rheology model previously proposed by Guy et al. 13 to analyze yielding in adhesive and frictional non-Brownian suspensions. The authors stress the differences between this class of materials where adhesion acts as a rolling constraint between contacting particles (sticky contacts) and attractive Brownian suspensions where particles are subject to attractive interaction without contacting. Papadopoulou et al. 22 showed, by using non-Brownian silica particles with different morphologies, surface functional groups and suspending media that shear-thinning can originate either from adhesion forces between particles or from load-dependent interparticle friction. This last scenario was first evidenced by Chatté et al. 10 who proved that the shear-thinning which takes place for higher stresses than the DST characteristic stress was correlated to the decrease of the friction coefficient of the particles as the contact force between particles increases. The authors used a quartz-tuning fork based atomic force microscope to measure the friction coefficient of two types of PVC particles immersed in a plasticizer, as a function of normal load. They also measured the viscosity of suspensions made of both types of particles and showed that there was a very strong correlation between the value of the jamming fraction (extrapolated from the viscosity measurement at a given volume fraction) and the value of the friction coefficient measured by the AFM. This was the first experimental proof that shear-thinning and friction could be related. Around the same time, Tanner et al. 24 proposed a model based on a bootstrap mechanism of friction where the friction coefficient is a decreasing function of the sliding speed between particles and Lobry et al. 25 showed that it was possible to quantitatively capture the shear-thinning by properly modeling the contact between particles. This model is based on the assumption that, in a shear flow, the particles come into contact via one or a few asperities. This idea that the contact between particles is promoted by surface asperities is now widely accepted. It has been particularly confirmed by experimental observations which evidenced the signature of such contacting asperities on particle pair trajectories 26 and pair distribution functions 27 in sheared suspensions. As a consequence, in the frame of a contact through one or a few asperities, it is expected that the friction coefficient will no longer be constant, in contrast to what happens for macroscopic contact involving many asperities 28 . In this model, roughnesses are represented by hemispheres whose radius is of the order of one thousandth of the radius of the particles. This roughness height is standard for polymer particles of a few tens of microns 29,30 . More and Ardekani 31 used the model proposed by Lobry et al. to study the influence of roughness size on the viscosity of a suspension. These two numerical studies use the single-contact friction model proposed by Brizmer et al. 32 where the contact between a plane and a perfectly smooth sphere (mimicking an asperity) is considered. Two regimes are identified: for small loads F n , the contact is assumed to be elastic and can then be modeled by the Hertz contact theory where the contact area varies as F 2/3 n . The friction coefficient which is given by the ratio between the tangential force that is proportional to the contact area and the normal force thus varies as the normal force to the power -1/3. For higher loads, the contact enters the plastic regime and the friction coefficient levels off at a constant value. The transition from the elastic regime to the plastic regime occurs at a given value of the normal load, L c , which depends on the radius of the asperities and on the material characteristics of the particles (Young modulus, yield strength and Poisson's ratio). As a consequence, this microscopic friction model introduces of a stress scale into the rheology of suspensions and thus opens up the possibility that the viscosity depends on the applied stress.

Finally, since the contact force experienced by the particles is related to the shear stress σ and the particle radius a through F n ∼ 6πa 2 σ , the elastic-plastic contact law results in two marked rheological behaviors, from shear-thinning at low stress, where the mean friction coefficient decreases, to Newtonian at high stress , where the friction coefficient levels off. Lobry et al. also proposed a correlation law, based on their simulation data, allowing to compute the whole viscosity curve from the microscopic friction law alone, making it possible to account for any particular friction law. In this model, the relevant interparticle force intensity that controls the mean interparticle friction coefficient is simply F n = 6πa 2 σ /1.69 where the scaling factor 1.69 is introduced for the model to fit with the earlier constant friction coefficient simulations of Gallier et al. 3 . The aim of this paper is to present a quantitative experimental validation of the model proposed by Lobry et al. 25 , hereafter called STIVF model (for Shear-Thinning Induced by Variable Friction model), that links viscosity to friction and, in particular, shear-thinning to load-dependent friction coefficient. For this purpose we have measured by using an AFM the pairwise friction coefficient of polystyrene particles, immersed in a Newtonian liquid, for different normal loads. The friction coefficient obtained is then introduced into the STIVF model. The viscosity of suspensions made of the same particles dispersed in the same liquid as the ones used for AFM measurements is then measured for several particle volume fractions. The very good agreement between the measured viscosity values and those predicted by the STIVF model with the friction coefficient measured by AFM as input data shows unambiguously the relevance of the scenario proposed by Lobry et al..

In section 2, we present the rheological measurements. The AFM measurement are described in Section 3 while section 4 is devoted to the comparison of the experimental results with the viscosity model that includes a load-dependent friction coefficient between particles.

Viscosimetric measurements

Materials

The suspensions are made of polystyrene particles (Dynoseeds TS40 from Microbeads) suspended in a mixture of water(90 vol%), Ucon oil (75-H-90,000, Dow) and a 1.5% aqueous solution of Zinc Bromide. Ucon oil is added to water to increase the viscosity of the suspending fluid in order to be able to perform rheology measurements over a wide range of shear stresses and still have a low Reynolds number. Zinc Bromide is added to precisely match the density of the liquid to that of the particles, which has been measured to be 1045 ± 4 kg/m 3 . The viscosity has been measured for various stresses comprised between 0.1 and 100 Pa. In this range, the suspending liquid is Newtonian with a viscosity equal to 0.055 ± 0.002 Pa.s at T = 25 o C. The nominal diameter of the particles is 40 µm. They are supposed to be monodisperse but Dai et al. 33 determined precisely the size distribution of the same kind of particles and showed a slight polydispersity with a standard deviation in size of about 4.3 µm around a mean diameter of 2a = 40.3 µm. By AFM imaging, we characterized the roughness of the particles that has been found to be around 36 nm (RMS value) 25 .

We studied the variation of viscosity with shear stress for five particle volume fractions, φ : 0.405, 0.45, 0.47, 0.49 and 0.507 and we compare the behavior of suspensions whose particles have been previously washed (WP) or not (RP). WP particles are washed by several steps (≈ 5) of filtration and redispersion in deionized water after which they are dried in a vacuum oven during 12 hrs at T = 60 o C. Dynoseeds TS 40 (MICROBEADS) particles are manufactured by using a technique of aqueous dispersion polymerisation that involves the presence of a stabilizer in the dispersion [START_REF] Nyhus | Process for the production of particular polymers[END_REF] . Thus, by washing the particles, we anticipate that the colloidal interactions between particles will be modified. The role of colloidal forces in suspension rheology has been much less studied for non-Brownian suspensions than for Brownian suspensions, probably because it is much weaker. However, recent studies have shown that changing the physico-chemical properties of non-Brownian suspensions has a noticeable impact on rheology 14,22,23 .

Rheometric experiments

The viscosity has been measured in rotating parallel plate geometry. The radius of the disks, R, is 30 mm and the gap height is set to 1.2 mm, i.e. 30 particle diameters, so that wall ordering effects can be neglected 35 . Rotating plate geometry provides the advantage that no migration 36 or very slow migration 37 takes place. The inconvenience of using such a geometry to study non-Newtonian materials is obviously the variation of the shear rate with radial position. To account for shear rate variation from 0 at the center to γR = ΩR/h at the rim (Ω being the angular velocity of the rotating plate), the Mooney-Rabinovitch correction is used to correct the shear stress:

σ R = σ rheo 4 3 + d ln σ rheo d ln γR (1)
where σ rheo = 2Γ/πR 3 with Γ the torque applied by the rheometer.

The viscosity is then obtained :

η = σ R γR (2)
Viscosity is measured by imposing stress steps whose duration varies depending on the intensity of the stress. Indeed, the steady viscosity of a suspension is only reached when its microstructure is at equilibrium, which requires that the suspension has been sheared over a deformation of several units. In the following, we will use the normalized viscosity defined by the ratio of the suspension viscosity to the suspending fluid viscosity: η S = η/η 0 . The variation of η S with particle volume fraction is shown in Fig. 1. In this figure, two different behaviors can be observed, depending on whether the particles have been washed (WP) before the suspension was made, or if raw particles (RP) where used.

As said before, washed particle suspensions are likely to be less stable than non-washed suspensions since washing the particles removes all or part of the surfactant. The particles are then possibly subject to adhesive forces. An order of magnitude of these forces can be obtained by using the JKR theory 38 :

F adh = 3 2 πwr c ( 3
)
where w is the energy per contact area. Following Hodges et al. 39 , for r c , we take the curvature radius of an asperity which is of the order of a hundred nanometers 25,40,41 . The value of w is found in the literature 42,43 and is of the order of 60 to 80 mJ/m 2 . These values lead to F adh ≈ 30 nN which corresponds to a typical shear stress of σ adh = F/6πa 2 ≈ 4 Pa. This order of magnitude is consistent with what is observed in Fig. 1 where the viscosity values of the two types of suspensions are shown to be different for stresses typically lower than 5 Pa. In this low stress range (σ < 5Pa), the viscosity of the WP suspensions increases sharply as the shear stress decreases, regardless of the volume fraction value. On the contrary, the viscosity of the less concentrated RP suspensions (i.e. φ = 0.405, 0.45 and 0.47) tends towards a plateau as the shear stress decreases. Even a relatively low but still visible shear-thickening regime can be observed for stresses below typically 1 Pa. The explanation for this shear-thickening is presumably the same as that proposed by Seto et al. 12 to explain DST: the observed continuous shear-thickening should originate from a competition between the repulsive forces resulting from the presence of the surfactant which tends to prevent contact between particles and the forces associated with the flow which, on the contrary, promote frictional contacts between particles. The characteristic shear stress for which shear-thickening takes place is smaller than the one usually reported in studies dedicated to shear-thickening. However, one must keep in mind that the particles used in the present study are as large as 40 µm in diameter. Guy et al. 13 have shown, based on their own data and using the data collected by Barnes in his review 44 , that the critical shearthickening stress varies as a -2 . The values of critical stresses reported by these authors for particles of several tens of microns in size are of the order of 0.1 Pa which is consistent with the findings of Fig. 1. As a last argument supporting that the flow regime is that of a suspension in the thickened state, we would like to stress that the relative viscosity values that we measure above σ ≈ 1 Pa are consistent with those that are measured -or calculated numerically-indifferently in concentrated non-Brownian frictional suspensions 3,7,21,25 or in shear-thickening suspensions beyond the shear-thickening transition 4,13,45 . Finally, note that for the two highest particle volume fraction, no peak viscosity is observed anymore in Fig. 1 for the raw particle suspensions. This finding is consistent with the recent work of Singh et al. 46 who studied how shear-thickening transition can be obscured by the presence of adhesive forces. On the one hand, the onset stress of shear-thickening does not vary with particle volume fraction 1,4,13 while, on the other hand, the yield stress increases with increasing particle volume fraction 22,23 . It is therefore likely that the residual adhesion forces cause a yield stress whose value increases with increasing φ and which, for large enough volume fractions, conceals the shear-thickening transition. Above 5 Pa, both suspensions behave the same. This means in particular that the shear-thinning behavior in the range of stresses above 5 Pa cannot be attributed to adhesive properties of the particles. This suggests, as proposed by Papadopoulou et al. 22 , that shear-thinning may originate from several mechanisms, including adhesion or friction. In the following, we will show how the shear-thinning observed for shear stresses larger than a few Pa -which is almost the same for the two types of suspensions-can be explained by a load-dependent inter-particle friction as proposed by Lobry et al. 25 in the STIVF model. To this aim, we first measure the frictional properties of the TS40 polystyrene particles immersed in the mixture (water +Ucon oil + ZnBr). This will be described in the next section.

AFM measurements

The experiment was performed using Atomic Force Microscope AFM (Bioscope, Bruker). The system is equipped with a liquid cell (DTFML-DD-HE, Bruker) that allows dynamic mode working in liquid environment. A Dynoseed TS40 sphere of radius 20 µm is attached to the end of an AFM rectangular cantilever (NSG-11/tipless, TipsNano, length 130 µm, width 35 ± 3 µm and thickness 2 ± 0.5 µm) using epoxy glue (Araldite)(see Fig. 2.a). Other particles are deposited on a glass substrate coated with a thin layer of glue and a jet of filtered air is used to remove particles that were not stuck on the surface. The particle-seeded substrate is mounted on a 3-axis piezoelectric actuation stage (NanoT series, Mad City Labs) that allows lateral displacement (along y-axis) of

Calibration

Since the vertical and lateral force are deduced from the photodiode signal, it is necessary to determine the conversion factors, k n and k y that enable to convert the photodiode output voltage into force values:

F n = k n V z F t = Q/2a = k l V y /2a = k y V y (4)
where F n and Q are the normal force and the torque applied to the cantilever. V z and V y are the output signals of the photodiode and F t is the tangential frictional force.

Many calibration methods exist (Cleveland method 47 , Sader method 48 , thermal noise method 49 , hydrodynamic method 50 ) and among them we chose the hydrodynamic method. To determine k n , we followed Craig and Neto 50 and calibrated the vertical stiffness of the cantilever by measuring the hydrodynamic drag on the sphere attached to the cantilever while a flat stiff substrate is approaching at constant velocity, U (Fig. 2.c). For a perfectly smooth sphere, in the lubrication approximation (i.e. h 0 << a), a linear relationship between the normal force, F n , and the inverse of the separation distance, h 0 , is expected:

F n = k n V z = 6πηa 2 U h 0 (5)
To account for probe roughness, Eqn. ( 5) is modified by adding to the separation distance a constant, h c , which is related to rough-

)): U/V z = k n (h 0 + h c )/(6πηa 2 ) 6
where h c is a constant free parameter that accounts for probe roughness and for possible error on the determination of the zero of separation. From the fitting of the experimental data, we obtain k n = 105 ± 5 nN/V and h c = 37 nm.

ness morphology (density, extension, height...) 51 :

F n = k n V z = 6πηa 2 U h 0 + h c (6)
In addition, inferring the exact position of the zero of separation raises important difficulties. So, for these two reasons, we estimated that it was appropriate to consider h c as a free parameter. Fig. 3 shows the variation of the approaching velocity divided by the photodiode signal as a function of the tip-substrate separation (see Supplementary Materials, for other calibration curves obtained with an other value of the approaching velocity and with a withdrawing velocity and details on the analysis). A linear variation of 1/V z versus (h 0 + h c ) is observed and a linear fitting provides a measurement of k n (see Eqn 6). The fitting is performed for separation distances lower than typically 1000 nm (region marked by the continuous curve in Fig. 3) and provides k n = 105 nN/m. For higher separation distances, the scatter of the data is rather important but extending the fitting range hardly modifies the fitting results. The uncertainty on the value of k n is dominated by the uncertainty in the viscosity value which is estimated to 5%. Thus we obtain k n = 105 ± 5 nN/V . For the constant h c , we obtained 37 nm which is also the order of magnitude of the roughness height on the probe surface. Note that, for any purpose, we also calibrated the photodiode sensitivity (i.e. the conversion factor between the photodiode output voltage and the cantilever deflection, C n ) with a force curve on a stiff substrate and obtained C n = 0.037 V /nm which gives a spring constant K spring = 3.9 ± 0.1 N/m.

To calibrate the lateral stiffness, we used the method proposed by Ryu and Franck 52 : the substrate is positioned at a given position, h 0 to the colloid probe and is oscillated horizontally: Fig. 4 Hydrodynamic calibration of the lateral stiffness. The probe is immersed in a mixture of water and glycerol with a viscosity of 0.049Pa.s and the substrate is oscillated at a frequency of 120 Hz while the probe is at a distance h 0 . The torque calculated from Eqns. (7,8) and ( 9) is plotted against the amplitude of the photodiode output voltage for four different values of h 0 . The slope of the lines (s) varies with h 0 but tends towards a constant value when h 0 becomes sufficiently small. The value of the slope obtained for the two smallest separation distances is used to determine the lateral conversion factor: k l = 1.30 10 -11 N.m/V . y s = y 0 cos ωt (Fig. 2.d). The horizontal motion of the substrate gives rise to a Couette flow between the flat and the particle which is subjected to both a viscous drag force and a torque that are given by:

F visc = 6πηaU f * τ visc = 8πηa 2 Uτ * (7)
where f * and τ * are functions of h 0 /a. To compute f * and τ * , we use the asymptotic formulas given by Goldman et al. 53 :

f * = - 8 15 ln h 0 a + 0.9588 τ * = 1 10 ln h 0 a + 0.1895 (8)
Assuming that the torque applied on the cantilever itself is negligible, the total torque exerted on the ensemble (cantilever + micro-particle) can be written as:

Q = F visc a + τ visc (9)
Fig. 4 shows the relation between the torque amplitude calculated with Eqns. (7, 8) and ( 9) where U = y 0 ω and the amplitude of the photodiode output voltage, V y0 . It is observed that V y0 increases linearly with Q and k l is given by the slope. The slope varies with h 0 which was expected since the expressions of the force and the torque as a function of the separation distance between the particle and the plane taken from 53 are valid only for small separation distances. Furthermore, we neglected the effect of the cantilever itself which is likely to modify the flow around the colloid probe. The relative contribution of this effect on Q becomes smaller as the separation distance decreases. For these two reasons, the calibration accuracy increases as h 0 decreases. But, on the other hand, the determination of the distance h 0 is quite tricky. Therefore, reducing this distance to too low values is likely to introduce an important error in the evaluation of the torque. h 0 ≈ 6 µm appeared to be a good compromise value and we found k l = 1.30 10 -11 N.m/V .

Friction measurements

To measure the friction force between TS40 particles, the flat substrate is replaced by the particle-seeded substrate whose horizontal position is controlled so as to vertically align one of the particles of the substrate with the particle attached to the cantilever. Then the substrate is vibrated in the y-direction with variable amplitudes (∼ 1 µm) and frequencies between 25 and 100 Hz. It is simultaneously moved vertically at constant speed (∼ 0.1µm/s). An example of force displacement curve is presented in Fig. 5 normal force is the combination of a continuous component that increases as soon as the particles touch (the vertical velocity of the substrate is sufficiently small for the vertical lubricating force to be negligible) and an oscillating force whose period is that of the horizontal vibration of the substrate (see insert of Fig. 5):

F n (t) = F mean n (t) + F osc n cosωt.
The quasi-periodic component of F n comes from geometry and is caused by the curvature of the particles. Fernandez et al. 16 studied in detail the effect of particle curvature on the AFM measurement of normal and tangential forces and showed that it could be very important if the particle size is not very large compared to the oscillation amplitude of the substrate. But, it is not the case in our study where the particles have a radius of 20 µm while the vibration amplitude of the substrate is of the order of one micron, which implies that, during a cycle, the height of the particle probe varies by about δ h ≈ 25 nm to remain in contact with the particle stuck on the substrate, which corresponds to very small angles 2δ h/a. In particular, this evaluation is fully compatible with the amplitude of oscillation of the normal force (F osc n ≈ K spring δ h ≈ 100 nN). We can get rid of curvature effects by considering only the measurements of normal force and tangential force over a sufficiently small interval of y so that the variation of the normal force is small in comparison with its mean value. To this aim, for each quasi-cycle (y, F n ), we consider only the central region whose extension is chosen so that the normal force varies by less than 10% (see blue squares in Fig. 6). In this region, the friction coefficient is defined by (10) where F up t and F down t are the tangential forces recorded when the substrate moves to the right (y > 0) or to the left (y < 0) and the symbol < ... > denotes the average over the chosen interval of y. This is illustrated in Fig. 6 where the typical variations of the normal force and tangential force during a cycle are shown. During a cycle, both the normal force and the absolute value of the tangential force vary. Therefore, to calculate the coefficient, we only select the data for which the normal force deviates from its mean value by lower than 10% (blue squares).

µ = < F up t > -< F down t > 2 F mean n
The friction coefficient is measured for six different particle pairs, for three different frequencies (25, 50 and 100) and different lateral displacement amplitudes (comprised between 600 and 2000 nm). Fig. 7 shows the measured friction coefficient for 12 realizations (small symbols) and the averaged friction coefficient (solid line). It should first be noted that, as expected in the frame of a mono-contact model, the friction coefficient decreases as the normal load increases. Secondly, the measured coefficients of friction vary quite significantly from one realization to another. These differences are not due to experimental uncertainties, but rather to a spatial variation of the local friction coefficient which is expected, especially taking into account the surface roughness. We recall that the characteristic extension of the surface roughness is of the order or slightly less than one micron, which is also the order of magnitude of the lateral displacement on which the friction is probed. It is therefore not surprising that, from one realization to another, the local friction coefficient changes. On the contrary, we verified that neither the frequency nor the amplitude of vibration of the substrate affected the measurement of the friction coefficient. In Fig. 7, some results obtained in quasi-static mode are also represented (large pink stars). These measurements were performed using a cantilever ten times softer (vertical stiffness 0.364 N/m) and vibrating the substrate at low frequency (1 Hz) with an amplitude of 1.5µm. The values of the quasistatic friction coefficients and their evolution with the load are comparable to what is obtained in dynamic mode. However, these quasistatic measurements can hardly lead to reliable quantitative results because they require a rather long duration of experiment during which the normal force is likely to vary due to the drift in time of the photodiode response. This issue may probably explain why the friction coefficient measured for the largest of the applied normal force -and which in the experiment was the last measurement-deviates from the general trend of the variation of µ with F n (see Fig. 7). Furthermore, as seen before, the friction coefficient measured for such small contact forces is essentially a statistical quantity since the result of the measurement depends on the local geometry of the asperities at contact. Therefore, in order to obtain a reliable average µ value, a large number of measurements are necessary. For these reasons, we decided to use the measurements made in quasistatic mode only to validate the dynamic measurements. Fig. 7 shows unambiguously that the inter-particle friction coefficient decreases with the load, contrarily to what is expected for macroscopic bodies whose friction is determined by a constant friction coefficient. Such a load-dependent friction coefficient has been described by Brizmer et al. 32 who considered the contact between a perfectly smooth sphere and a flat. The resulting friction coefficient is given by: µ = 0.27 coth 0.27

F n L c 0.35 (11)
This expression accounts for the transition from elastic deformation of the sphere (L < L c ) to plastic deformation (L > L c ).

In the elastic deformation regime, the Hertz contact law is approximately recovered. In this regime, the contact area varies as F 2/3 n and, since the tangential force is proportional to the contact area, the friction coefficient decreases as F -1/3 n . Eqn. (11) leads to µ ∝ F -0.35 n in this regime, which is consistent. In the plastic regime, the friction coefficient levels off at the value 0.27. L c is a function of the Young modulus, Poisson ratio and Yield stress of the material, and of the contacting sphere radius as well 32 . Coming back to the case of the microspheres whose friction coefficient is measured here, it is likely that, given the magnitude of the normal forces involved in the AFM experiment, the particles come into contact through only one or a few asperities present on their surface § . These correspond to the smooth sphere of the model of Brizmer et al., which consequently predicts a decrease in the friction coefficient with the load.

The friction coefficient computed from Eqn. ( 11) is displayed in Fig. 8, together with the mean friction coefficient from the experiments (symbols). At very small normal load, our results are in fairly good agreement with Brizmer's predictions but then Brizmer model predicts a slower decrease in the friction coefficient than the one we have measured. Also note the discrepancy between the experimentally measured µ values and those predicted by Brizmer et al. for high values of F n . However, these differences are not surprising since Eqn. ( 11) considers an idealized contact between a plane and a hemisphere whereas, in our case, the roughness geometry is probably more complicated. There is therefore no reason why the friction law should strictly obey Brizmer model. Finally, an equation close to that proposed by Brizmer et al. is fitted to the measured mean friction coefficient (see orange line in Fig. 8). The friction coefficient at infinite normal load is taken as 0.18 and L c and n are free parameters:

µ = 0.18 coth 0.18 F n L c n (12)
It is recalled here that L c is to be considered as the critical load that marks the transition from the elastic deformation of one asperity at low load to the plastic deformation at higher load. The best fit (in the range 12) which gives the following fitting parameters: L c = 33.2 nN and n = 0.54. For the sake of comparison the original result of Brizmer et al. 32 is also shown (dashed green line).

Viscosity deduced from friction

In this section, we will test the scenario used in the STIVF model to explain the shear-thinning observed in non-Brownian suspensions by a friction coefficient that varies with load. Let us first recall the main features and results of the model. The authors performed discrete simulations of sheared suspensions, where the friction coefficient between particles obeys Eqn. (11). Suspensions are bidisperse (with particle radii a 1 and a 2 = 1.4a 1 in equal number) in order to avoid cristallization. They found that the variation of the computed viscosity with shear rate could be accurately described using a very simple and physically sound model. The main idea is that the overall shear stress controls the normal force between particles, hence the relevant friction coefficient and consequently the jamming volume fraction of the suspension. More precisely, the viscosity may be computed from the self-consistent relation:

η S = α (µ(Σ)) 1 - φ φ m (µ (Σ)) 2 ( 13
)
where Σ is the normalized shear stress:

Σ = 6πa 2 1 σ 1.69 L c = F n L c (14)
and the relevant friction coefficient is given by Eqn. (11), in the STIVF model. Note that a 1 is 0.8 times the mean radius, so in the following we will take a 1 = 0.8 a = 16 µm.

α(µ) and φ m (µ) are deduced from numerical simulations per-formed for various constant values of µ and are given by 25 :

α(µ) = α ∞ + (α 0 -α ∞ ) exp (-X α atan(µ)) -exp (-πX α /2)) 1 -exp (-πX α /2)) φ m (µ) = φ ∞ m + (φ 0 m -φ ∞ m ) exp (-X p atan(µ)) -exp (-πX p /2)) 1 -exp (-πX p /2)) (15)
with the fitting parameters displayed in Table 1 (first line). The striking feature here is that the relevant friction coefficient in Eqn. ( 13) is actually the microscopic friction coefficient used in the discrete simulations. Therefore the influence of any particular microscopic friction law on the suspension viscosity may be easily probed without the need of further discrete simulations. Using Eqns. ( 13)-( 15), together with the experimental friction law (Eqn. ( 12)) and the values in table 1 (second line), we compute the relative viscosity. As shown in table 1, the values of the rheological parameters from the STIVF model were taken from Lobry et al. 25 , except for φ 0 m and φ ∞ m that we take here equal to 0.65 and 0.55 instead of 0.7 and 0.546, respectively.

The results are presented in Fig. 9 together with the experimental results. Note that the friction coefficient has been measured in the same ranges of force and contact velocity as those explored in the rheometric experiments. In rheometry, the interparticle force range is given by: F = 6π a 2 σ /1.69 ≈ 10 -1500 nN and the relative velocity experienced by the particles by v ≈ a γ ≈ 10 -1500 µm/s. In AFM experiments, normal force has been varied between 10 nN (the estimated resolution) and 1 µN and the relative velocity of particles (v = 2π f y 0 ) between 10 and 1200 µm/s.

The agreement between the experimental results obtained for five different volume fractions and the STIVF model 25 is excellent. This quantitative agreement shows that the rheology of concentrated non-Brownian suspensions is actually controlled by the frictional properties of the particles, as proposed by several authors 10,14,24 . It appears that the numerical simulations of Lobry et al. carried out with a constant friction coefficient have made it possible to propose reliable fitting laws for linking the viscosity to the friction coefficient. Moreover, the good agreement validates many of the hypotheses of the STIVF model and, in particular, the choice to describe the contact between particles by a Brizmer-type single-asperity contact model. The variation of the friction coefficient with normal load measured by AFM for normal forces of the same magnitude as those experienced by particles in sheared suspensions is qualitatively close to that predicted by the Brizmer model, which justifies the use of the latter to predict the shearthinning of concentrated non-Brownian suspensions, as proposed in the STIVF model. Finally, the agreement between the experimental results and the model also highlights the interest and the relevance of the model in Eqn. (13-15), which establishes a direct connection between the suspension viscosity and the microscopic friction coefficient, by means of the relevant shear induced contact force 6πa 2 σ /1.69.

Conclusions

The comparison of the rheology of washed or raw particle suspensions -i.e. in the absence or presence of surfactant -allowed us first to distinguish the respective effects of adhesion and loaddependent friction properties on the shear-thinning of concentrated non-Brownian suspensions. This rapid comparison echoes the detailed study conducted by Papadopoulou et al. 22 on these respective roles and is not the essential result of our paper.

The core of our message is to confirm the scenario proposed by Lobry et al. in the STIVF model that variable friction between particles is at the origin of shear-thinning. To this purpose, we first carried out AFM measurements in order to determine the friction coefficient of the T S40 microspheres and its variation with the normal force in the force range experienced by particles in a sheared suspension. The friction coefficient has been shown to decrease with increasing load, which contrasts to the constant friction coefficient usually observed in the case of contacting macroscopic bodies. This behavior suggests that the contact between particles involves a small number of asperities, as proposed in the STIVF model. The measured friction law slightly differs from the law proposed by Brizmer et al. 32 to describe the friction that occurs between a single smooth hemisphere (the asperity) and a plane. The discrepancy may be explained by the geometry of the surface asperities which is certainly more complicated than the ideal hemisphere considered by Brizmer et al. Nevertheless the experimental friction law displays the same trends as the mono-asperity contact model: at small load, in the elastic deformation regime, the friction coefficient decreases, and it levels off as the normal force increases sufficiently for plastic deformation to occur. The characteristic force L c that marks the transition from a fully elastic to a plastic contact has been found quite close to the value deduced from the roughness size and the mechanical properties (Young modulus and yield strength) of polystyrene. Despite these slight differences, the correlation observed between the decrease in friction coefficient and the decrease in viscosity supports the scenario used in the STIVF model and highlights the close links between the microscopic friction properties of the particles and the macroscopic rheological behavior of the suspension. Furthermore, the quantitative agreement between the predictions of the rheological model where the experimental microscopic friction law has been introduced and the experimental measurements of the viscosity should be stressed. This makes the STIVF model an accurate tool able to quantitatively predict the rheology of moderately concentrated non-Brownian suspensions from the microscopic contact law.

To go further, it would be interesting to reproduce this study with other particles that could be made of a different material or have particular surface properties (electrically charged or grafted with polymer brushes). Finally, in the case of adhesive suspensions, it would be interesting to characterize the adhesion forces using an AFM and to jointly study the rheology of these suspensions in order to establish quantitative links between adhesion and rheology (yield stress, shear-thinning), as it has been done for friction.

a

  Laboratoire Onde et Matière d'Aquitaine, UMR 5798 CNRS-Université Bordeaux, 33405 Talence cedex, France b Institut de Physique de Nice, UMR 7010 CNRS-Université Côte d'Azur, 06108 Nice cedex 2, France; E-mail: elemaire@unice.fr.

Fig. 1

 1 Fig. 1 Variation of the viscosity with shear stress for suspensions made of either washed (hollow squares) or non-washed (stars) TS40 particles.

Fig. 2

 2 Fig. 2 AFM measurement device. a) microsphere glued on a rectangular cantilever. b) sketch of friction measurement. c) and d) sketches of vertical and lateral stiffness calibrations.

Fig. 3

 3 Fig.3Ratio of the approaching velocity (U = 39 µm/s) to the output signal of the photodiode corresponding to the vertical deflection as a function of tip-substrate separation. For this calibration the probe is immersed in a mixture of water and glycerol with a viscosity of 0.050 Pa.s. The conversion factor, k n is obtained by fitting the data in the linear region (dark blue line) with the lubrication equation (Eqn. (6)): U/V z = k n (h 0 + h c )/(6πηa 2 ) where h c is a constant free parameter that accounts for probe roughness and for possible error on the determination of the zero of separation. From the fitting of the experimental data, we obtain k n = 105 ± 5 nN/V and h c = 37 nm.

Fig. 5

 5 Fig.5Normal force as a function of the vertical position of the substrate. The blue curve is the measured force while the orange line is the value of the force averaged over a period of oscillation of the substrate. Insert: zoom of the variation of the normal load with vertical position.

Fig. 6

 6 Fig.6Variation of the normal (up) and tangential (down) force versus the lateral position of the substrate during one period of oscillation. During a cycle, both the normal force and the absolute value of the tangential force vary. Therefore, to calculate the coefficient, we only select the data for which the normal force deviates from its mean value by lower than 10% (blue squares).

Fig. 7

 7 Fig.7The measured friction coefficient for six different particle pairs, different frequencies (25, 50 and 100 Hz) and different lateral displacement amplitudes (ranging from 600 to 2000nm) is depicted against the normal force. The pink stars correspond to quasi-static measurements (frequency 1 Hz, amplitude 1500 nm). The orange line is the mean friction coefficient averaged over all the realizations.

F n ∈ [ 10 ,Fig. 8

 108 Fig.8Variation of the friction coefficient with normal load. The blue squares correspond to the average of all the measurements shown in Fig.7and obtained for different particle pairs, amplitudes and frequencies. The orange line is the result of the data fitting with Eqn. (12) which gives the following fitting parameters: L c = 33.2 nN and n = 0.54. For the sake of comparison the original result of Brizmer et al.32 is also shown (dashed green line).

Fig. 9

 9 Fig.9Normalized viscosity versus shear stress. The symbols represent the experimental measurements performed with the suspensions of nonwashed TS40 microparticles for different particle volume fractions while the orange lines are the results of the model with only two free parameters that are φ 0 m and φ ∞ m , the jamming fractions of suspensions made of frictionless and infinitely frictional particles respectively. Note that only viscosity data obtained for stresses that correspond to the force range in which the friction coefficient has been measured are shown.

Table 1

 1

	fitting parameters
	STIVF model 25	α 0 1	α ∞ 0.64	X α 1.85	φ 0 m 0.70	φ m ∞ 0.546	X p 2.43
	This study	1	0.64	1.85	0.65	0.55	2.43

[1][2][3][4][5][6][7][8][9][10]

1-10 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

This work was supported by the French National Agency (ANR) under the program Blanc AMARHEO (ANR-18-CE06-0009-01). Muhammad Arshad acknowledges financial support from the Higher Education Commission of Pakistan.

Notes and references

Author contributions

Author contributions are defined based on the CRediT (Contributor Roles Taxonomy) and listed alphabetically. Conceptualization: EL, LL, FP. Formal analysis: EL, AM. Funding acquisition: EL. Investigation: MA, CC, EL. Methodology: EL, AM, FP, . Project administration: EL. Supervision: EL, AM. Writing original draft: EL. Writing review and editing: MA, CC, EL, LL, AM, FP.

Conflicts of interest

There are no conflicts to declare.


Mircea Moscu

Ricardo Augusto Borsoi

Cédric Richard

Online Graph Topology Inference with Kernels for Brain Connectivity Estimation

Keywords: Topology inference, reproducing kernel, graph signal processing, adaptive algorithm, brain connectivity estimation

In graph signal processing, there are often settings where the graph topology is not known beforehand and has to be estimated from data. Moreover, some graphs can be dynamic, such as brain activity supported by neurons or brain regions. This paper focuses on estimating in an online and adaptive manner a network structure capturing the non-linear dependencies among streaming graph signals in the form of a possibly directed, adjacency matrix. By projecting data into a higher-or infinite-dimension space, we focus on capturing nonlinear relationships between agents. In order to mitigate the increasing number of data points, we employ kernel dictionaries. Finally, we run a series of tests in order to experimentally illustrate the usefulness of our kernel-based approach on biomedical data, on which we obtain results comparable to state-of-the-art methods.

INTRODUCTION

Graphs represent powerful mathematical objects able to model and analyse any kind of network. Due to their inherently distributed nature, they are also suited for big data analysis, since a distributed solution can be easily applied on such an object. In applications such as gene regulation systems [START_REF] Richard | A second generation human haplotype map of over 3.1 million snps[END_REF], socio-economical interactions [START_REF] Heiberger | Predicting economic growth with stock networks[END_REF], or brain activity [START_REF] Kramer | Emergent network topology at seizure onset in humans[END_REF], information about network structure is paramount. Most graph signal processing algorithms assume the graph structure as known beforehand. Among graph topology inference methods, most assume linear dependencies between the agents (brain regions, genes in a network, sectors of a market economy). Therefore, here resides the need of developing algorithms capable of modeling nonlinear relationships, which are naturally present in real-world applications. Using reproducing kernels, through their ability to model non-linear relationships between nodal signals, is hence motivated.

In this article, we consider a setting where online nodal measurements are acquired and used in order to infer the topology of the underlying network. In the developed approach, it is desired to model non-linear dependencies, to be able to adapt to changes in the network structure, and to be distributed over the different agents, aiming at estimating a possibly directed adjacency matrix. Since many real world examples, such as social graphs, show considerable edge sparsity, the block 1-norm regularisation is used in this paper due to its block sparsity-inducing properties. Definitions: A graph G consists of a set N of N nodes, and a set E of edges such that if nodes m and n are linked, then (m, n) ∈ E. For undirected graphs, these node pairs are unordered. At the node level, we collect a signal y [y1, . . . , yN ] , assumed to be real-valued, where yn is the sample of the signal y at node n. The adjacency matrix A [START_REF] Biggs | Algebraic Graph Theory[END_REF][START_REF] Sandryhaila | Big data analysis with signal processing on graphs: Representation and processing of massive data sets with irregular structure[END_REF], is defined as an N × N matrix whose entries anm are zero if (m, n) / ∈ E and set to one otherwise. Prior works: In topology identification, several works have been developed. An early proposition is in [START_REF] Dempster | Covariance selection[END_REF], where a covariance estimation based method of inferring links is introduced. On the same line, in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical Lasso[END_REF] the graphical Lasso is employed in order to estimate the inverse covariance matrix from data. In [START_REF] Segarra | Network topology inference from spectral templates[END_REF], the authors advocate that connectivity can be recovered from estimated spectral templates, while the authors in [START_REF] Vlaski | Online graph learning from sequential data[END_REF] propose an adaptive algorithm for learning the topology from streaming graph signals driven by a diffusion process. Under a graph signal smoothness assumption, a so-called pairwise distances matrix whose entries are pairwise squared-norm distances between nodal measurements is introduced in [START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF] to estimate a weighted adjacency matrix.

In modeling non-linear phenomena, works such as [START_REF] Harring | A comparison of methods for estimating quadratic effects in nonlinear structural equation models[END_REF][START_REF] Finch | Modeling nonlinear structural equation models: A comparison of the two-stage generalized additive models and the finite mixture structural equation model[END_REF] focus on polynomial structural equation models, while the authors of [START_REF] Lim | Operator-valued kernel-based vector autoregressive models for network inference[END_REF] use their non-linear counterparts. They, however, have some limitations, such as assuming knowledge of certain connections or the form of the non-linear functions. Reproducing kernels have seen widespread use in topology inference problems. One of these works is [START_REF] Shen | Topology inference of directed graphs using nonlinear structural vector autoregressive models[END_REF] where kernels, chosen to best fit the data, model nonlinear relationships between nodes based on measurements at successive time instants. The authors present an auto-regressive framework that allows to track graph connectivity over time, proving useful in providing insights on brain connectivity. The multi-kernel approach in [START_REF] Zhang | Going beyond linear dependencies to unveil connectivity of meshed grids[END_REF] uses partial correlations to encode graph topology and p-norm regression to enhance performance. In [START_REF] Giannakis | Topology Identification and Learning over Graphs: Accounting for Nonlinearities and Dynamics[END_REF], a thorough analysis of the kernel-based topology inference problem is given. This work focuses on capturing non-linear and dynamic links. A review of the state of the art methods for graph topology inference is given in [START_REF] Dong | Learning graphs from data: A signal representation perspective[END_REF].

Encoding the graph topology with the adjacency matrix A is ubiquitous in graph signal models. This operator describes the interactions between entities and, by extension, it can be considered as a tool for representing relationships between data. However, to the best of our knowledge, a kernel-based online solution to the topology inference problem has not yet been considered. In this work, we propose an online approach that can sensibly reduce computational stress, as well as adapt to slow changes in the topology that occur in dynamical environments. For the particular case of brain connectivity estimation, the developed method has the advantage of adaptability concerning data availability: due to its online nature, the data acquisition process can be stopped exactly when the estimate is obtained. This, in turn, can render the medical process of signal acquisition less strenuous for both the patient and medical personnel. Notations: Normal font letters denote scalars, while boldface lowercase and uppercase letters stand for column vectors and matrices, respectively. Operator • 0 returns the number of their non-zero entries. Uppercase calligraphic letters denote sets, of cardinality | • |.

NON-LINEAR MODEL AND CENTRALIZED PROBLEM

Consider an N -node graph with adjacency matrix A which models a system such as the brain network. In this setting, the brain activity in every considered region 1 through N can be measured at different time instants i, thus obtaining a signal y(i). Each of these regions influences and is influenced by the other regions, and these links are encoded in the matrix A. For the particular case of brain connectivity, the existence of non-linear connections have been reported by studies such as [START_REF] Freeman | EEG analysis gives model of neuronal template-matching mechanism for sensory search with olfactory bulb[END_REF][START_REF] De Zwart | Hemodynamic nonlinearities affect BOLD fMRI response timing and amplitude[END_REF]. This motivates modeling certain systems, such as the brain, with non-linear connections. With these remarks, we consider the following data model:

y(i) = Af (i) + v(i), (1)
where A is the adjacency matrix of the graph. This matrix models how entries of f (i) col{fm(y Lm (i))} N m=1 influence every node. We consider fm : R Lm → R a non-linear function whose argument is y Lm (i) = [ym(i), . . . , ym(i -Lm + 1)] . The signal v(i) models innovation noise. The output at every node m is therefore non-linearly dependent of all the other signals from the other nodes, including the past. Given nodal measurements y(i) acquired online, the goal is estimating the adjacency matrix A.

Under the least-mean-squares criterion, the optimization problem can be expressed as:

A * = argmin A 1 2 E y(i) -Af (i) 2 + Ψ(A) subject to anm ∈ {0, 1}, (2)
where Ψ(A) is a regularization term to account for some prior knowledge of A such as symmetry or sparsity. The constraint aims at forcing the entries of the adjacency matrix to be binary. A relaxation of the constraint may be employed, by enforcing anm ∈ [0, 1].

Focusing on a single node n, model (1) becomes

yn(i) = m∈N anmfm(y Lm (i)) + vn(i) (3)
and the local optimization problem is now

a * n = argmin an 1 2 E yn(i) - m∈N anmfm(y Lm (i)) 2 + ψ(an) subject to anm ∈ {0, 1}, (4)
where an is the n th row of A, with entries anm, m = 1, . . . , N . Problem (4) has to be solved based only on local measurements y( ), ≤ i that are available at a certain instant i.

RKHS AND KERNEL DICTIONARIES

Reproducing Kernel Hilbert Spaces

Several solutions in modeling non-linearities exist, such as nonlinear and polynomial Structural Equation Models [START_REF] Jöreskog | Nonlinear structural equation models: The Kenny-Judd model with interaction effects[END_REF], as well as function selection from an existing function set [START_REF] Song | A bayesian modeling approach for generalized semiparametric structural equation models[END_REF]. The focus in this work is on kernel methods, which are able to deal with nonlinearities in classification or regression problems by applying linear algorithms over a high-dimensional representation of the input data on a Reproducing Kernel Hilbert Space (RKHS) H k associated with a positive definite reproducing kernel k(•, •).

For the sake of clarity, consider Lm = 1 from now on, making y1(i) ≡ ym(i). Let us denote φnm = anmfm, which allows us to incorporate the binary variable anm and turn (4) into a problem that is linear in φnm. Assuming that φnm belongs to a RKHS H km , for m ∈ N , and approximating the expected value in (4) by empirical averages computed over the available measurements for ≤ i, a non-parametric version of the local optimization problem for node n and time instant i can be written as:

{φ * nm } N m=1 = argmin φnm∈H km m=1,...,N 1 2i i =1 yn( ) - m∈N φnm(ym( )) 2 + m∈N ψH km ( φnm H km ), (5)
where ψH km : R → [0, ∞[ are non-decreasing functions. Since (5) employs a convex loss function, the conditions of the linear representer theorem are satisfied [START_REF] Schölkopf | A generalized representer theorem[END_REF]. Thus, the solution to (5) admits a finite-dimensional representation of the form:

φ * nm (•) = i p=1 αnmpkm(•, ym(p)), m = 1, . . . , N, (6)
where αnmp ∈ R and {ym(i), . . . , ym(1)} is the set of data available at node m and time instant i.

An important question is how to introduce sparsity in the graph connections now that the problem is formulated in terms of φnm. Since φnm = anmfm, anm = 0 implies that function φnm ≡ 0. Thus, promoting sparsity over A is equivalent to promoting sparsity over the functions φnm, for m = 1, . . . , N . Fortunately, the coefficient-based representation [START_REF] Dempster | Covariance selection[END_REF] means that this can be performed equivalently by promoting sparsity of groups of variables {αnmp} i p=1 , for m ∈ N . This can be done very efficiently by using a block-sparse regularization over the coefficients, which leads to the following optimization problem for each time instant i:

α * n = argmin αn 1 2i i =1 yn( ) -α n k( ) 2 + ηn αn B,1 , (7)
where the N i × 1 block vectors αn and k( ) are defined as:

k( ) = k 1 ( ), . . . , k N ( ) , αn = α n1 , . . . , α nN , (8)
with

km( ) = col{km(ym( ), ym(p))} i p=1 , αnm = col{αnmp} i p=1 . (9)
Constant ηn > 0 is a regularization parameter. Also, a block sparsity-inducing regularisation on αn was added through the term αn B,1 , the block-wise 1-norm, i.e., αn B,1 = N m=1 αnm 2. This norm is known to promote group sparsity [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], favoring solutions with entire blocks of variables αnm equal to zero, from which it can be inferred that anm = 0 and thus there will be no connection from node m towards n.

Optimization and kernel dictionaries

Solving [START_REF] Friedman | Sparse inverse covariance estimation with the graphical Lasso[END_REF] in batch mode is, however, impractical and computationally costly. This is why we propose a stochastic gradient descent-based solution in order to update αn every instant i. Remark that a subgradient of the block-1 regularization αn B,1 = m∈N αnm 2 is given by the block vector [START_REF] Jin | Model-driven online parameter adjustement for zero-attracting LMS[END_REF], where each block Γnm is:

Γn = [Γ n1 , . . . , Γ nN ] [24,
Γnm =    αnm αmn 2 if αmn 2 = 0 0 if αmn 2 = 0 . ( 10
)
This entails the use of group zero-attracting LMS (GZA-LMS) [START_REF] Jin | Adaptive parameters adjustement for group reweighted zero-attracting LMS[END_REF], leading to the following update rule, for every time instant i:

αn(i + 1) = αn(i) + µn[rk y -Rkk αn(i) -ηnΓn(i)], (11)
where rk y = E{ k(i)yn(i)}, Rkk = E{ k(i) k (i)}. Estimating these second order moments can however prove to be unattainable or computationally costly. This warrants for the use of approximations, such as those based on instantaneous realisations:

rk y ≈ yn(i) k(i), Rkk ≈ k(i) k (i). ( 12
)
The use of ( 12) leads to the stochastic GZA-LMS update:

αn(i + 1) = αn(i) + µn k(i)[yn(i) -k (i) αn(i)] -µnηnΓn(i). ( 13
max ω j ∈I i m |km(ym(i), ym(ωj))| ≤ ξm, (14)
where ξm ∈ [0, 1[ determines the level of sparsity and coherence of the dictionary [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. Vector k(i) now only stores functions which satisfy [START_REF] Shen | Topology inference of directed graphs using nonlinear structural vector autoregressive models[END_REF], leading to a rewriting of (9):

k(i) = [k 1 (i), . . . , k N (i)] , km(i) = col{km(ym(i), ym(ωj))} ω j ∈I i m . (15)
It is worth noting that using this approach, every time one kernel function is added to a Dm, all the blocks αnm increase in size by one new entry, ∀n ∈ N . Also, at each instant i, αn and k(i) are of size m∈N |I i m | × 1. Most importantly, the number of entries in the dictionary satisfy |I i m | < ∞ when i → ∞ [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. Algorithm 1 summarizes the developed method. In the last step, τn acts as an edge identification threshold. This parameter is used in order to identify the topology from the estimated coefficients αn(i), determining whether there exist links from each node m ∈ N towards n. When processing real data, τn can be set as to obtain an estimated topology which realistically explains the studied process.

EXPERIMENTS

Multiple experiments have been conducted, with the goal of showcasing different characteristics of the developed algorithm. Firstly, a simple 3-node graph is considered, where non-linearities are present in node interactions. Secondly, on real biomedical data, the obtained results are coherent with the results obtained in other works. Thirdly, we take into consideration a dynamic graph setting in order to test the adaptive capabilities of the algorithm. Along our experiments, we either used the Gaussian kernel k Update αn(i) using ( 13) Set ânm(i) to 1 if αnm(i) ≥ τn, to 0 otherwise Fig. 1. Performance, as well as estimates at i = 50, 150, 250.

Discretized Lorenz Attractor

Consider the discretized version of the Lorenz attractor [START_REF] Lorenz | Computational chaos-a prelude to computational instability[END_REF][START_REF] Tobar | Multikernel least mean square algorithm[END_REF]:

  y i+1 1 y i+1 2 y i+1 3   =   y i 1 y i 2 y i 3   + 0.01    10(y i 2 -y i 1 ) y i 1 (28 -y i 3 ) -y i 2 y i 1 y i 2 - 8 3 y i 3    , i ≥ 0, (16)
where exponents, for ease of notation, denote time instants, with initial conditions [y 0 1 , y 0 2 , y 0 3 ] = [0.01, 0.01, 0.01] . We used the Gaussian kernel and set µn = 0.1, σn = 8, ξn = 0.8. Parameters ηn and τn were set as to achieve the best performance in terms of edge identification rate, EIER Agt -Â 0/(N (N -1))•100%, where Agt is the ground truth, whose binary entries encode direct influence, based on [START_REF] Giannakis | Topology Identification and Learning over Graphs: Accounting for Nonlinearities and Dynamics[END_REF], between node couples, excluding self-loops.

Fig. 1 shows the EIER, and the estimates at iterations 50, 150, 250. These results, obtained after less than 300 samples with the dictionaries containing only between 6 and 8 kernel functions per node, show that the proposed method is able to infer links in a distributed and online manner, even if they are based on non-linear interactions.

Tests on epilepsy seizure data

In this experiment, we aimed to show how the estimated topology using the presented method is consistent with results presented in other works. The used data come from a 39-year-old female subject suffering from intractable epilepsy. The data acquisition and preprocessing information is provided in [START_REF] Kramer | Emergent network topology at seizure onset in humans[END_REF]. The data-set contains 8 instances of electrocorticography (ECoG) time series, each instance representing one seizure and contains voltage measurements from 76 different regions on and inside the brain, during the 10 seconds before the epilepsy seizure (preictal interval) and the first 10 seconds during the seizure (ictal interval). The Gaussian kernel was used, and we set µn = 5 • 10 -5 , σn = 90, ξn = 0.9 for each n. Concerning the choice of the ηn and τn, since there was no ground truth, they were set as to obtain coherent results.

In Fig. 2 we show the estimated connectivity of the brain, during both these intervals, averaged over the 8 instances. Fig. 3 depicts the degree (sum of in-and out-degree), encoded in the radii of the circles, relative to each interval. Interestingly, our online estimate reveals roughly the same behavior before and during the seizure as the estimate obtained using the method developed in [START_REF] Shen | Nonlinear structural vector autoregressive models for inferring effective brain network connectivity[END_REF]. More precisely, the number of total connections decreases from one interval to the other, especially due to the variation of in-degrees for nodes 30 to 50. Further analyzing the connections, nodes 75 and 76 have a small in-degree, however they present a more important outdegree. Observe the decrease of the degree of node 26 or the major increase for node 73. This behavior is consistent with the findings of the aforementioned paper. The algorithm is therefore able to obtain results similar to those obtained in previous works, while based on kernel dictionaries. For reference, the number of kernel functions inserted in the dictionaries varied between 15 to 30, after 4000 samples. These results show how only a reduced of kernel funcare actually needed order to obtain a topology estimate. This fact, alongside the online approach, can translate in reduced computational complexity, depending on the solver, due to the drastically reduced number of needed kernel functions.

Real dynamic setting

The goal of this experiment is to analyze how the proposed method adapts to slight dynamic changes in topology. Once again, we use real data whose details are found in [START_REF] Ford | Did I do that? Abnormal predictive processes in schizophrenia when button pressing to deliver a tone[END_REF]. They represent electroencephalography (EEG) measurements taken from a group of 81 subjects in total, some of which are healthy and some of which suffer from schizophrenia. A simple button-pressing task is set up, in three separate settings where subjects either: 1 -pressed the button and a tone was immediately played, 2 -listened to the tone without the button press, and 3 -pressed the button and the tone was not played. The goal of the experiment was to check how the subjects' brains respond to sensory consequences of their own actions, in healthy and unhealthy subjects. This behavior arises when, for example, one voluntarily moves their eyes from side to side and their brain knows that the environment is not actually shifting. Patients suffering from schizophrenia have difficulties in differentiating between internally and externally generated stimuli. During our experiment, we used the measurements pertaining to three of the healthy subjects, namely subjects 1, 2, 3, and three of the schizophrenia-suffering, namely 67, 68, 69. A total of 5000 measurements per task and per subject were selected and fed to our algorithm, in order: task 1, task 2 and task 3, as if done one after the other. For both the healthy and unhealthy sets we used the exponential kernel with σn = 1, µn = 9 • 10 -2 , ξn = 0.1. In order to obtain comparable results, the same manner of choosing parameters ηn and τn was used in both cases. The estimated topologies, averaged over the three subjects in the healthy and unhealthy groups respectively, are depicted in Fig. 4 for each of the three tasks. For both the healthy subjects (first row) and the unhealthy schizophrenia-suffering (bottom row), no important changes appear in topology while moving from task to task. However, comparing the average healthy and unhealthy subjects, a rather different network structure arises for each case. Some graph metrics are given in Table 1. For each subject in the healthy group, 30 to 80 kernel functions were chosen for the dictionary, while for each subject in the unhealthy group between 20 and 60.

CONCLUSION

An online, kernel-based, and distributive graph topology inference method was devised, which advocates the use of kernel dictionaries as a sparsifying solution. The use of kernels allows for the inference of connections when non-linear links are presumed. While most state of the art methods rely on batch methods, which usually come with a high computational cost, the developed online algorithm comes with advantages such as adaptability. On the considered biomedical data, the method proved effective, paving the way to further work and research, such as a thorough theoretical analysis and the use of a multi-kernel approach.

  ) An immediate observation concerning update (13) is that the size of k(i) can become prohibitive as i increases, since each acquired measurement increases the number of kernel functions. A solution to this problem are kernel dictionaries which admit a new candidate kernel function only if the candidate function passes a certain sparsification rule. Under this framework, each node m in the network creates, updates, and stores a dictionary of kernel functions, Dm = {km(•, ym(ωj)) : ωj ∈ I i m ⊂ {1, . . . , i -1}}, where I i m represents the set of time indices of elements selected for the dictionary, before instant i. A candidate kernel function km(•, ym(i)) is added in Dm if the following sparsification condition holds [26]:

  G n (a, b) = exp -ab 2 /2σ 2 n or the exponential kernel k E n (a, b) = exp (-ab /σn).

Algorithm 1 :

 1 Kernel-based online topology inference For every node n: Inputs: Parameters µn, ηn, kn(•, •), ξn, and τn Initialization: Set all entries of αn(0) to 0 Algorithm: At each time instant i ≥ 1 Introduce kn(•, yn(i)) in Dn if (14) is verified Receive and store the updated dictionaries Dm, m = n Compute k(i) defined in[START_REF] Zhang | Going beyond linear dependencies to unveil connectivity of meshed grids[END_REF]

  (a) Preictal interval (b) Ictal interval

Fig. 2 .

 2 Fig. 2. Estimated adjacency matrices for each interval

Fig. 3 .

 3 Fig.3. Summed in-and out-degrees for the estimated graphs for both preictal and ictal intervals. The larger the radius corresponding to node n, the larger the summed degree of node n.

Fig. 4 .

 4 Fig. 4. Estimated topologies per task, averaged per group.

Table 1 .

 1 Metrics for the estimated topologies concerning the average healthy and unhealthy subject

	Metric	Task	Healthy Unhealthy
		1	0.356	0.410
	Network density	2	0.286	0.451
		3	0.310	0.395
		1	22.41	25.84
	Average in-(and out-)degree	2	17.98	28.44
		3	19.50	24.86