# Heterogeneous Multiscale Multivariate Autoregressive Model: existence, sparse estimation and application to functional connectivity in neuroscience

Stefano Spaziani, Gabrielle Girardeau, Ingrid Bethus, Patricia Reynaud-Bouret

<figure><img src="image_1.png" /></figure>

## To cite this version:

Stefano Spaziani, Gabrielle Girardeau, Ingrid Bethus, Patricia Reynaud-Bouret. Heterogeneous Multiscale Multivariate Autoregressive Model: existence, sparse estimation and application to functional connectivity in neuroscience. Journal of Mathematical Biology, 2025, 90 (6), pp.62. 10.1007/s00285-025-02224-x. hal-05355797

HAL Id: hal-05355797

https://hal.science/hal-05355797v1

Submitted on 9 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Heterogeneous Multiscale Multivariate
Autoregressive Model: existence, sparse
estimation and application to functional
connectivity in neuroscience

Stefano Spaziani<fnref n="1" />, Gabrielle Girardeau<fnref n="2" />,<fnref n="3" />, Ingrid Bethus<fnref n="4" />,
Patricia Reynaud-Bouret<fnref n="1" />*

*Corresponding author(s). E-mail(s):

Patricia.Reynaud-Bouret@univ-cotedazur.fr;

Contributing authors: Stefano.Spaziani@univ-cotedazur.fr;

gabrielle.girardeau@inserm.fr; Ingrid.Bethus@univ-cotedazur.fr;

Abstract

In neuroscience, functional connectivity is a concept that can be mathematically
formulated as a graph of interactions between oscillatory brain rhythms and indi-
vidual neuronal activity. Typically, each graph corresponds to a cognitive state
and provides insight into higher cognitive processes, such as learning. However, a
model or method to assess directed interactions within and across brain oscilla-
tions and individual neuronal activity is lacking. In this article, we propose a new
model called HM-MVAR (Heterogeneous Multiscale Multivariate Autoregressive)
that represents linear combinations of neural interaction patterns. These patterns
include phase-locking, where the amplitude of one signal is modulated by the
phase of another, and power-triggered phenomena, where the amplitude of one

<footnote n="1">¹Université Côte d'Azur, CNRS, LJAD, 28 Avenue Valrose, Nice, 06100,
France.</footnote>

<footnote n="2">²Sorbonne Université, CNRS, Inserm, Center of Neuroscience Neuro-SU,
F-75005, Paris, France.</footnote>

<footnote n="3">³Sorbonne Université, CNRS, Inserm, Institut de Biologie Paris-Seine,
IBPS, F-75005, Paris, France.</footnote>

<footnote n="4">⁴Université Côte d'Azur, CNRS, IPMC, 660 Rte des Lucioles, Valbonne,
06560, France.</footnote>
signal is modulated by the amplitude of another. Due to the multiscale struc-
ture, we use a block version of stationarity to exhibit conditions under which the
process obeying the HM-MVAR model exists and is stationary. We also propose
a data-driven weighted LASSO estimator based on martingale exponential devi-
ation inequalities, which may have intrinsic interest in probability theory. We
prove that our estimator satisfies an oracle inequality and present its strong per-
formance in realistic simulations. Finally, applying our model and method on a
publicly available multiscale electrophysiological data set, we recover interactions
described in the literature but also uncover new phenomena of potential interest.

**Keywords:** Autoregressive process, Hawkes process, Multiscale approach, Functional connectivity, Wavelet

MSC Classification: 62M10, 62P10, 62G05, 62H12, 60G10

# 1 Introduction

The brain continuously processes various external and internal stimuli to generate and adapt behavior. These processes are captured through electrophysiological recordings from different brain areas in freely moving animals. The raw recordings consist of continuous signals sampled at a high frequency. From these recordings, two types of data are extracted and analyzed: spike trains and local field potentials (LFPs). Typically, multiple electrodes are used, allowing for the capture of one LFP along with a collection of spike trains for each recorded brain region.

Spike trains represent the times at which individual neurons emit action potentials. Usually, there are tens to hundreds of spike trains, each of them corresponding to a neuron in close proximity to the electrode. Spike trains are encoded as sequences of 0s and 1s -also called point processes in discrete time- with a '1' indicating the occurrence of an action potential of a given neuron at a given time. In contrast, a local field potential is a continuous compound signal that emerges primarily from the summed activity of synapses in a larger area surrounding the electrode (Buzsáki, Gy and Anastassiou, Costas A and Koch, Christof). The LFP is represented as a real-valued time series, with the signal typically sampled at around 1,000 Hz.

LFPs are similar to electroencephalogram (EEG) signals; however, while tradi-
tional EEGs are recorded from the scalp, typically in humans, LFPs are captured from
within the brain using intracranial electrodes, often in rodents. Both LFP and EEG
signals can be decomposed into various frequency components associated with differ-
ent brain states and cognitive functions. These brain rhythms, such as hippocampal
theta (6–12 Hz) and cortical gamma (30–150 Hz), are believed to play a crucial role
in synchronizing neural activity within and across brain areas to efficiently process
information (Buzsáki, 2009).

The interactions within and between these two types of signals—spike trains and
local field potentials (LFP)—are believed to underlie cognitive processes in the brain.
For example, interactions between the spike trains of different neurons are involved
in plastic processes related to learning (Kobayashi et al., 2019). Interactions between
rhythms of higher and lower frequencies, such as the nesting of gamma oscillations within specific phases of theta, support the selective transfer of information between brain areas (Buzsáki and Wang, 2012; Lisman and Jensen, 2013). Additionally, oscillatory patterns can trigger significant changes in the spiking activity of neurons, such as the increased firing rate observed during large gamma oscillations or fast ripples in the hippocampus (Jensen, 2005). Conversely, spiking activity can also modulate brain rhythms (Teleńczuk et al., 2017).

All these interactions—Spike-Spike, Spike-LFP, and LFP-LFP—are not constant over time and often emerge in association with specific events (e.g., stimuli or behaviors) or brain states (e.g., attention or sleep). The association between a particular interaction pattern and a cognitive state is referred to as "functional connectivity" in neuroscience, representing a unique and statistics-oriented approach to understanding brain operations (Lisman and Buzsáki, 2008). While various tools exist for estimating interactions among homogeneous data—such as spike trains on one side or LFP/EEG rhythms on the other side—the interplay between spikes and rhythms remains more elusive. Most studies have focused on specific pairs of data types and patterns, such as the phase-locking of neuronal activity with gamma rhythms, deemed relevant to the cognitive processes under investigation.

To move beyond a pairwise view, the community needs a versatile tool capable of handling heterogeneous data, jointly analyzing all spike trains and LFP rhythms across different brain regions. This tool should produce a visualization in the form of an interaction graph that encompasses all recorded variables. To achieve this, a mathematical model that corresponds to the graph of interactions is essential. The purpose of this article is to provide such a model and the associated estimation method, enabling a more precise understanding of the dynamics from a holistic perspective.

Graphs of dependence between heterogeneous data have become increasingly common in statistics (De Leon and Chough, 2013). However, these models do not account for temporal evolution or frequency decomposition and are therefore poorly suited for capturing interactions between spike trains and LFP rhythms.

Let us focus on homogeneous data. For point processes, various works aimed at estimating graphs of interactions, even when these graphs are not constant over time (Kreiss et al., 2019). Regarding models of spike trains, extensive research has been conducted on Hawkes processes and their variants, such as the Galves-Löcherbach model, to estimate functional connectivity graphs (Galves and Löcherbach, 2015; Lambert et al., 2018). Conditions are established to determine when the model is stationary or not (Galves and Löcherbach, 2015; Brémaud and Massoulié, 1996), and sparse estimation procedures exist for graph reconstruction in this context (Lambert et al., 2018). Additionally, tools are available to test for the existence of interactions within the graph (Carstensen et al., 2010). It has also been proved that Lasso estimators satisfy an oracle inequality (Hansen et al., 2015; Ost and Reynaud-Bouret, 2018). The present work aims to extend this framework to the more general case of heterogeneous data, where spike trains are combined with time series data, such as LFP or EEG.

When only LFP or EEG data are considered, multivariate autoregressive regression (MVAR) models are extensively used and developed (Murtagh et al., 2004; Baranowski and Fryzlewicz, 2020; Kwon et al., 2007; Zhang et al., 2017). In particular
the conditions for stationarity are known (Baranowski and Fryzlewicz, 2020) and the reconstruction of the functional connectivity is often done via Partial Directed Coherence (PDC) and tests for the existence of an interaction (Baccala and Sameshima, 2001; Faes et al., 2012). In (Baranowski and Fryzlewicz, 2020), the convergence of the Least Squares estimator is proved; however, to our knowledge, no oracle inequality or sparse estimator for the graph is currently available. Moreover, when examining rhythms in LFP/EEG data, the existing PDC approach can analyze each rhythm individually but fails to identify interactions between different rhythms directly. Some attempts assess coherence or partial coherence; however, the interpretability of these results is acknowledged to be limited (Albo et al., 2004) (see also (Vinck et al., 2011; Zhu and Wang, 2015) for an application between rhythms and spike trains).

Interactions between spike trains and LFP are also studied, but without defining a priori rhythms of interest. In (Gong et al., 2019), the authors present a model-free framework in which interactions are detected through non-parametric tests using Granger causality. In (Lepage et al., 2013), the LFP can influence the spike trains but it is considered only as exogenous data and is not modeled. Lepage et al. (2013) prove the consistency of the least square estimator for spiking probabilities, as well as the consistency for the variance estimator for spike-LFP modulation. In (Wang and Shanechi, 2019), a similar approach is presented, but the maximal interaction lag is a supplementary parameter to be estimated. In (Abbaspourazad et al., 2019), with a similar model, the parameters are fitted through an unsupervised learning algorithm. Lastly, in (Song et al., 2022), the model also includes a latent brain state variable.

Let us also mention works that mix heterogeneous data (Dion and Lemler, 2020; Bonnet et al., 2022), where spike trains—modeled as a Hawkes process—directly influence the voltage of a given neuron. The primary purpose of these studies is to adaptively estimate the functions that parameterize the voltage itself, rather than to explore the full range of possible interactions among the data.

The model we propose here, called Heterogeneous Multiscale MultiVariate AutoRegressive model (HM-MVAR), combines point processes and wavelet coefficients derived from time series. It relies on multiple sort of interaction functions: from wavelet coefficients to wavelet coefficients, from spike train to spike train, from wavelet coefficients to spike trains, and from spike trains to wavelet coefficients.

In the literature, numerous models exist that capture the dynamics of discrete-time wavelet coefficients. In particular, Markov models are used, especially for the dependency of one rhythm on another, but only the directed interactions from the coarser wavelet coefficient to the finest wavelet coefficients were considered (Willsky, 2002; Frese et al., 2002; Zhang et al., 2004). In (Zhang et al., 2004) the convergence rate of the Kalman Filter Least Squares estimator is proved, along with the corresponding algorithm complexity. However, in all of these works, the relationships between the variables are always vertical, generally from coarser to finer scales and time is not taken into account, which means the evolution over time of the process is not modeled.

Thus, to the best of our knowledge, a model able to estimate directed interactions in time between wavelet coefficients modeling the LFP rhythms and point processes modeling the spike trains is still lacking. Here, (i) we propose a generative probabilistic model, called HM-MVAR, which accounts for various types of interactions between
all LFP rhythms and spike trains. In our framework, an interaction corresponds to a nonzero function between a source and a target—regardless of whether they are spike trains or LFP rhythms—that quantifies the impact of the source on the behavior of the target as a function of the time lag between both quantities. Therefore the interaction is directed and has a causal interpretation in Granger sense (Granger, 1980). (ii) We show that such a process exists and has a block stationary version under some constraints. (iii) We propose a Lasso estimator of the interaction functions that satisfies an oracle inequality and therefore provides a sparse directed graph, which is supported by theoretical guarantees, to summarize the interactions between the data. And finally, (iv) we show on simulation and on a representative open data set (Petersen et al., 2020; Girardeau et al., 2017b) that our method is able to provide a meaningful heterogeneous functional connectivity graph in neuroscience.

Throughout the article, the mathematical framework is accompanied by "Practical perspective" paragraphs, to explain our modeling choices and to provide concrete examples. In Section 2, we present the model equations and the variety of phenomena that it can reproduce. In Section 3, we prove the existence of a block stationary solution. In Section 4, we propose a Lasso estimator of the interaction functions, with weights issued from martingales concentration results that may have an interest per se in mathematics. We also prove that this Lasso estimator satisfies an oracle inequality, which provides theoretical conditions for retrieving the best network that describes the observed activity. In Section 5, we present a simulation study and in Section 6, we apply the method on real data and discuss the relevance of the results.

# 2 The HM-MVAR model

## 2.1 Notation

We are focusing on a discrete time setting where the time index $t$ lies in $\mathbb{Z}$. We use the notation $a : b$ to represent all integers $i$ such that $a \le i \le b$.

We want to model the interactions of $P$ time series, $(X_t^p)_{p \in 1:P, t \in \mathbb{Z}} \in \mathbb{R}^{P \times \mathbb{Z}}$ (e.g. the local field potentials of $P$ different brain areas) and $M$ discrete-time point processes, $(S_t^m)_{m \in 1:M, t \in \mathbb{Z}} \in \{0, 1\}^{M \times \mathbb{Z}}$ (e.g. spike trains of $M$ different neurons).

**Practical perspective 1 (Process notation).** Mathematically speaking, we investigate stationarity properties, so one needs to assume that our processes $(X_t^p)_{p \in 1:P, t \in \mathbb{Z}}$ and $(S_t^m)_{m \in 1:M, t \in \mathbb{Z}}$ are defined from $-\infty$ to $+\infty$. This is why the processes are indexed by $\mathbb{Z}$ even if in practice we observe them only on a finite interval (see Section 4.2).

Also to simplify mathematical formulas, the time index $t = 0, 1, ...$ has no unit. In practice the signal is observed with a certain sampling frequency $f_S$ in Hz (for instance 1250 Hz) and time (index) $t$ correspond to the real time $t/f_S$ seconds.

We consider a discrete, eventually bi-orthogonal wavelet basis $\phi_{j,k}(\cdot)$, where $j \in 0 : J$ is the scale parameter and $k \in \mathbb{Z}$ the localisation parameter (see (Cohen et al., 1992) for a continuous time framework for wavelets and (Gubner and Chang, 1995) for a discrete time one). More precisely, the scale $j = 0$ corresponds to the father wavelet with $\phi_{0,k}(t) = \phi_{0,0}(t - k2^J)$. For $j \ge 1$, we have that the mother wavelets satisfy
$$ \phi_{j,k}(t) = \phi_{j,0}(t - k2^{J+1-j}), \quad \forall t \in \mathbb{Z}, k \in \mathbb{Z}, j \in 1:J. \qquad (1) $$

We assume that we have the following wavelet decomposition of the time series:

$$ \forall t \in \mathbb{Z}, \forall p \in 1: P, \quad X_t^p = \sum_{j=0}^{J} \sum_{k \in \mathbb{Z}} W_{j,k}^p \phi_{j,k}(t). \qquad (2) $$

We also assume in the sequel that the bi-orthogonal wavelet basis $\phi_{j,k}$ has a finite support that we denote $[\ell_{j,k}, r_{j,k}]$ and because of the translation properties before, we always have that

$$ \forall k \in \mathbb{Z}, \quad \left\{ \begin{array}{ll} \ell_{j,k} = \ell_{j,0} + (k-1)2^{J+1-j} & \text{if } j \ge 1 \\ \ell_{0,k} = \ell_{0,0} + (k-1)2^J & \text{if } j = 0. \end{array} \right. $$

and $r_{j,k} = l_{j,k+1} - 1$ for $j \ge 0$ and $k \in \mathbb{Z}$.

**Practical perspective 2 (General wavelet filters).** Wavelet functions $\phi_{j,k}$ are time filters that are applied on the time series $(X_t^p)_{t \in \mathbb{Z}}$ that correspond to the recording of a brain area p. Each filter is indexed by two parameters, the scale j and the localization k. The scalar product in time between the wavelet filter $\phi_{j,k}$ and the signal results in the wavelet coefficient $W_{j,k}^p$:

$$ W_{j,k}^p = \sum_{t \in \mathbb{Z}} X_t^p \phi_{j,k}(t). \qquad (3) $$

In the case of orthogonal wavelets, the same filters are also used in the wavelet decomposition formula (see (2)). The biorthogonal wavelet set-up is the most general set-up where the transform and anti-transform filters may differ. Whatever the case, these filters $\phi_{j,k}$ are obtained by translation and approximate dilation of the mother wavelet $\phi_{1,1}$ (see (1) for an example), which leads to a physical interpretation of change in frequency as the scale j changes.

Remark that because of this interpretation of the wavelet coefficients in terms of application of local filters at a certain frequency, one can consider

$$ X_{t,j}^p = \sum_{k \in \mathbb{Z}} W_{j,k}^p \phi_{j,k}(t), \qquad (4) $$

as the filtered signal $X_t$ on the frequency corresponding to the scale j.

Finally, the father wavelet $\phi_{0,k}$ filter gives the basic trend of the signal at the coarser scale. This takes into account the fact that the signal may have a certain drift through time.

**Practical perspective 3 (The particular example of the Haar filter).** A typical example of such a wavelet basis is the Haar basis, illustrated in Figure 1, for different values of j and k. For this basis, we have that for all $j \ge 1, k \in \mathbb{Z}$,
$$ \phi_{0,k} = \frac{1}{2^{J/2}} \mathbf{1}_{t \in [1+(k-1)2^J, k2^J]} \quad \text{and} $$

$$ \phi_{j,k}(t) = \frac{1}{2^{(J-j+1)/2}} \mathbf{1}_{t \in [1+(k-1)2^{J-j+1}, k2^{J-j}]} - \frac{1}{2^{(J-j+1)/2}} \mathbf{1}_{t \in [1+k2^{J-j}, k2^{J-j+1}]} $$

In particular, with the notation above, $\phi_{j,k}$ has a finite support on $[\ell_{j,k}, r_{j,k}]$ with $\ell_{j,k} = 1 + (k-1)2^{J-j+1}$ and $r_{j,k} = k2^{J-j+1}$.

The shape of the Haar oscillation - upper phase in the first half of the support and lower phase in the second half - results in the following biological interpretations:

* Due to the finite and non-overlapping supports of the wavelets $\phi_{j,k}$, the corresponding wavelet coefficient $W_{j,k}^p$ is the signed strength of the oscillation starting at point $\ell_{j,k}$ and ending at point $r_{j,k}$;

* Since there are $f_S/2^{J+1-j}$ such oscillations in 1s, this corresponds to the frequency $f_S2^{j-1-j}$ Hz. These frequencies can be put in correspondence with well-known rhythms in the neuroscience literature (see Table 1);

* In practical applications, given a sampling frequency $f_S$ Hz, $W_{j,k}^p$ is the signed strength at $\frac{f_S}{2^{j-1}}$ Hz between times $\frac{\ell_{j,k}}{f_S}$ and $\frac{r_{j,k}}{f_S}$ s.

<figure><img src="image_3.png" /><figcaption>Fig. 1 Illustration of the wavelet filtering with the Haar basis. Left: in red, some Haar wavelet filters $\phi_{j,k}$ represented with their starting and ending points $\ell_{j,k}$ and $r_{j,k}$. The x axes represent the discrete time indices in $\mathbb{Z}$, here shown starting from 1 to 64 and the related physical times in seconds with a sampling frequency of 1250Hz. Right: Wavelet decomposition of a signal composed by a superposition of three different rhythms. In the upper figure the signal is shown, again with the discrete and physical time x axes, and in the lower figure, a heatmap represents the values of the wavelet coefficients obtained through (3). As j increases, the temporal resolution and the corresponding frequency interpretation increase. Strong values of the coefficient correspond to the position of the strong oscillations at different frequencies from the original signal. Note also that the larger j, the more coefficients we have in the same interval.</figcaption></figure>
<table><caption>**Table 1** Discrete set of frequencies employed in the simulations (Section 5) and in the analysis of experimental data (Section 6), corresponding to a maximum scale $J = 10$ and to an LFP sampling frequency of 1250Hz. The names of the frequencies (low theta $\theta$, high theta $\Theta$ etc) broadly follow the known literature (Buzsáki, 2009).</caption><thead><tr><th>Rhythm name</th><th>Slow 1</th><th>δ</th><th>Δ</th><th>θ</th><th>Θ</th><th>β</th><th>γ</th><th>Γ</th><th>Ripples</th><th>N/D</th></tr></thead><tbody><tr><th>Scaling (j)</th><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td></tr><tr><th>Frequency (Hz)</th><td>1.22</td><td>2.44</td><td>4.88</td><td>9.76</td><td>19.53</td><td>39.06</td><td>78.12</td><td>156.25</td><td>312.5</td><td>625</td></tr><tr><th>Graph Label</th><td>S1</td><td>de</td><td>DE</td><td>th</td><td>TH</td><td>be</td><td>ga</td><td>GA</td><td>RI</td><td>N/D</td></tr></tbody></table>

As seen above in Practical perspective 3, we can now say more generally that the coefficient $W_{j,k}^p$ appears at time location $\ell_{j,k}$, because it represents the strength of the oscillation that started at time $\ell_{j,k}$, and that it represents the oscillatory power at that time location, at a given frequency.

Our model describes the evolution of $(X_t^p)_{p \in 1:P, t \in \mathbb{Z}}$ and $(S_t^m)_{m \in 1:M, t \in \mathbb{Z}}$ by describing directly the evolution of the coefficients themselves $(W_{j,k}^p)_{p,j}$ as a function of $k$, jointly with the evolution of $(S_t^m)_{m=1,...,M}$, as a function of $t$. In particular, since $W_{j,k}^p$ represents the strength of the oscillation of the $p$th LFP at scale $j$ and time location $\ell_{j,k}$, we want to describe the co-evolution of processes that happen at different time scales. From a modeling point of view, we want to keep $X_t^p$ adapted to the filtration $\mathcal{F}_t$ that we use. This is done by defining the following filtration,

$$ \forall t \in \mathbb{Z}, \quad \mathcal{F}_t = \sigma \left( (S_u^m)_{\substack{m \in 1:M \\ u \le t}}, (W_{j,k}^p)_{\substack{p \in 1:P \\ j \in 0:J \\ \ell_{j,k} \le t}} \right), $$

so that with (2), we can compute $X_t^p$ thanks to all the variables whose values are revealed before time $t$.

**Practical perspective 4.** (Finite support of the model) Even if infinite support wavelet like the Morlet wavelet (Cohen, 2019) is used since a long time to represent EEG signals, here we need a finite support of the basis to ensure interpretability of the resulting reconstruction. Starting from an initial condition, our model is going to produce the spike times as well as the strength of the oscillation at scale $j$ and time location $\ell_{j,k}$, step by step and as a function of the past. This means reciprocally that the LFP signals $(X_t^p)_{p \in 1:P}$ can be built step by step by progressively discovering the corresponding coefficients through time and using (2). If we authorize infinite support for the wavelet, then one would need to observe all the past and future coefficients at once before being able to even assess the value of $(X_t^p)_{p \in 1:P}$ for $t = 1$ and this would therefore prevent us to access any form of causality interpretation between the data of interest. Note that in Granger sense, causality is defined as a significative statistical dependence between two variables with the source event happening before the target event (Ding et al., 2006).

## 2.2 Evolution of the point process part

The model of evolution of the point processes part, i.e. the spiking probability of the single neurons, is inspired by discrete Hawkes processes (Ost and Reynaud-Bouret,
2018), where we incorporate the dependency on the LFP wavelet coefficients. This leads to this set of equations for the point process part $(S_t^m)_{m \in 1:M, t \in \mathbb{Z}}$:

$$ \mathbb{P}(S_t^m = 1 | \mathcal{F}_{t-1}) = \nu^m + \sum_{\substack{m' \in 1:M \\ u < t}} h_{m'}^m(t-u)S_u^{m'} + \sum_{\substack{p \in 1:P \\ \kappa_{pj}^S \in \mathbb{K}_{pj}^S \\ j, k/\ell_{j,k} < t}} h_{p,j,\kappa_{pj}^S}^m(t-\ell_{j,k})\kappa_{pj}^S(W_{j,k}^p) \\ = \nu^m + \sum_{\substack{m' \in 1:M \\ d \in 1: +\infty}} h_{m'}^m(d)S_{t-d}^{m'} + \sum_{\substack{p \in 1:P \\ \kappa_{pj}^S \in \mathbb{K}_{pj}^S \\ j, k/\ell_{j,k} < t}} h_{p,j,\kappa_{pj}^S}^m(t-\ell_{j,k})\kappa_{pj}^S(W_{j,k}^p) \quad (5) $$

In the previous equation, the interaction function $h_{m'}^m$ represents the interaction of spike train number $m'$ on spike train number $m$, whereas $h_{p,j,\kappa_{pj}}^m$ represents the interaction of the brain rhythm of scale $j$ in region $p$, modulated by the mask function $\kappa_{pj}^S$ on spike train number $m$. By forcing $h_{p,j,\kappa_{pj}}^m = 0$ we recover the classical Hawkes process equation (see (Ost and Reynaud-Bouret, 2018)). Indeed, the spontaneous part, $\nu^m$, corresponds to the baseline firing rate of neuron $m$ if no interaction takes place. Moreover if $h_{m'}^m(d) > 0$, neuron $m'$ excites neuron $m$ after a delay $d$, whereas if $h_{m'}^m(d) < 0$, neuron $m'$ inhibits neuron $m$ after a delay $d$. An example of such interactions is provided in the following example.

**Example 1** (Excitation and inhibition of one neuron by another neuron). Figure 2 shows the spike trains of two neurons. If the two neurons are independent, they exhibit a Poisson behaviour which is modulated only by their spontaneous firing rate $\nu^m$. The excitation case corresponds to a positive $h_2^1$ at short range: shortly after a spike generated by Neuron 2, there are more spikes produced by Neuron 1. The modeling of inhibition corresponds to a negative $h_2^1$, so it is less probable that Neuron 1 spikes shortly after a spike from Neuron 2.

The LFP wavelet contribution in the above equation (5) is more intricate. First, (see below) the $W_{j,k}^p$'s are in general unbounded, so we cannot make use of a linear transformation of $W_{j,k}^p$ in (5), as we could possibly end up with a probability of spiking which is negative or greater than 1. This is why we use a mask $\kappa_{j,k}^S(W_{j,k}^p)$ with $\kappa_{j,k}^S(0) = 0$, assumed to be 1-Lipschitz and bounded by two absolute constants $\underline{w}_{p,j} \le 0$ and $\bar{w}_{p,j} \ge 0$ so that

$$ \underline{w}_{p,j} \le \kappa_{pj}^S(W_{j,k}^m) \le \bar{w}_{p,j}. $$

A typical example is (with fixed $w > 0$)

$$ \kappa_{pj}^S(W_{j,k}^p) = \kappa_w^S(W_{j,k}^p) := \begin{cases} w & \text{if } W_{j,k}^p > w \\ W_{j,k}^p & \text{if } -w \le W_{j,k}^p \le w \\ -w & \text{if } W_{j,k}^p < -w \end{cases} $$
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Excitation and inhibition of one neuron by another one. The time occurrences of spikes are represented by dots (blue for Neuron 2, black for Neuron 1). The red stripes of length 250 which appear after the spikes of Neuron 2 represent the time interval for which the spiking probability of Neuron 1 is modified if $h_2^1(t)$ is nonzero. (Top): Independent neurons case, so $h_2^1(t) = 0$. The parameters characterizing the two spiking probabilities (see (5)) are: $v^1 = v^2 = 2e - 3$. (Bottom left): Excitation case. $h_2^1(t) = 0.011_{t \in [1,250]}$. Some new spikes for Neuron 1 appear inside the red stripes, thus showing the excitatory effect of Neuron 2 in action. (Bottom right): Inhibition case. In this case $h_2^1(t) = -0.011_{t \in [1,250]}$. The negative interaction lowers the spiking probability of Neuron 1 inside the red bands, thus Neuron 2 inhibits Neuron 1. Simulations are performed with the same background noise (fixed random seed).</figcaption></figure>

Because we can envision many of these masks, one might want to incorporate several of them to mimic different LFP→Spikes phenomena. This is why we can use a finite ensemble $K_{pj}^S$ of at most K masks, that might be tuned for each $(p, j)$, thus for every brain area $p$ and scale $j$.

For each of these masks $\kappa_{pj}^S(W_{j,k}^p)$ and each $(p, j)$, one has an interaction function $h_{p,j,\kappa_{pj}^S(t-\ell_{j,k})}^m$. Differently from the interaction functions $h_{m'}$, its impact is not just modulated by a 0/1 phenomenon but by the value of $\kappa_{pj}^S(W_{j,k}^p)$. Moreover, the time lag between the source and the target variables cannot take arbitrary values, since the coefficient $W_{j,k}^p$ appears at time location $\ell_{j,k}$.

To fix ideas, one can envision the following two examples.

**Example 2** (Phase lock of spikes on a given rhythm). Let us define two masks $\kappa_+(W) = |\kappa_w^S(W)|1_{W>0}$ and $\kappa_-(W) = |\kappa_w^S(W)|1_{W<0}$. If we use the Haar wavelet
basis, one can take for a fixed positive parameter a,

$$h_{p,j,\kappa_{+}}^{m} = a \mathbb{1}_{1:2^{J-j-1}} \text{ and } h_{p,j,\kappa_{-}}^{m} = a \mathbb{1}_{2^{J-j-1}:2^{J+1-j-1}}.$$

This models a phase locking phenomenon, where spikes of m are more likely to appear in the upper oscillation of the jth rhythm of the pth LFP (see Figure 3). Indeed if $W_{j,k}^p > 0$, the oscillation of length $2^{J+1-j}$ presents the up phase followed by the down phase, so in the first half of the period, the neuron is excited. If $W_{j,k}^p < 0$, the oscillation of length $2^{J+1-j}$ presents the down phase followed by the up phase, so it is in the second half of the period that the neuron is excited. Of course this phenomenon is stronger when $|W_{j,k}^p|$ is large, as one can see on Figure 3.

<figure><img src="image_5.png" /><figcaption>Fig. 3 Phase lock of spikes on a given rhythm. The LFP filtered at scale j = 2 (see (4)) is represented in black and the spikes of Neuron 1 are represented by blue dots. The red stripes represent the upper phase of the j = 2 oscillations of the LFP. (Top): Independent case, which means $h_{1,2,\kappa_{+}}^{1} = 0$. The other parameters characterizing the equations (5) and (7) are: $\nu^{1} = 0.003, \mu^{1,2} = 0, \sigma^{1,2} = 0.1$. (Bottom left): Positive effect case, in particular $h_{1,2,\kappa_{+}}^{1} = 0.05\mathbb{1}_{1:2^{8-1}}, h_{1,2,\kappa_{-}}^{1} = 0.05\mathbb{1}_{2^{8:2^{9}-1}}$. New spikes appear during the upper phases, caused by the positive phase lock effect. (Bottom right): Negative effect case, corresponding to $h_{1,2,\kappa_{+}}^{1} = -0.05\mathbb{1}_{1:2^{8-1}}, h_{1,2,\kappa_{-}}^{1} = -0.05\mathbb{1}_{2^{8:2^{9}-1}}$. Some of the spikes that in the independent case were present in the upper phase get inhibited. Simulations are performed with the same background noise (fixed random seed)</figcaption></figure>
**Example 3 (Power locking of the spiking activity).** The spiking frequency of a neuron is affected also by the overall power of a LFP oscillations, regardless of the phase. We can call this phenomenon "power modulated spiking" or "power lock". In order to model this interaction we can make use of the mask $\kappa(W) = |\kappa_w^S(W)|$ and, if we consider the Haar basis to locate the single oscillations, the corresponding interaction function can be defined, for a real parameter a, as

$$h_{p,j,\kappa}^m = a \mathbf{1}_{0<t-\ell_{j,k}<2^{J-j+1}}.$$

For a simulation of this phenomenon, see Figure 4.

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Powerful oscillations triggering neuronal spikes. The LFP filtered at $j = 3$ is represented in black and the spike train of Neuron 1 in blue. (Left): independent case, $h_{1,3}^1(t) = 0$. The parameters characterizing the model equations are: $\nu^1 = 0.004$, $\mu^{1,3} = 0$, $\sigma^{1,3} = 5$. (Right): Positive effect of LFP oscillations on spikes, corresponding to $h_{1,3}^1(t) = 0.005 \mathbf{1}_{t-\ell_{j,k}}$. The oscillations $W_{3,k}^P$ increase the spiking probability for all the duration of the oscillation proportionally to their modulus. Simulations are performed with the same background noise (fixed random seed)</figcaption></figure>

Note that (5) can only hold if a condition of the following type is satisfied, to ensure that the sum of the contributions can indeed be interpreted as a probability and thus takes value in $[0, 1]$:

$$ \left\{ \begin{array}{l} 1 \ge \nu^m + \sum_{\substack{m' \in 1:M \\ u < t}} \max(h_{m'}^m(t-u), 0) + \sum_{\substack{p \in 1:P \\ \kappa_{pj}^S \in \kappa_{pj}^S \\ j, k/\ell_{j,k} < t}} \max(h_{p,j,\kappa_{pj}^S}^m(t-\ell_{j,k}), 0) \bar{w}_{p,j} \\ 0 \le \nu^m + \sum_{\substack{m' \in 1:M \\ u < t}} \min(h_{m'}^m(t-u), 0) + \sum_{\substack{p \in 1:P \\ \kappa_{pj}^S \in \kappa_{pj}^S \\ j, k/\ell_{j,k} < t}} \min(h_{p,j,\kappa_{pj}^S}^m(t-\ell_{j,k}), 0) \underline{w}_{p,j} \end{array} \right. \qquad (6) $$

## 2.3 Evolution of the wavelet coefficients

The model of evolution for the time series of the LFP/EEG signals is inspired by Multivariate Autoregressive models for the continuous processes (Faes et al., 2012), except that it is written as an evolutionary process for the wavelet coefficients, keeping
in mind that a given $W_{j,k}^p$ is the strength of the oscillation of scale $j$ starting at the time location $\ell_{j,k}$. This leads to the following equation.

$$W_{j,k}^p = \mu^{p,j} + \sum_{\substack{m \in 1:M \\ u < \ell_{j,k}}} h_m^{p,j}(\ell_{j,k}-u)S_u^m + \sum_{\substack{p', j' \in 1:P \\ \kappa_W^{p',j'} \in \mathbb{R}^W \\ j', k' \in \ell_{j',k'} < \ell_{j,k}}} h_{p',j',\kappa_W^{p',j'}}^{p,j}(\ell_{j,k}-\ell_{j',k'})\kappa_{p',j'}^W(W_{j',k'}^{p'}) + \varepsilon_{j,k}^p \quad (7)$$

Here, all the $\varepsilon_{j,k}^p$ are independent Gaussian variables with mean 0 and variance $(\sigma^{p,j})^2$. The parameter $\mu^{p,j}$ gives the baseline trend of the wavelet coefficient $W_{j,k}^p$. If we use the Haar basis, set $J = 0$ and discard the interactions with spikes, we have that $W_{0,k}^p = cX_t^p$ for some real constant $c$ and with $t = k$. In this case, we recover a classical MVAR model where $X_t^p$ is piecewise constant and depends as a linear regression on the past values of $X_s^{p'}$, for $p' \in 1:P$ and $s < t$.

The Spikes→LFP phenomena, which are interactions of the wavelet coefficients with the point process part in equation (7) can be easily understood in the light of the next example.

**Example 4** (Spikes triggering a given rhythm). If the function $h_m^{p,j}$ is large over $1: 2^{J+1-j}-1$ and if neuron m spikes a lot, the next wavelet coefficient increases. This is what we see in Figure 5.

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Spikes triggering a given rhythm. The LFP filtered at $j = 6$ (4) is represented in black and the spike train of Neuron 1 in blue. The red stripes of length 250 appear after the spikes of Neuron 1 and represent the times intervals in which they can have an effect on the LFP oscillations in case of a nonzero interaction function. (Left): independent case, $h_1^{1,6}(t) = 0$. The parameters characterizing equations (5) and (7) are: $\nu^1 = 0.004$, $\mu^{1,6} = 0$, $\sigma^{1,6} = 0.1$. (Right): Positive effect of spikes on LFP oscillations, corresponding to $h_1^{1,5=6}(t) = 0.61 \mathbb{I}_{t \in [1,250]}$. The oscillations $W_{6,k}^p$ which start inside the red stripes are positively affected by the spikes and their power is enhanced. If more stripes overlap, effects from multiple spikes add up. Simulations are performed with the same background noise (fixed random seed)</figcaption></figure>

Interaction between LFP rhythms can be interpreted in the same way. Again one might want to use masks $\kappa_{p'j'}^W$ such that $\kappa_{p'j'}^W(0) = 0$. These masks are not forced to be bounded, as 1-Lipschitz is a sufficient constraint that is assumed in the sequel. This is because any wavelet coefficient $W_{j,k}^p$ can potentially assume any real value. A typical
mask is the identity $\kappa_{p'j'}^W(W_{j',k'}^{p'}) = W_{j',k'}^{p'}$, but one might want to use a different set of masks to describe a phase-amplitude modulation phenomenon among rhythms, as in the two following examples.

**Example 5** (A fast rhythm power that is in power-phase modulation with slower oscillations). Imagine $j > j'$ and use $\kappa^1(w) = \max(w, 0)$, $\kappa^2(w) = \min(w, 0)$, then with

$$h_{p',j',\kappa^1}^{p,j} = a1_{1:2^{j-j'}-1} \quad \text{and} \quad h_{p',j',\kappa^2}^{p,j} = a1_{2^{j-j'}:2^{j+1-j'}-1},$$

we obtain that

$$h_{p',j',\kappa^1}^{p,j}(\ell_{j,k} - \ell_{j',k'})\kappa^1(W_{j',k'}^{p'}) + h_{p',j',\kappa^2}^{p,j}(\ell_{j,k} - \ell_{j',k'})\kappa^2(W_{j',k'}^{p'}),$$

is large when $|W_{j',k'}^{p'}|$ is large and when $\ell_{j,k}$ is in the upper part of the oscillation corresponding to $W_{j',k'}^{p'}$, whatever the sign of $W_{j',k'}^{p'}$. A fast rhythm modulated by the upper phase of a slower oscillation is given in Figure 6.

<figure><img src="image_8.png" /><figcaption>**Fig. 6** A fast rhythm power that is in phase-power modulation with slower oscillations. The slow rhythm (LFP filtered at $j = 3$, see (4)) is presented in black whereas the fast rhythm (LFP filtered at $j = 6$, see (4)) is in blue. The red stripes correspond to the up phases of the slower oscillation. (Left): Independent case, thus $h_{1,3,\kappa_1}^{1,6} = h_{1,3,\kappa_2}^{1,6} = 0$. The parameters of the model (7) are: $\mu^{1,3} = \mu^{1,6} = 0$, $\sigma^{1,3} = \sigma^{1,6} = 6$. (Right): Positive effect case with $h_{1,3,\kappa_1}^{1,6} = 31_{2^7:2^8-1}h_{1,3,\kappa_2}^{1,6} = 31_{1:2^7-1}$. The fast oscillations get enhanced whenever they correspond to the times of the upper phase of the slow oscillation. Simulations are performed with the same background noise (fixed random seed).</figcaption></figure>

**Example 6** (Power modulation between two LFP rhythms). A common phenomenon in the neuroscientific literature and more broadly, in signal analysis, is the power modulation between two rhythms. Maintaining the notation of our framework, we can model this interaction by defining the mask $\kappa(W) = |W|$. If we employ the Haar wavelet basis, one can use for a fixed positive parameter a,

$$h_{p',j',\kappa}^{p,j} = a1_{0<\ell_{j,k}-\ell_{j',k'}<2^{J-j'+1}}.$$

With this modeling choice in the wavelet equation, any target oscillation $W_{j,k}^p$ is affected by the power of the past source oscillation $W_{j',k'}^{p'}$ which is the closest to the start of the target one. For a simulation of this phenomenon, see in Figure 7.
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Couple of LFP rhythms with a power modulation. The LFP filtered at *j* = 7 is represented in black LFP filtered at *j* = 8 is represented in blue. (Left): independent case, $h_{1,7}^{1,8}(t) = 0$. The parameters characterizing the model equations are: $\mu^{1,7} = 0, \sigma^{1,7} = 0.1, \mu^{1,8} = 0, \sigma^{1,8} = 0.25$. (Right): Positive modulation of the LFP oscillations from frequency *j* = 7 targeting frequency *j* = 8, corresponding to $h_{1,7}^{1,8}(t) = 1.51_{t \in [1,250]}$. The oscillations $W_k^{1,7}$ increase the power of the oscillations $W_k^{1,8}$ which start at the same time or shortly after, in particular the ones such that $0 < \ell_{8,k'} - \ell_{7,k} < 2^{J-7+1} = 16$. The effect is proportional to the modulus $|W_k^{1,7}|$. Simulations are performed with the same background noise (fixed random seed).</figcaption></figure>

**Definition 1 (Connectivity graph).** For a HM-MVAR model given by (5) and (7), the corresponding connectivity graph is defined as a graph where vertices correspond to either neurons (from 1 to M) or couples (p, j) standing for the scale j of the pth LFP. A directed arrow from one vertex i' to another i means that there is a non null interaction function $h_{i'}^i$. When the graph is represented in a figure, the colors of the edges can stand for instance for excitation or inhibition, style for the type of interaction (masks) and width for the strength of the interaction function. Some examples are given in Sections 5 and 6.

**Practical perspective 5 (Functional connectivity).** Functional connectivity graphs are often constructed in neuroscience at different stages of the experiment. Neuroscientists can then discuss whether and how these graphs can serve as markers of certain brain states or cognitive functions (like Slow-Wave sleep - see Section 6). Typically, these functional connectivity graphs are obtained by drawing a (non directed) link between two variables if they are correlated (Faes et al., 2012). Our HM-MVAR graph provides more precise information. For instance, if the HM-MVAR graph contains the arrows 1 $\rightarrow$ 2 $\rightarrow$ 3 between 3 different variables, it means that 1 (Granger) causes 2, which in turns (Granger) causes 3. A simple graph of correlation would only show that 1-2, 2-3 and 3-1 are correlated pairs. In this sense, HM-MVAR is proposing a deeper understanding of the dependencies between the variables. Moreover, beyond the graph, HM-MVAR provides interaction functions that contain valuable information on the time lag between sources and targets and the intensity of the interaction. Note that if the graph is always complete, meaning all edges exist, it would be difficult to find relevant neurophysiological interpretations. This is why we will seek solutions with few nonzero interactions, resulting in sparse graphs. We will discuss this topic in detail in Section 4.
From now on, we are interested in a model where both (5), (6) and (7) hold. The next section focuses on giving conditions for the existence of such processes on $\mathbb{Z}$, whereas Section 4 focus on the estimation of the unknown parameters of the model.

# 3 Existence of a block stationary solution

The following construction is technical, but the main point is to rewrite the HM-MVAR process in a shape where Doukhan and Wintenberger’s result (Doukhan and Wintenberger, 2007) can be applied to show block stationarity of the process. Also the vision by block of the process is useful for the theoretical statistical results of Section 4. For readers with a more biological background, this section might be skipped in a first approximation except Practical perspective 7.

## 3.1 Block notation

Let us vectorize the previous equations. For all time indices $t$, indices $p$ and $j$, we denote $W_{t}^{p,j}$ the process defined by $W_{j,k}^{p}$ when $t = \ell_{j,k}$ and 0 otherwise.

We denote by $Y_t$ the vector of $\mathbb{R}^{M+P(J+1)}$ such that its first $M$ coordinates are the $S_t^m$ and the next are the $W_t^{p,j}$ (using the lexicographic order on $(p, j)$).

We also denote by $\kappa^S(Y_t)$ (resp. $\kappa^W(Y_t)$) the finite vector formed as $Y_t$ but where a given coefficient $W_t^{p,j}$ is replaced by a family of coefficients $\kappa_{pj}^S(W_t^{p,j})$ (resp. $\kappa_{pj}^W(W_t^{p,j}))$, for $\kappa_{pj}^S \in \mathbb{K}_{pj}^S$ (resp. $\kappa_{pj}^W \in \mathbb{K}_{pj}^W$). Note that $\kappa^S(Y_t)$ (resp. $\kappa^W(Y_t)$) is a vector with at most $M+KP(J+1)$, where $K$ is the maximal number of masks (cardinality of $\mathbb{K}_{pj}^S \cup \mathbb{K}_p^W j$). Note that if $K=0$, meaning that no interaction originating from wavelets are considered, then $\kappa^S(Y_t) = \kappa^W(Y_t) = (S_t^m)_{m \in 1:M}$.

The matrix formed of columns $Y_t$ for $a \le t \le c$, is denoted $Y_a^c$. In this notation, $a$ might be $-\infty$. In the same spirit, $\kappa^S(Y_a^c)$ (resp. $\kappa^W(Y_a^c)$) is a matrix which is the juxtaposition the columns defined by $\kappa^S(Y_t)$ (resp. $\kappa^W(Y_t))$, for $a \le t \le c$.

For a given matrix $Z$ that has a potentially infinite number of columns labeled from $-\infty$ to $-1$ and $Q$ rows and a given matrix $H$ with $Q$ columns and $D$ rows of functions with support in $1: +\infty$, we denote $H \ast Z$ the vector of dimension $D$ given by

$$H \ast Z = \sum_{k=1}^{+\infty} H(k)Z_{-k},$$

where the product $H(k)Z_{-k}$ is the classical matrix product between the matrix $H$ evaluated at time $k$ and the column of $Z$ with label $-k$ if it exists - one multiplies by 0 if $Z_{-k}$ does not exist in $Z$.

In particular, let $H^S$ be the matrix formed at row $m$ by the $h_{m'}^m(.)$’s and $h_{p,j,\kappa_{p,j}^S}$’s and let $H^W$ be the matrix formed at row $(p, j)$ (in the lexicographical order) by the
$h_m^{p,j}(.)$'s and $h_{p',j',\kappa_{p',j}^W}$'s. We have that for all indices $m$ and $p, j$,

$$[H^S * \kappa^S(Y_{-\infty}^{t-1})]_m = \sum_{m' \in 1:M \atop u < t} h_{m'}^m(t-u)S_u^{m'} + \sum_{p \in 1:P \atop \kappa_{pj}^S \in \kappa_{pj}^S \atop j, k/\ell_{j,k} < t} h_{p,j,\kappa_{pj}}^m(t-\ell_{j,k})\kappa_{pj}^S(W_{j,k}^p)$$

and

$$[H^W * \kappa^W(Y_{-\infty}^{t-1})]_{p,j} = \sum_{m \in 1:M \atop u < \ell_{j,k}} h_m^{p,j}(\ell_{j,k}-u)S_u^m + \sum_{\substack{p' \in 1:P \atop \kappa_{p'}^W \in \kappa_{p'}^W \atop j', k' / \ell_{j',k'} < \ell_{j,k}}} h_{p',j',\kappa_{p',j'}^W}(\ell_{j,k}-\ell_{j',k'})\kappa_{p',j'}^W(W_{j',k'}^{p'})$$

Note that it is straightforward to prove the following control of the $\ell_1$ norm:

$$\forall H, Z, \quad \|H * Z\|_1 \le \sum_{k=1}^{+\infty} \|H(k)\|_{1,\infty} \|Z_k\|_1, \qquad (8)$$

where for any matrix $A$, $\|A\|_{1,\infty} = \sum_i \max_j |A_{i,j}|$.

Finally let us denote for all $b \in \mathbb{Z}$, $i_b = b2^J$, $j_b = (b+1)2^J - 1$ and $X_b = Y_{ib}^{jb}$.

Let us formalize a trivial statement that is useful to define properly block stationarity.

**Definition 2** (Definition of Block Stationarity). If the sequence of blocks $(X_b)_{b \in \mathbb{Z}}$ is stationary, then for all $s, t \in \mathbb{Z}$ such that $s = t \mod (2^J)$, $Y_s^{s+2^J-1}$ has the same distribution as $Y_t^{t+2^J-1}$. This is what we call block stationarity.

In particular, it means that

$$\left( \left( S_u^m \right)_{\substack{m \in 1:M \\ u \in s : s+2^J-1}} , \left( W_{j,k}^p \right)_{\substack{p \in 1:P \\ j \in 0:J \\ \ell_{j,k} \in s : s+2^J-1}} \right) \quad \text{and} \quad \left( \left( S_u^m \right)_{\substack{m \in 1:M \\ u \in t : t+2^J-1}} , \left( W_{j,k}^p \right)_{\substack{p \in 1:P \\ j \in 0:J \\ \ell_{j,k} \in t : t+2^J-1}} \right)$$

have the same distribution if $s = t \mod (2^J)$.

Therefore, in the sequel, we want to derive under which condition the blocks $X_b$ of length $2^J$ may form a stationary sequence. Note that we cannot hope for full stationarity of the vector $Y_t$, because even the number of non zero coefficients in it depends on the pattern of the $\ell_{j,k}$'s and is therefore $2^J$ periodic. To prove block stationarity, we want to apply Theorem 3.1 of Doukhan and Wintenberger (Doukhan and Wintenberger, 2007) and to do so we need to express $X_b$ as a function of the previous blocks. This is what we do in the following section.

## 3.2 Recursion formulas

Let us introduce the vector of innovation $\xi_t = \begin{bmatrix} \mathbf{U}_t \\ \varepsilon_t \end{bmatrix}$ of size $M + P(J + 1)$ such that the first $M$ coordinates, given by $\mathbf{U}_t$, are i.i.d. Uniform variables on $[0, 1]$ and
the other $P(J+1)$ coordinates, $\epsilon_t$, represent the independent noises on the wavelet coefficients, the coordinate $(p, j)$ (in the lexicographical order) being $\epsilon_t^{p,j}$, a centered Gaussian variable of variance $(\sigma^{p,j})^2$.

At time $s$, it is quite clear that if (5), (6) and (7) are satisfied, then one can represent $Y_s$ as a function $g_s$ of the past $Y_{-\infty}^{s-1}$ and of the innovations:

$$Y_s = g_s(Y_{-\infty}^{s-1}, \xi_s) = \left[ \mathbb{1}_{U_s \le \nu+H^S*\kappa^S(Y_{-\infty}^{s-1})} (\mu + H^W * \kappa^W(Y_{-\infty}^{s-1}) + \epsilon_s) \cdot \delta(s) \right],$$

where the first indicator function has to be understood coordinate by coordinate. The product is the product coordinate by coordinate and $\delta(s)$ is the vector of size $P(J+1)$ which at coordinate $(p, j)$ is $\mathbb{1}_{\exists k, \ell_{j,k}=s}$.

Note that $\delta(s)$ is $2^J$ periodic, so that $g_s(\cdot, \cdot) = g_{s+2^J}(\cdot, \cdot)$. This means that the dependence of a block $X_b$ on its past only depends on its past and the innovation vectors. To find the exact function, let us now focus on the first block $X_0$, with index $b=0$, for simplicity. The second column of the block is depending on $Y_{-\infty}^0 = [Y_{-\infty}^{-1}, g_0(Y_{-\infty}^{-1}, \xi_0)]$ etc. So we can recursively define for a given past $y = Y_{-\infty}^{-1}$ of the block $X_0$

$$\begin{cases} \text{if } s=0, & \mathbb{G}_0(y, \xi_0) = g_0(y, \xi_0) \qquad \in \mathbb{R}^{M+P(J+1)} \\ \text{if } s>0, & \mathbb{G}_s(y, \xi_0^s) = [\mathbb{G}_{s-1}(y, \xi_0^{s-1}), g_s([y, \mathbb{G}_{s-1}(y, \xi_0^{s-1})], \xi_s)] \qquad \in \mathbb{R}^{(M+P(J+1))\times(s+1)} \end{cases} \tag{9}$$

At last, if (5), (6) and (7) hold, we always have that for all $b \in \mathbb{Z}$,

$$X_b = \mathbb{G}_{2^J-1}(X_{-\infty}^{b-1}, \zeta_b), \tag{10}$$

with the following notations:

* $X_{-\infty}^{b-1} = Y_{-\infty}^{jb}$ is the past of block $b$ and consists of all the blocks that have happened before block $X_b$,

* $\zeta_b = \xi_{i_b}^{jb}$ is the matrix of innovations for the block $X_b$.

### 3.3 Block stationarity

Now we are ready to apply Doukhan and Wintenberger's result (Doukhan and Wintenberger, 2007).

**Theorem 6.** Let $K$ be the maximal number of masks in $K_{pj}^S$ and $K_{pj}^W$ and let

$$e_b = (K \lor 1) \left( \max_{k \in i_{b-1}+1:j_b} \|H^S(k)\|_{1,\infty} + \max_{k \in i_{b-1}+1:j_b} \|H^W(k)\|_{1,\infty} \right)$$

and

$$d = (K \lor 1) \left( \max_{k \in 1:2^J-1} \|H^S(k)\|_{1,\infty} + \max_{k \in 1:2^J-1} \|H^W(k)\|_{1,\infty} \right).$$
Assume that (6) holds and that

$$
\left\{
\begin{array}{ll}
\displaystyle \sum_{b=1}^{+\infty} e_b < \frac{d}{(1+d)^{2J}-1} & \text{if } d > 0 \\
\\
\displaystyle \sum_{b=1}^{+\infty} e_b < 2^{-J} & \text{if } d = 0
\end{array}
\right.
$$

then there exists a block τ-weakly dependent stationary solution to (5) and (7), which
is integrable.

If in addition all the masks $\kappa_{pj}^W$'s are bounded, then this block stationary solution
has moment of every order.

Note that the condition on the interaction functions is satisfied as soon they are
summable and small enough. In particular the condition does not depend on $\sigma_{p,j}$.

When $J = 0$, meaning that we do not use the wavelet expansion and end up with
a spike/time series interaction model, we have $d = 0$. In this case, the term $e_b$ is just

$$
e_b = (K \lor 1) (\|H^S(b)\|_{1,\infty} + \|H^W(b)\|_{1,\infty})
$$

In particular, if there is only spikes and no LFP in the system (and therefore
$K = 0$), we recover the condition

$$
\sum_{b=1}^{\infty} \sum_{m \in 1:M} \max_{m' \in 1:M} |h_{m'}^m(b)| < 1.
$$

This condition is therefore slightly more reductive than the classical Hawkes condition
saying that the spectral radius of $(\sum_b |h_{m'}^m(b)|)_{m,m' \in 1:M}$ is strictly smaller than 1
(Brémaud and Massoulié, 1996) (see also (Ost and Reynaud-Bouret, 2018) for more
restrictive conditions).

Therefore our condition is not necessary but only sufficient in this case. We do not
know if it is necessary when considering the whole multiscale approach.

To apply Doukhan and Wintenberger to our problem, we need to rely on a norm
on $\mathbb{R}^{M+P(J+1)}$ that we chose as the $l_1$-norm and an Orlicz norm linked to a convex
function $\Phi$, that we chose equal to identity. This is why we end up proving the existence
of a $L^1$ block stationary solution.

We need to prove two conditions to apply Theorem 3.1 of Doukhan and Wintenberger (Doukhan and Wintenberger, 2007). This is detailed in the two following lemmas. The corresponding proofs are in Section 1 of the Supplementary Material (Spaziani, Girardeau, Reynaud-Bouret, Bethus (2023)).

**Lemma 1.** Let $m > 0$ be a fixed integer. The function $\mathbb{G}_{2^J-1}$ defined in Section 3.2 by (10) and the associated innovation matrix $\zeta_0$ for block 0, satisfy, for all $s \ge 0$

$$
\mathbb{E} (\| \mathbb{G}_{2^J-1}(\mathbf{0}, \zeta_0) \|_1^m) < \infty,
$$
where 0 is the matrix corresponding to null past before 0.

**Lemma 2.** Let $X_{-\infty}^{-1}$ and $Z_{-\infty}^{-1}$ be two different (block) infinite pasts that are possible for the block $X_0$, then

$$ \mathbb{E}(\|G_{2^J-1}(X_{-\infty}^{-1}, \zeta_0) - G_{2^J-1}(Z_{-\infty}^{-1}, \zeta_0)\|_1) \le \sum_{b=1}^{\infty} a_b \mathbb{E}(\|X_{-b} - Z_{-b}\|_1), $$

where

$$ a_b = \begin{cases} ((1+d)^{2^J} - 1) \frac{e_b}{d} & \text{if } d > 0 \\ 2^J e_b & \text{if } d = 0 \end{cases} $$

**Practical perspective 7 (Stationarity interpretation).** Theorem 6 shows that if the interaction are small enough, the HM-MVAR process may be maintained in a stationary distribution indefinitely, with patterns that are shaped by blocks of wavelets. Indeed the pattern of the wavelet basis supports enforces a block structure of the data, that allows the coexistence of variables of different scales. In the case of the Haar wavelet, the length of a wavelet block corresponds to the support of the father wavelet (see Figure 1).

From a more biological point of view, stationarity does not exist in the brain and states are transient even if their length (like the sleep stages, see Section 6) can last up to minutes or hours. However, for the HM-MVAR model, the interaction lags usually range from few milliseconds up to seconds. It is therefore possible that ranges of parameters of HM-MVAR would make this model divergent, by either blow-up or extinction in a few minutes, phenomena which would not be realistic for the modeling purpose we have in mind. Therefore it is of utmost importance to ensure that our model possesses ranges of parameters which do not allow for these phenomena. Theorem 6 provides ranges of strength for the global set of interaction functions h, that ensures block stationarity and therefore prevents explosions or extinctions.

# 4 Sparse estimation by Lasso criterion

In this section, we want to estimate $\mu, \nu$ and the interaction functions h from Spike-LFP data, to reconstruct a sparse connectivity graph between neurons and LFP rhythms. Because formula (5) and (7) are linear in these unknown functions and parameters, and since these interaction functions only depend on a certain lag with respect to the data at hand ($S_t^m$ or $W_{j,k}^p$), we approximate formulas (5) and (7) by a linear combination of cylindrical functions of the past, i.e. functions only depending on a small past neighborhood of the data. Let us present in more details this dictionary of cylindrical functions.

Note that Section 4 proves theoretical results that give guarantees on the whole estimation procedure. Readers with biological background might want to go directly Practical perspective 9.
## 4.1 Dictionary of cylindrical functions

Let us denote

$$ \mathcal{I} = 1 : M \cup (1 : P \times 0 : J), $$

so that an index $i \in \mathcal{I}$ can be either referring to neuron $m$ or to a couple $(p, j)$, standing for the scale $j$ of the $p$th LFP.

**Definition 3 (Neighborhood).** A neighborhood $v$ is a finite subset of $\mathcal{I} \times 1 : +\infty$. It represents a collection of indices of data and time lags.

Using the block notation of the previous Section, we consider real valued functions $\phi$ of the past, that is $\phi : Y_{-\infty}^{t-1} \to \mathbb{R}$, that only depends on the past configuration and not on the precise value of $t$. Note however that because of the pattern of zeros in $Y_{t-2,j}$, which is the same as the one at time $t$, this function $\phi$ has access to the information telling if $t = \ell_{j,k}$ for some $(j, k)$ or not.

**Definition 4 (Cylindrical function).** Such a function $\phi$ of the past is cylindrical on the neighborhood $v$ if $\phi(Y_{-\infty}^{t-1})$ only depends on any $i$th coefficient of $Y_{t-s}$, if $(i, s) \in v$. Informally, it means that $\phi(Y_{-\infty}^{t-1})$ only depends on the data that are in the neighborhood $v$ of $Y_t$. Its maximal lag is given by $\max\{d|\exists i \in \mathcal{I}, (i, d) \in v\}$.

Let us give some illustrative examples of such functions with the Haar basis or any basis that would have an up phase in $1 : 2^{J-j}$ and a down phase in $2^{J-j+1} + 1 : 2^{J-j+1}$ (see also the Supplementary Material). Let us fix some $t \in \mathbb{Z}$. Since we are approximating (5) and/or (7) by a linear combination of these functions, we also give an interpretation of the multiplicative coefficient in the expansion.

1. $\phi(\cdot) = 1$ is cylindrical on the empty neighborhood $v = \emptyset$. Its maximal lag is therefore 0. The coefficient in front of this $\phi$ is typically a candidate to be a $\nu^m$ (if we approximate (5)) or $\mu^{p,j}$ (if we approximate (7)).

2. $\phi(Y_{-\infty}^{t-1}) = S_{t-d}^m$, for some index $m \in 1 : M$ and lag $d > 0$. This $\phi$ is cylindrical on $v = \{(m, d)\}$. The coefficient in front of this $\phi$ is typically a candidate to be a $h_m^{m'}(d)$ (if we approximate (5)) or $h_m^{p,j}(d)$ (if we approximate (7)).

3. For $j > 0$, one can also represent the last coefficient of the $j$th rhythm on the $p$th LFP before time $t$ by

$$ \phi(Y_{-\infty}^{t-1}) = \sum_{k \in \mathbb{Z}} W_{j,k}^p 1_{t-2^{J+1-j} \le \ell_{j,k} < t} = \sum_{s=1}^{2^{J+1-j}} W_{t-s}^{p,j}, $$

since in both sums, only one term is not null, the one equal to $W_{j,k}^p$ with $\ell_{j,k} = s$. This $\phi$ is cylindrical on $v = \{(p, j)\} \times 1 : 2^{J+1-j}$. The coefficient in front of this $\phi$ is typically a candidate to be the constant value of $h_{p,j}^{m'}(d)$ for $d \in 1 : 2^{J+1-j}$ (if we approximate (5)) or $h_{p,j}^{p',j'}(d)$ (if we approximate (7)). This coefficient corresponds to the dependency w.r.t. to the previous oscillation.
4. One can also want to look at the first half of the oscillation, which is especially useful for phase lock phenomena:

$$ \phi(Y_{-\infty}^{t-1}) = \sum_{s=1}^{2^{J-j}} W_{t-s}^{p,j}, $$

is the last non null $W_{j,k}^p$ if its corresponding oscillation starts less than half a period away and 0 otherwise. It is cylindrical on $v = \{(p, j)\} \times 1 : 2^{J-j}$.

5. One can also consider the second half period with

$$ \phi(Y_{-\infty}^{t-1}) = \sum_{s=2^{J-j+1}}^{2^{J+1}-j} W_{t-s}^{p,j}, $$

which is the last non null $W_{j,k}^p$ if its corresponding oscillation starts between half a period and one complete period away and 0 otherwise.

6. In Cases 3-5, one can replace $W_{t-s}^{p,j}$ by masked version, for instance $\kappa_w(W_{t-s}^{p,j})$, $\kappa^1(W_{t-s}^{p,j})$ or $\kappa^2(W_{t-s}^{p,j})$, depending if the function $\phi$ is used to approximate the interactions in (5) or (7).

Also the oracle inequality proved in the sequel only works for bounded dictionaries, which means that we need to use bounded masks also in the wavelet equation. With our previous definition of the mask $\kappa_w$, one can also take $|\kappa_w(W_{t-s}^{p,j})|$ to have access to (an approximation of) the power of the rhythm itself.

7. To avoid having too many coefficients to estimate, one can also want to group several times or cycles in one single $\phi$. For instance, with spikes, one might take, for fixed positive lags $L < R$

$$ \phi(Y_{-\infty}^{t-1}) = \sum_{s=L}^{R} S_{t-s}^{m}, $$

which is in fact the number of spikes in the window $t-R : t-L$. This $\phi$ is cylindrical on the neighborhood $v = \{m\} \times L : R$. Its corresponding coefficient $a_\phi$ would correspond to an approximation of $h_m^{m'}(d)1_{L \le d \le R}$ in (5) or of $h_m^{p,j}(d)1_{L \le d \le R}$ in (7) by a constant value $a_\phi 1_{L \le d \le R}$.

8. The same thing can be done on rhythms, specifically to merge the contribution of a certain number of past oscillations. To do so, we can consider as time lag a multiple of the oscillation duration. So let $L > R$ be fixed positive integers and define for instance

$$ \phi(Y_{-\infty}^{t-1}) = \sum_{s=R2^{J-j+1}+1}^{L2^{J-j+1}} W_{t-s}^{p,j}, $$

which is the sum of all the $L-R$ coefficients of the $j$th rhythm of the $p$th LFP before time $t - R2^{J-j+1}$. This is cylindrical on $\{(p, j)\} \times R2^{J-j+1} + 1 : L2^{J-j+1}$. This might be especially useful if we think that large oscillations in a fast rhythm
for a certain time might trigger a large oscillation at a slower rhythm a few cycles later. This could make also sense if we replace $W_{t-s}^{p,j}$ by $|W_{t-s}^{p,j}|$, to have access to the total past power and avoid cancellations between wavelet coefficients with different signs.

9. Finally, and especially for influence of spikes on LFP, one might want to aggregate the activity of the whole neuronal population and estimate its effect on the oscillations, so that one could consider

$$ \phi(Y_{-\infty}^{t-1}) = \sum_{m=1}^{M} \sum_{s=L}^{R} S_{t-s}^{m}, $$

which is the total number of spikes produced on all recorded, in the window $t-R: t-L$. The corresponding coefficient measures the excitation effect of a population of neurons on a particular rhythm. In the real data application (see Section 6), we will in particular aggregate over all the neurons in a given brain region.

With the previous examples, formula (5) and (7) can be reinterpreted as fixed linear combination of elementary cylindrical functions (Cases 1 to 6 typically). Moreover, if we want to approximate the interaction functions by piecewise constant functions over a coarser grid than the one given by $t$, we might want to use the grouped version (Case 7 and 8 for instance). One might also want to group in particular neurons (Case 9). These dimension reductions (at least the ones on time) are common in the Hawkes setting (Ost and Reynaud-Bouret, 2018), much less in the MVAR one.

In the sequel, we fix the index $i$ we are looking at (either a neuron $m$ or a scale $j$ of a LFP $p$), and fix a finite dictionary $\Phi^i$ of cylindrical functions to explain the conditional expectation $f_t^i(Y_{-\infty}^{t-1}) = \mathbb{E}(Y_{t,i}|\mathcal{F}_{t-1})$, where $Y_{t,i}$ is the $i$th index of $Y_t$.

Note that $f_t^i$ is, similarly to $g_t$ in Section 3, a periodic function of the past $Y_{-\infty}^{t-1}$ and that in the most general case $f_{t+2}^i = f_t^i$. In particular, $f_t^i$ is null when $\delta^i(t)$ is null, where $\delta^i(t)$ is defined by

$$ \delta^i(t) = \begin{cases} 1 & \text{if } i = m \in 1: M \\ 1_{\exists \ell_{j,k}, \ell_{j,k}=t} & \text{if } i = (p, j) \in 1: P \times 0: J. \end{cases} $$

The approximation of $f^i = (f_t^i)_{t \in 0:T}$ on the dictionary $\Phi^i$ is therefore given, for a $a \in \mathbb{R}^{\Phi^i}$ by $f_a^i = (f_{a,t}^i)_{t \in 0:T}$, with

$$ f_{a,t}^{i}(Y_{-\infty}^{t-1}) = \sum_{\phi \in \Phi^{i}} a_{\phi} \phi(Y_{-\infty}^{t-1}) \delta^{i}(t). $$

We work with the following random norm on the vectors of functions of the past, $q = (q_t)_{t \in 0:T}$

$$ \|q\|^2 = \sum_{t \in 0:T} [q_t(Y_{-\infty}^{t-1})]^2, $$

and associate scalar product $\langle q, q' \rangle = \sum_{t \in 0:T} q_t(Y_{-\infty}^{t-1})q'_t(Y_{-\infty}^{t-1})$. For intuition, in the sequel, note that if $q$ is $2^J$ periodic, which is the case for $f^i$ and $f_a^i$, these quantities
should be growing with $N$ the number of blocks. If the sequence $Y$ is block stationary, we can even hope by ergodicity that this grows linearly with $N$.

On the other hand, we also use the classical euclidian norm on the vector $a \in \mathbb{R}^{\Phi^i}$ given by

$$ \|a\|^2 = \sum_{\phi \in \Phi^i} a_{\phi}^2 $$

and associated scalar product also denoted $\langle \cdot, \cdot \rangle$. We use the same notation for both norms but note that depending on the objects it has different meaning.

## 4.2 Contrast and Lasso penalization

Let $A$ be the maximal lag of all the cylindrical functions $\phi$ used in dictionary $\Phi$. Let $B$ be the largest integer such that $B2^J \ge A$. Let us also set $N$ and $T$ two positive integers such that $T = N2^J - 1$.

Up to a shift in time of the data, let us assume that we observe the $P$ time series and the $M$ spike trains on $-(B+1)2^J : (N+1)2^J - 1$. Thanks to the wavelet transformation (Gubner and Chang, 1995), it means we that we have access to the $W_{j,k}^p$ for all $\ell_{j,k} \in -A : T$, which means that all $f_a^i(Y_{-\infty}^{t-1})$ are computable for all $i \in \mathcal{I}$ and $t \in 0 : T$.

We consider the following least-square contrast.

**Definition 5 (Least-square contrast).** The least-square contrast is defined by:

$$ C(f_a^i) = -2 \sum_{t=0}^{T} f_{a,t}^i(Y_{-\infty}^{t-1})Y_{t,i} + \sum_{t=0}^{T} [f_a^i(Y_{-\infty}^{t-1})]^2 = -2 \langle f_a^i, Y^i \rangle + \|f_a^i\|^2, $$

with $Y^i = (Y_{t,i})_{t \in 0:T}$. This can be rewritten as $C(f_a^i) = -2a^\top b_i + a^\top G_i a$, where $a = (a_\phi)_{\phi \in \Phi^i}$ is the vector of unknown coefficients,

$$ b_i = \left( \sum_{t=0}^{T} \phi(Y_{-\infty}^{t-1})Y_{t,i} \right)_{\phi \in \Phi^i} \quad \text{and} \quad G_i = \left( \sum_{t=0}^{T} \phi(Y_{-\infty}^{t-1})\phi'(Y_{-\infty}^{t-1}) \right)_{\phi, \phi' \in \Phi^i}. $$

We want sparse solutions, which is why we want to solve this weighted Lasso problem (Hansen et al., 2015):

$$ \hat{a}_i = \arg \min_{a \in \mathbb{R}^{\Phi^i}} -2a^\top b_i + a^\top G_i a + \gamma w_i^\top |a|, \quad (11) $$

where $\top$ denotes the transpose, $w_i$ is a weight vector, which adapts the penalty to each $\phi$ (Weighted Lasso), $\gamma > 0$ is a tuning parameter and $|a| = (|a_\phi|)_{\phi \in \Phi^i}$. This will lead us to a sparse reconstructed connectivity graph.

**Definition 6 (Reconstructed connectivity graph).** This connectivity graph has vertices given by $\mathcal{I}$. A directed edge from $i' \in \mathcal{I}$ to $i \in \mathcal{I}$ means that

* there is a non zero $\hat{a}_{i,\phi}$ in the Lasso problem of $i$,
• for a $\phi$ cylindrical on a neighborhood $v$ such that there exists some positive lag $d$ for which $(i', d) \in v$.

The edges of the graph are coded to represent which dictionary element $\phi$ is used so that one can have multiple arrows of different type from on node to another. The delays of interaction lie also in the neighborhood $v$ on which the selected $\phi$ is cylindrical.

To evaluate the full connectivity graph between spike trains and rhythms, we solve successively all the Lasso problems for each $i \in \mathcal{I}$. For each of these Lasso problems, we might want to tune the dictionary $\Phi^i$ and the weights $w_i$: in particular they depend on the nature of $i$ (neurons or rhythm).

## 4.3 Choice of minimal weights $d$

In this section, we derive a set of minimal weights for which our method has theoretical guarantees of achieving good performance. All the $w_{i,\phi}$ used in the sequel will be such that $w_{i,\phi} \ge d_{i,\phi}$, where the values of $d_{i,\phi}$ are derived in this section. This choice of minimal weights $d_i$ is based on exponential deviation inequalities, derived from the martingale structure of our model. The first result holds for the spike data and is inspired by similar results in continuous time (Hansen et al., 2015).

**Lemma 3.** For any bounded function $\phi: Y_{-\infty}^{t-1} \to \mathbb{R}$, and any index $m$, let

$$M_T^\phi = \sum_{t=0}^{T} \phi(Y_{-\infty}^{t-1}) [S_t^m - \mathbb{E}(S_t^m | \mathcal{F}_{t-1})] \quad \text{and} \quad \hat{V}_T^\phi = \sum_{t=0}^{T} [\phi(Y_{-\infty}^{t-1})]^2 S_t^m.$$

Then for any $\epsilon, x > 0$,

$$\mathbb{P} \left( M_T^\phi \ge \sqrt{3(1+\epsilon)\hat{V}_T^\phi} x + \left(\frac{1}{3} + \sqrt{6(1+\epsilon)}\right) \|\phi\|_\infty x \right) \le 2 \left( \frac{\log\left(\frac{T+1}{2x} + 1\right)}{\log(1+\epsilon)} + 1 \right) e^{-x}.$$

As a consequence, for any family $\Phi$ of such functions $\phi$, with finite cardinality $\#\Phi$, if $T > 1$ and $\alpha \in (0, 1)$, the event

$$\Omega_{\Phi}^{\alpha} = \{\forall \phi \in \Phi, |M_{T}^{\phi}| \le d_{\phi}\},$$

with

$$d_{\phi} = \sqrt{3(1+\epsilon)\hat{V}_{T}^{\phi} \log\left(\frac{4\#\Phi\left(\frac{\log(T)}{\log(1+\epsilon)} + 1\right)}{\alpha}\right)} + \left[\frac{1}{3} + \sqrt{6(1+\epsilon)}\right] \log\left(\frac{4\#\Phi\left(\frac{\log(T)}{\log(1+\epsilon)} + 1\right)}{\alpha}\right) \|\phi\|_{\infty},$$

is of probability at least $1 - \alpha$. In particular, one can take

$$d_{\phi} = \sqrt{\frac{9}{2} \hat{V}_{T}^{\phi} \log \left( \frac{16 \# \Phi \log(T)}{\alpha} \right)} + \frac{10}{3} \log \left( \frac{16 \# \Phi \log(T)}{\alpha} \right) \|\phi\|_{\infty}, \quad (12)$$
**Lemma 4.** For any bounded function $\phi: Y_{-\infty}^{t-1} \to \mathbb{R}$, and any index $p, j$, let

$$M_T^\phi = \sum_{k/0 \le \ell_{j,k} \le T} \phi(Y_{-\infty}^{t-1}) [W_{j,k}^p - \mathbb{E}(W_{j,k}^p | \mathcal{F}_{\ell_{j,k}-1})] \quad \text{and} \quad V_T^\phi = \sigma_{p,j}^2 \sum_{k/0 \le \ell_{j,k} \le T} [\phi(Y_{-\infty}^{t-1})]^2.$$

Then for any $x, \epsilon, \eta > 0$,

$$\mathbb{P} \left( M_T^\phi \ge \sqrt{2(1+\epsilon) \left[ V_T^\phi + \eta \sigma_{p,j}^2 \|\phi\|_\infty^2 \right]} x \right) \le \left( \frac{\log\left(\frac{T+1}{\eta}\right) + 1}{\log(1+\epsilon)} + 1 \right) e^{-x}.$$

As a consequence, for any family $\Phi$ of such functions $\phi$, with finite cardinal $\#\Phi$, if $T > 1$ and $\alpha \in (0, 1)$, the event

$$\Omega_{\Phi}^{\alpha} = \{\forall \phi \in \Phi, |M_{T}^{\phi}| \le d_{\phi}\},$$

with

$$d_{\phi} = \sqrt{2(1+\epsilon)[V_{T}^{\phi} + \eta\sigma_{p,j}^{2}\|\phi\|_{\infty}^{2}]\log\left(\frac{\#\Phi\left(\frac{\log((T+1)\eta^{-1}+1)}{\log(1+\epsilon)}\right) + 1}{\alpha}\right)}$$

is of probability at least $1 - \alpha$. In particular one can take

$$d_{\phi} = \sqrt{3[V_{T}^{\phi} + 0.5\sigma_{p,j}^{2}\|\phi\|_{\infty}^{2}]\log\left(\frac{36\#\Phi\log(T)}{\alpha}\right)} \quad (13)$$

Note that to use this $d_{\phi}$, we need to know $(\sigma_{p,j})^2$, but an upper bound can work as well and we can take in practical applications an overall estimation of the variance of $W_{j,k}^p$.

## 4.4 Oracle inequality

Before proving oracle inequalities, let us state the property we assume on the matrices $G_i$'s. In the following, for a subset $J \subset \Phi^i$ and $x \in \mathbb{R}^{\Phi^i}$, we denote $x_J = (x_\phi)_{\phi \in J}$, $|x_J|_1 = \sum_{\phi \in J} |x_\phi|$, $\|x_J\|^2 = \sum_{\phi \in J} |x_\phi|^2$ and $J^c = \Phi^i \setminus J$. We also note $S(x)$ the support of the vector $x$ that is $S(x) = \{\phi \in \Phi^i / a_\phi \neq 0\}$.

**Definition 7 (Property RE(κ, c, s) for Gᵢ/N).** For $κ, c$ real positive numbers and $s$ positive integer, $G_i/N$ has the Restricted Eigenvalue property RE(κ, c, s), if for all $x \in \mathbb{R}^{\Phi^i}$ and $J$ of cardinal less than $s$, such that $|x_{J^c}|_1 \le c|x_J|_1$, we have

$$x^\top G_i x \ge N \kappa \|x_J\|^2.$$

First of all, we have a factor $N$ the number of blocks because in the stationary regime, we can hope that by ergodicity, $G_i/N$ converges. This has been proved for
instance in the pure Hawkes case (Ost and Reynaud-Bouret, 2018). Moreover, by controlling both the distance to this limit and the limit itself, it has been proved (see (Ost and Reynaud-Bouret, 2018)) that if the models have some fast decaying Kalikow decomposition, with large probability, the $G_i/N$ matrix satisfies $\mathbf{RE}(\kappa, c, s)$ for some parameters $\kappa, c$ and $s$. The same argument is not possible here because Kalikow decomposition do not exist for Gaussian processes (up to our knowledge). On the other hand, there are many works on Gaussian matrices to show under which property $\mathbf{RE}(\kappa, c, s)$ is satisfied (see (Raskutti et al., 2010) and the references therein). We do not know how to prove such properties when the process is a mix of both Hawkes and Gaussian cases. This is why the following theorem is stated on the event where $\mathbf{RE}(\kappa, c, s)$ holds.

**Theorem 8.** Let $i \in \mathcal{I}$ be a fixed index, with associated finite dictionary of cylindrical functions $\Phi^i$ and let $\gamma$, the tuning parameter, be such that $\gamma > 2$.

Let us fix some weights $w_i$ such that $w_{i,\phi} \ge d_{i,\phi}$, for all $\phi \in \Phi^i$ where the $d_i$'s and the corresponding $\Omega_{\Phi^i}^\alpha$ are given in Lemma 3 or 4, depending on the type of $i$ (spike train or rhythm). Let also $\gamma > 2$.

On the event

* where $\Omega_{\Phi^i}^\alpha$ of probability larger than $1 - \alpha$ holds

* and where $G_i/N$ satisfies $\mathbf{RE}(\kappa, c, s)$ with

$$c = \left( \frac{\gamma + 2}{\gamma - 2} \right) \frac{\max_{\phi \in \Phi^i} w_{i,\phi}}{\min_{\phi \in \Phi^i} w_{i,\phi}},$$

one has that $\hat{a}_i$ defined by (11) satisfies

$$\|f_{\hat{a}_i}^i - f^i\|^2 \le \inf_{a/\#S(a) \le s} \left\{ \|f_a^i - f^i\|^2 + \frac{(\gamma+2)^2}{4\kappa} \sum_{\phi \in S(a)} \frac{w_{i,\phi}^2}{N} \right\}.$$

Moreover if there exists $a_i^*$ such that $f^i = f_{a_i^*}^i$ and $\#S(a_i^*) \le s$, then

$$|\hat{a}_i - a_i^*|_1 \le \frac{2\gamma(\gamma+2)}{2\kappa(\gamma-2)\sqrt{N}} \sum_{\phi \in S(a^*)} \frac{w_{i,\phi}^2}{\sqrt{N} \min_{\phi \in \Phi^i} w_{i,\phi}}.$$

Note first that if we take $w_{i,\phi} = d_{i,\phi}$, the leading term in $w_{i,\phi}^2$ is of order $\hat{V}_T^\phi$ in Lemma 3 and $V_T^\phi$ in Lemma 4, times a logarithmic factor of the form $\log(\#\Phi^i) + \log\log(T)$. In particular, it increases with the number of blocks $N$ and we can hope by ergocidicity that $\hat{V}_T^\phi/N$ and $V_T^\phi/N$ converges to a constant. Also in both cases these terms are estimating a variance. Therefore one can see roughly $w_{i,\phi}^2/N$ as a constant with respect to $N$ up to logarithmic term. Hence the first equation is an oracle inequality stated on a norm to the square that should increase linearly in $N$ and for which the bias term is counterbalanced by a variance term that should remain roughly constant and proportional to the number of non zero coefficients in $a$.
The convergence when $N$ increases is more clear in the second equation. If there is indeed a sparse representation of $f^i$ on the dictionary, then one is able by the previous procedure to retrieve this set of coefficients and the error on the coefficients decrease like $[\log(\#\Phi^i) + \log \log(T)]^{1/2}/\sqrt{N}$, which is what one would expect for such convergence results in Lasso setting, except the factor $\log(\log T)$, which comes from a very sharp control of the martingales in our dependent setting.

Note however that we do not need for the previous theorem to assume that the sequence of observations $Y$ is block stationary to prove the theorem. The stationary case just gives us the intuition of the order of magnitude of the several terms.

The constant $c$ in $\mathbf{RE}(\kappa, c, s)$ depends both on $\gamma$ and the weights. In particular if $\gamma$ is large or close to 2 or if the weights have a very large amplitude, $c$ becomes very large and $\mathbf{RE}(\kappa, c, s)$ is much less likely to be fulfilled. One could in particular have a smaller $c$ by taking constant weights $w_{i,\phi} = \max_{\phi \in \Phi^i} d_{i,\phi}$, so that $c = (\gamma+2)/(\gamma-2)$.

However, in practice the choice of varying weights makes the estimation much more robust to heterogeneity of the inputs and of the dictionary. This kind of comparison on simulation has already been made in (Hansen et al., 2015; van de Geer et al., 2011).

**Practical perspective 9** (Interpretation of the oracle inequality and the role of the dictionary). Informally, the first equation of Theorem 8 guarantees that, even if the data dynamics are not created strictly by the model equations ((5),(7)), we can retrieve the best functional connectivity graph within the model that fits such dynamics. Indeed, we have, up to multiplicative constants, the best compromise between the bias of the model and a price to pay for each coefficient that is non zero. This price is given by $w_{i,\phi}^2/N$, which is, in the case, of the minimal choices of Lemma 3 or 4 is exactly the size of the fluctuation of each individual estimator and is therefore unavoidable. The second equation tells us that in an ideal world where indeed the dynamics were generated by the model, then the longer the recording the better becomes our precision in retrieving the true underlying connectivity.

What is important for the biological interpretation of our estimator, is that our method is only as good as the dictionary that is filled in entry. This dictionary needs to be broad enough so that it can catch most of the complex patterns that may exist. But on the other hand, if we filled this with dictionary elements that are not interpretable per se (like masks that would not be monotonous) or if we fill it with too few elements (for instance only phase lock phenomenon and no power phenomenon), the resulting graph might be poor. Trying to improve the dictionary is a task for neurobiologists, who can iteratively add meaningful new interactions thanks to the framework presented in this work. In practice, it is impossible to create a perfect dictionary which catches all the complex brain dynamics. Also, the weights $w_{i,\phi}$ depends on the global size of the dictionary and searching in a gigantic dictionary also deteriorates the performance (typical convergence in $\log(\#\Phi^i)/\sqrt{N}$). Still we can have the guarantee, thanks to Theorem 8, that whatever the dictionary, the method makes the best possible approximation of the data by the possible pattern of interactions that are in the dictionary.
# 5 Simulations

In this section we test the performance of the previous LASSO estimation procedure on simulated data.

## 5.1 Overall simulation plan

We are only using here the Haar basis both for simulation and reconstruction. To make biologically relevant numerical experiments, we do as if we have a sampling frequency $f_S$ and all the time are measured in seconds. All the following experiments follow the same schedule:

1. We performed 5 different numerical settings, where each time a different sparse true connectivity graph is used. Each setting corresponds therefore to the choice of the number of neurons $M$ and LFP signals, $P$, as well as the finer scale $J$. Then each neuron $m$ and LFP rhythm $(p, j)$ is characterized by its spontaneous value $\nu^m$ $(\mu^{p,j})$ and by its interactions with the other variables $h_m^j(t)$ $(h_{j,m}^j(t))$. To simplify, we chose each time to be in the case where we can decompose our process on a specific dictionary, presented in the next section. The dictionary is chosen to mimic patterns of cerebral activity that are known in the literature (see also the examples of Section 2). The parameters definitions are given in Table 2. The values of the parameters which remain constant throughout all simulations are also written in the table, they are of the order of magnitude of classical experiments in neuroscience. (Note: when we set the $\nu^m$ in the tables of parameters, we are writing the biology-friendly $\nu^m \cdot f_S$ - in Hz). The remaining values are given in Figures 9 and 10 for each setting.

<table><caption>**Table 2** Definition of the simulation parameters.</caption><thead><tr><th>Parameter</th><th>Definition</th></tr></thead><tbody><tr><td>Simu# = 100</td><td>number of simulations for a given set of parameters.</td></tr><tr><td>f<sub>S</sub> = 1250</td><td>sampling frequency for the spikes and the LFPs (in Hz).</td></tr><tr><td>M</td><td>number of neurons.</td></tr><tr><td>P = 1</td><td>number of LFP signals.</td></tr><tr><td>J = 10</td><td>finest scale.</td></tr><tr><td>N = 70</td><td>number of blocks of 2<sup>J</sup> LFP points (see Section 4.2).<br>The real duration of each block is 2<sup>J</sup>/f<sub>LFP</sub> seconds.</td></tr><tr><td>Υ = 57.344</td><td>observation duration, equal to 2<sup>J</sup>/f<sub>LFP</sub> N in seconds.</td></tr><tr><td>γ<sub>S</sub></td><td>LASSO penalization parameter when estimating the parameters for a spiking probability.</td></tr><tr><td>γ<sub>W</sub></td><td>LASSO penalization parameter when estimating the parameters for a wavelet coefficient.</td></tr><tr><td>ν<sup>m</sup></td><td>spontaneous firing rate for neuron m (in Hz).</td></tr><tr><td>μ<sup>p,j</sup> = 0</td><td>spontaneous value for LFP p, rhythm j (in mV).</td></tr><tr><td>σ<sup>p,j</sup> = 0.25</td><td>standard deviation of the Gaussian noise ϵ<sub>j,k</sub><sup>p</sup> (see 7).</td></tr><tr><td>γ<sub>S</sub> = 0.6</td><td>penalization parameter (see 11) when estimating the parameters of 5.</td></tr><tr><td>γ<sub>W</sub> = 0.3</td><td>penalization parameter (see 11) when estimating the parameters of 7.</td></tr></tbody></table>

2. Thanks to equations (5) and (7) we can simulate the evolution of the spiking activity of the neurons and of the wavelet coefficients of the LFPs.
3. The synthetic data is then employed for the estimation of the original connectivity graph. The choice of a dictionary $\Phi$ permits the computation of the vectors $b_i$, the Gram matrices $G_i$, and of the weights $w_i$. The latter are chosen for the neuron part as in (12) and for the wavelet part as in (13), with the slight difference that $(\sigma^{p,j})^2$ is over-estimated by the empirical variance of each of the wavelet coefficients.

4. Finally we studied 3 different methods: Least Squares without penalization, denoted LS, Lasso (11), and Lasso followed by a Least Square step on the nonzero parameters found by the Lasso, denoted L+LS.

5. The estimation performance is evaluated with the help of various measures. The differences between the two connectivity graphs can be analyzed by counting missing or extra edges, and in general how many times the graph is perfectly reconstructed. Then, the precision of the estimation is assessed by measuring the errors on the spontaneous values and on the coefficients. The definition of all the performance measures we used is given in Table 3 and more precise mathematical definitions of the errors are given in Section 2 of the Supplementary Materials. Note that in the sequel, we give the average values of all the performance measures over the 100 simulations that are made.

<table><thead><tr><th>Names</th><th>Meaning</th></tr></thead><tbody><tr><td>$\nu_e$</td><td>total error for the spontaneous firing rates</td></tr><tr><td>$\mu_e$</td><td>total error for the LFP spontaneous values</td></tr><tr><td>$e_{sS}$</td><td>total error for spike-spike interactions</td></tr><tr><td>$e_{sW}$</td><td>total error for spike-LFP interactions</td></tr><tr><td>$e_{wS}$</td><td>total error for LFP-spikes interactions</td></tr><tr><td>$e_{wW}$</td><td>total error for LFP-LFP interactions</td></tr><tr><td>DG</td><td>=1 if graph is perfectly reconstructed</td></tr><tr><td>S</td><td># of spontaneous values detected</td></tr></tbody></table>

<table><caption>**Table 3** Definition of the various measures of performance. More precise definitions are available in the Supplementary Material.</caption><thead><tr><th>Names</th><th>Meaning</th></tr></thead><tbody><tr><td>$c+sS$</td><td># false positive spike-spike parameters</td></tr><tr><td>$c-sS$</td><td># false negative spike-spike parameters</td></tr><tr><td>$c+sW$</td><td># false positive spike-LFP parameters</td></tr><tr><td>$c-sW$</td><td># false negative spike-LFP parameters</td></tr><tr><td>$c+wS$</td><td># false positive LFP-spike parameters</td></tr><tr><td>$c-wS$</td><td># false negative LFP-spike parameters</td></tr><tr><td>$c+wW$</td><td># false positive LFP-LFP parameters</td></tr><tr><td>$c-wW$</td><td># false negative LFP-LFP parameters</td></tr></tbody></table>

## 5.2 Precise tuning of the dictionary for neuroscience purpose

Throughout all the simulations we used a dictionary that is specially tuned to detect interaction that are well-known in the biological literature. In particular, we focused on the interactions that can appear among and between the hippocampus and the amygdala in a mouse during a specific sleep phase called "slow wave sleep" (see Section 6 on experimental data). The scales $j$ are associated with classical brain rhythms (see Table 1) to design a dictionary that is meaningful from this neuroscientific point of view.

### 5.2.1 Dictionary for a neuron

Let us first describe the dictionary to reconstruct the interactions targeting a neuron:
1. The Spontaneous firing rate $\nu^m$ corresponds to the coefficient of

$$ \phi_1(Y_{-\infty}^{t-1}) = 1 \qquad (14) $$

2. Effect of other neurons. This function is needed to detect the increase or decrease of the spiking probability due to the firing of other neurons. Given a time $t \in T$, and lag $r \in \mathbb{N}$ and a bin size $\delta \in \mathbb{N}$, the function counts the spikes of other neurons in the interval $t - r\delta : t - (1 + (r - 1)\delta)$.

$$ \phi_{r,\delta}^{m}(Y_{-\infty}^{t-1}) = \sum_{1+(r-1)\delta \le s \le r\delta} S_{t-s}^{m} \qquad (15) $$

In the experiments, we took {$\phi_{r,\delta}^m, m \in 1 : M, r \in 1 : 2$} with $\delta = 37$, which corresponds to a physical time of about 0.03 seconds.

3. LFP Power influence. Powerful LFP rhythms can increase or decrease the spiking probability. In the slow wave sleep phase theta waves ($j = 4, 5$) are negligible (Levenstein et al., 2019), so we discarded the theta power influences on the neurons. For each of the other $j$ ($0 : 9 \setminus \{4, 5\}$), we have added to the dictionary:

$$ \phi_{Pow}^{p,j}(Y_{-\infty}^{t-1}) = \sum_{s=1}^{.2^{J-j+1}} \max(|W_{t-s}^{p,j}|, w), \qquad (16) $$

which in fact represents the last wavelet coefficient in absolute value, truncated at the fixed positive constant $w$. In practice we took $w = 2$ and this value is so large that the wavelet coefficients have never been truncated in our experiments.

4. Phase lock. Phase lock phenomena are common in the recordings of cerebral activity. For the same $j$ index as before, we need to add 2 functions in the dictionary, one for the up phase and one for the down phase:

$$ \phi_{\uparrow}^{p,j}(Y_{-\infty}^{t-1}) = \sum_{s=1}^{2^{J-j}-1} \max(|W_{t-s}^{p,j}|, w)\mathbb{1}_{W_{t-s}^{p,j}>0} + \sum_{s=2^{J-j}-1}^{2^{J-j+1}-1} \max(|W_{t-s}^{p,j}|, w)\mathbb{1}_{W_{t-s}^{p,j}<0}, \qquad (17) $$

which in fact is the absolute value of the last wavelet coefficient if $t$ is in the corresponding up phase, and

$$ \phi_{\downarrow}^{p,j}(Y_{-\infty}^{t-1}) = \sum_{s=1}^{2^{J-j}-1} \max(|W_{t-s}^{p,j}|, w)\mathbb{1}_{W_{t-s}^{p,j}<0} + \sum_{s=2^{J-j}-1}^{2^{J-j+1}-1} \max(|W_{t-s}^{p,j}|, w)\mathbb{1}_{W_{t-s}^{p,j}>0}, \qquad (18) $$

which is the absolute value of the last wavelet coefficient if $t$ is in the corresponding down phase.
### 5.2.2 Dictionary for a wavelet coefficient.

In addition to the previous functions, we want here to model the global effect of a neuronal population, so that instead of using $\phi_{r,\delta}^m$ for specified neurons $m$, we grouped them in regions $R_i = \{m \in 1 : M|\text{neuron } m \text{ is in region } i\}$ and considered the cumulative effect of all the spikes in a region:

$$ \phi_{r,\delta}^{R_i}(Y_{-\infty}^{t-1}) = \sum_{m \in R_i} \sum_{1+(r-1)\delta \le s \le r\delta} S_{t-s}^m, \quad (19) $$

still with $r \in 1:2$ and $\delta = 37$ (i.e. 0.03 s).

Concerning LFP-LFP interactions, the choice of the dictionary depends on $j$, because of its neurobiological interpretation.

* For all of the rhythms we use $\phi_1$ so that the corresponding coefficient is estimating $\mu^{p,j}$.

* Again for all of them we added $\phi_{r,\delta}^{R_i}$ for $r = 1, 2$ and $\delta/f_S = 0.03$ seconds, to leave open the possibility of estimating the effect of the spikes on each LFP frequency.

* To study the effect of the power of LFP frequencies on other frequencies, we employ the dictionary function $\phi_{Pow}^{p',j'}(Y_{\infty}^{t-1})$. We employed tactically this fraction of the dictionary by adding the effect of one cycle of gamma power ($j' = 7, 8$) on Ripples ($j = 9$). Indeed, ripples are bursts of oscillatory activity whose frequency may be above 200 Hz. Therefore they should be mainly caught by the rhythm $j = 9$ but it is also possible that its effect leaks to lower, close frequencies.

* Finally, we left the highest number of degrees of freedom for the power-phase modulations. Since by definition of this phenomenon slower frequencies target faster frequencies (Buzsáki, 2009) for all $j \in \{0:9\}/\{4,5\}$ and all $j' < j$ we added the two dictionary functions $\phi_{\uparrow}^{p',j'}(Y_{\infty}^{t-1})$ and $\phi_{\downarrow}^{p',j'}(Y_{\infty}^{t-1})$.

A schematic view of the complete dictionary is given in Figure 8.

## 5.3 Results

We envision 5 numerical settings.

1. *Excitation*: The spiking of Neuron 2 triggers the spiking of Neuron 1.

2. *Phase Lock*: The low gamma ($j = 7$) up phase of the LFP signal excites Neuron 1.

3. *LFP power-phase modulation*: low gamma ($j = 7$) coefficients are diminished in the down phase of the high delta ($j = 3$) frequency.

4. *Spike-LFP average*: In this setting, we mimic a negative effect of neuronal spiking onto a specific LFP rhythm (high gamma, $j = 8$).

5. *Close to real data*: In the final simulation we have a complex connectivity graph, which is the fitting result of the cerebral sleep activity of a rat (see Section 6). The difference is that for simplicity of exposition, the spontaneous values and interaction strengths are restricted to few possible values.

The description of the parameters of each setting, as well as the corresponding connectivity graph and the average performance measures are reported in Figures 9
## Dictionary of Interactions

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Schematic of the chosen interaction dictionary. An arrow from one node to another indicates inclusion of the corresponding effect: a power effect if the arrow is solid (Spike-Spike (15), Spike-LFP (19), or LFP-LFP (16)); a phase effect if dotted (see (17) and (18)); or both power and phase effects if dashed. A loop from a node to itself signifies interactions within the same category.</figcaption></figure>

and 10. The first 4 settings correspond to the classic examples presented in Section 2, but with smaller parameters to increase the difficulty of detection. The last setting is closer to a real data behavior.

As shown by Figures 9 and 10, the L+LS method is the one giving the smallest error. More importantly, with less than 1 minute of observation, the method is able to recover perfectly the simplest graphs (DG=1). For the close-to-real data, for which the graph is much more complex, in less than one minute of observation, the estimation method is able to have about 0.5% of false positive and 19% of false negative for the LFP → Spikes interactions, whereas the LFP-LFP and Spike-Spike interactions are even more robustly reconstructed (at most 0.2% of false positives and no false negatives).

In conclusion to this Section, we simulated and reconstructed the functional connectivity patterns using a Haar basis, showing good performance even with only one minute of recordings. The codes are available at https://neuromod.gitlabpages.inria.fr/nm-hmmvar/. Other basis, with a distortion in the choice of the basis between the simulation and the reconstruction are explored in the Supplementary Material.

# 6 Experimental Data

After the numerical validation of the LASSO procedure, we test its capacity with the exact same dictionary on real data. We have almost an hour of recordings per condition in the real data set and we are therefore quite confident that we make almost no false positives, whereas the percentage of false negatives should be smaller (if the model is adequate for these data).

The objective of this section is to recover, from a high-quality dataset, relationships between neuronal populations and LFP which are well documented in the literature,
to see the ability of our procedure to screen complex data and detect automatically relevant interactions.

The data comes from the Buzsaki Lab dataset (Petersen et al., 2020; Girardeau et al., 2017a) and consists of the spikes and LFP recorded into two different brain regions: the hippocampus (HPC) and the basolateral amygdala (BLA). The neural data was recorded during training on a spatial aversive task in which a specific location is associated with an averse stimulus, and during sleep episodes before and after training. The experiment is explained in detail in (Girardeau et al., 2017b). The data used here consists in the spikes and LFPs recorded during a particular sleep phase called non-REM sleep or "slow wave sleep", characterized by widespread slow-frequency (0,5-4 Hz) oscillations in the cortex and ripples in the hippocampus.

In this framework, we treat the data as follows:

1. We isolate the section of the experiment that we want to analyze. We consider the recordings performed before (preRUN) and after (postRUN) the aversive spatial task in Rat08 on the 6th day of the experiment. As the rat switched multiple times between slow wave sleep, REM sleep, and wakefulness, we isolate and then join the segments of slow wave sleep data.

2. We want to infer the functional connectivity via the above L+LS (Lasso + Least Squares Step) procedure, with the same set-up as in the simulations, that is considering $J = 10$ as the finest scale and employing the same dictionary of functions. The hyperparameters are similar to Section 5, that is $\gamma_S = 0.6$ and $\gamma_W = 0.3$.

3. Since the LFP sampling frequency of the dataset is 1250 Hz, every block of wavelet coefficients has a duration of $2^{10}/1250 = 1024/1250 = 0.8192$ seconds. Because of the different duration of each sleep period, we decided to consider 74 segments of $N = 70$ wavelet blocks for the preRUN slow wave sleep and 72 segments of $N = 70$ wavelet blocks for the postRUN slow wave sleep. This leads to an observation time of $\Upsilon = 57.344$ seconds for every segment and to aggregated observation times of $\approx 1\text{h}10\text{m}40\text{s}$ for the preRUN sleep and $\approx 1\text{h}8\text{m}45\text{s}$ for the postRUN sleep. The longer the observation time, the higher is the computational cost but we gain in estimator performance.

4. From this time interval we isolate the spikes and the LFP signals coming from the two brain regions of interest: the hippocampus (HPC) and baso-lateral amygdala (BLA).

The resulting connectivity graphs are displayed in Figure 11. In total, 156 neurons (with two additional global nodes for the HPC and BLA neuronal populations) and 10\*2 LFP rhythms (10 rhythms per region) are used for the graph reconstruction. Thus the complete graph would be on 178 nodes, but for readability, we chose to show only the variables which are sources or targets of a nonzero, i.e. the nodes connected by an arrow.

One of the main differences in the connectivity graphs of slow wave sleep preceding vs. following the aversive task is the emergence of a new positive connection between hippocampal ripples (R1-HPC) and a BLA neuron (neuron 88). This finding echoes the results of (Girardeau et al., 2017b), in which they show that during sleep following the task, there is an increase in the modulation of a small subpopulation of BLA neurons
by hippocampal ripples. The graph also shows a higher effect of beta frequency in the BLA on BLA neurons, in the postRUN, probably triggered by training on the task. It is a new finding that will be explored in more detail along with other connectivity changes within and across the two structures in follow-up works.

# 7 Conclusion

In conclusion, we proposed a model of spike and LFP rhythm interactions. We identified conditions under which the model exists in a block stationary version. We proposed a weighted Lasso estimator of the interactions, for which we proved an oracle inequity. In theory, our method can accept any dictionary, but the interpretability of the results depends on its choice. We decided to fine-tune the dictionary to unravel the existence and magnitude of classical LFP-LFP, Spike-Spike and LFP-Spike interactions of the literature in neuroscience.

In future works, by incorporating successively new type of interactions in the dictionary, our method will be able to assess whether these new interactions are statistically relevant and can improve the understanding of the functional dynamics of the network. In particular, we will further investigate the data of (Girardeau et al., 2017b) to establish statistical evidence of the connectivity changes within and across experimental days.

The method requires calibration of the hyperparameters $\gamma_S$ and $\gamma_W$. We performed this calibration for the Haar basis and the chosen dictionary (see also Supplementary Section 3 for other bases) through extensive, one-time simulations. Cross-validation is another popular method of calibration in statistics, that calibrates more finely the parameters for each data set and each dictionary, but with a huge computational cost at each data set analysis. It could be a future line of development of the method.

Finally, it would be interesting from a mathematical point of view, to see if we can adapt our model to continuous time, mixing point processes in continuous time and continuous Gaussian processes. Moreover the possibility to have a continuum of frequencies would be of interest to fine-tune the detection to meaningful rhythms of the neuroscience literature.

**Supplementary information.** The Supplementary Material contains the additional proofs and the precise formulas for the performance measures, as well as a numerical assessment of the performance of non Haar basis in reconstruction.

**Acknowledgements.** This research was supported by the French government, through CNRS (eXplAIn team), the UCA<sup>Jedi</sup> and 3IA Côte d'Azur Investissements d'Avenir managed by the National Research Agency (ANR-15 IDEX-01 and ANR-19-P3IA-0002), INSERM (Atip-Avenir), directly by the ANR project ChaMaNe (ANR-19-CE40-0024-02), by the interdisciplinary Institute for Modeling in Neuroscience and Cognition (NeuroMod) and Ville de Paris (Emergence(s)). We would like to thank the referee, whose comments help to greatly improve this manuscript.

# References

Albo, Z., Di Prisco, G.V., Chen, Y., Rangarajan, G., Truccolo, W., Feng, J., vertes,
R.P., Ding, M.: Is partial coherence a viable technique for identifying generators of neural oscillations? Biological Cybernetics **90**, 318–326 (2004)

Abbaspourazad, H., Hsieh, H.-L., Shanechi, M.M.: A multiscale dynamical modeling and identification framework for spike-field activity. IEEE Transactions on Neural Systems and Rehabilitation Engineering **27**(6), 1128–1138 (2019)

Bonnet, A., Dion, C., Gindraud, F., Lemler, S.: Neuronal network inference and membrane potential model using multivariate hawkes processes. Journal of Neuroscience Methods (2022)

Baranowski, R., Fryzlewicz, P.: Multiscale autoregression on adaptively detected timescales. In: Proceedings of the 12th International Conference of the ERCIM WG on Computational and Methodological Statistics (CMStatistics 2019) (2020)

Brémaud, P., Massoulié, L.: Stability of nonlinear hawkes processes. The Annals of Probability **24** (1996) https://doi.org/10.1214/aop/1065725193

Baccala, L.A., Sameshima, K.: Partial directed coherence: a new concept in neural structure determination. Biological Cybernetics **84**(6), 463–474 (2001)

Buzsáki, Gy and Anastassiou, Costas A and Koch, Christof, v.n.p.y.p. journal=Nature reviews neuroscience: The origin of extracellular fields and currents—eeg, ecog, lfp and spikes

Buzsáki, G.: Rhythms of The Brain, p. 448 (2009). https://doi.org/10.1093/acprof:oso/9780195301069.001.0001

Buzsàki, G., Wang, X.-J.: Mechanisms of gamma oscillations. Annual Review of Neuroscience **35**(1), 203–225 (2012) https://doi.org/10.1146/annurev-neuro-062111-150444 PMID: 22443509

Cohen, A., Daubechies, I., Feauveau, J.-C.: Biorthogonal bases of compactly supported wavelets. Communications on Pure and Applied Mathematics **45**(5), 485–560 (1992)

Cohen, M.X.: A better way to define and describe morlet wavelets for time-frequency analysis. NeuroImage **199**, 81–86 (2019)

Carstensen, L., Sandelin, A., Winther, O., Hansen, N.R.: Multivariate hawkes process models of the occurrence of regulatory elements. BMC Bioinformatics **11**, 1–19 (2010)

Ding, M., Chen, Y., Bressler, S.L.: Granger causality: basic theory and application to neuroscience. Handbook of time series analysis: recent theoretical developments and applications, 437–460 (2006)
Dion, C., Lemler, S.: Nonparametric drift estimation for diffusions with jumps driven by a hawkes process. Statistical Inference for Stochastic Processes (2020) https://doi.org/10.1007/s11203-020-09213-5

De Leon, A.R., Chough, K.C.: Analysis of Mixed Data: Methods & Applications. CRC Press, ??? (2013)

Doukhan, P., Wintenberger, O.: Weakly dependent chains with infinite memory. arXiv (2007). https://doi.org/10.48550/ARXIV.0712.3231 . https://arxiv.org/abs/0712.3231

Frese, T., Bouman, C.A., Sauer, K.: Adaptive wavelet graph model for bayesian tomographic reconstruction. IEEE Transactions on Image Processing 11(7), 756-770 (2002)

Faes, L., Erla, S., Nollo, G.: Measuring connectivity in linear multivariate processes: Definitions, interpretation, and practical analysis. Computational and Mathematical Methods in Medicine (ID 140513) (2012)

Gubner, J.A., Chang, W.-B.: Wavelet transforms for discrete-time periodic signals. Signal Processing 42(2), 167-180 (1995)

Girardeau, G., Inema, I., Buzsaki, G.: Simultaneous large-scale recordings in dorsal hippocampus, basolateral amygdala and neighbouring deep nuclei and structures in rats performing a spatial aversive task and sleeping. (2017). www.CRCNS.orghttp://dx.doi.org/10.6080/K0MS3QXD

Girardeau, G., Inema, I., Buzsáki, G.: Reactivations of emotional memory in the hippocampus-amygdala system during sleep. Nature Neuroscience 20(11), 1634-1642 (2017)

Galves, A., Löcherbach, E.: Modeling networks of spiking neurons as interacting processes with memory of variable length. arXiv (2015). https://doi.org/10.48550/ARXIV.1502.06446 . https://arxiv.org/abs/1502.06446

Gong, X., Li, W., Liang, H.: Spike-field granger causality for hybrid neural data analysis. Journal of neurophysiology 122(2), 809-822 (2019)

Granger, C.W.J.: Testing for causality: A personal viewpoint. Journal of Economic Dynamics and Control 2, 329-352 (1980) https://doi.org/10.1016/0165-1889(80)90069-X

Hansen, N.R., Reynaud-Bouret, P., Rivoirard, V.: Lasso and probabilistic inequalities for multivariate point processes. Bernoulli 21(1) (2015) https://doi.org/10.3150/13-bej562

Jensen, O.: Reading the hippocampal code by theta phase-locking. Trends in Cognitive
Sciences **9**(12), 551-553 (2005) https://doi.org/10.1016/j.tics.2005.10.003

Kobayashi, R., Kurita, S., Kurth, A., Kitano, K., Mizuseki, K., Diesmann, M., Richmond, B.J., Shinomoto, S.: Reconstructing neuronal circuitry from parallel spike trains. Nature Communications **10**(1), 4468 (2019)

Kwon, H.-H., Lall, U., Khalil, A.F.: Stochastic simulation model for nonstationary time series using an autoregressive wavelet decomposition: Applications to rainfall and temperature. Water Resources Research **43**(5) (2007) https://doi.org/10.1029/2006WR005258
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2006WR005258

Kreiss, A., Mammen, E., Polonik, W.: Nonparametric inference for continuous-time event counting and link-based dynamic network models. Electronic Journal of Statistics **13**, 2764–2829 (2019)

Lisman, J., Buzsáki, G.: A neural coding scheme formed by the combined function of gamma and theta oscillations. Schizophrenia Bulletin **34**(5), 974–980 (2008)

Levenstein, D., Buzsáki, G., Rinzel, J.: Nrem sleep in the rodent neocortex and hippocampus reflects excitable dynamics. Nature communications **10**(1), 2478 (2019)

Lepage, K.Q., Gregoriou, G.G., Kramer, M.A., Aoi, M., Gotts, S.J., Eden, U.T., Desimone, R.: A procedure for testing across-condition rhythmic spike-field association change. Journal of neuroscience methods **213**(1), 43–62 (2013)

Lisman, J., Jensen, O.: The theta-gamma neural code. Neuron **77**, 1002–16 (2013)
https://doi.org/10.1016/j.neuron.2013.03.007

Lambert, R.C., Tuleau-Malot, C., Bessaih, T., Rivoirard, V., Bouret, Y., Leresche, N., Reynaud-Bouret, P.: Reconstructing the functional connectivity of multiple spike trains using hawkes models. Journal of Neuroscience Methods **297**, 9–21 (2018)

Murtagh, F., Starck, J.L., Renaud, O.: On neuro-wavelet modeling. Decision Support Systems **37**(4), 475–484 (2004) https://doi.org/10.1016/S0167-9236(03)00092-7
Data mining for financial decision making

Ost, G., Reynaud-Bouret, P.: Sparse space-time models: Concentration Inequalities and Lasso. arXiv (2018). https://doi.org/10.48550/ARXIV.1807.07615. https://arxiv.org/abs/1807.07615

Petersen, P.C., Hernandez, M., Buzsáki, G.: The buzsaki lab databank – public electro-
physiological datasets from awake animals (2020) https://doi.org/10.5281/zenodo.
4307883

Raskutti, G., Wainwright, M.J., Yu, B.: Restricted eigenvalue properties for correlated
gaussian designs. Journal of Machine Learning Research **11**, 2241–2259 (2010)

Song, C.Y., Hsieh, H.-L., Pesaran, B., Shanechi, M.M.: Modeling and inference methods for switching regime-dependent dynamical systems with multiscale neural observations. Journal of Neural Engineering **19**(6), 066019 (2022)

Teleńczuk, B., Dehghani, N., Le Van Quyen, M., Cash, S.S., Halgren, E., Hatsopoulos, N.G., Destexhe, A.: Local field potentials primarily reflect inhibitory neuron activity in human and monkey cortex. Scientific Reports **7**(1), 40211 (2017)

Vinck, M., Battaglia, F., Womelsdorf, T., Pennartz, C.: Improved measures of phase-coupling between spikes and the local field potential. Journal of computational neuroscience **33**, 53–75 (2011) https://doi.org/10.1007/s10827-011-0374-4

Geer, S., Bühlmann, P., Zhou, S.: The adaptive and the thresholded lasso for potentially misspecified models (and a lower bound for the lasso). Electron. J. Stat. **5**, 688–749 (2011)

Willsky, A.S.: Multiresolution markov models for signal and image processing. Proceedings of the IEEE **90**(8), 1396–1458 (2002) https://doi.org/10.1109/JPROC.2002.800717

Wang, C., Shanechi, M.M.: Estimating multiscale direct causality graphs in neural spike-field networks. IEEE Transactions on Neural Systems and Rehabilitation Engineering **27**(5), 857–866 (2019)

Zhang, Y., Liu, B., Ji, X., Huang, D.: Classification of eeg signals based on autoregressive model and wavelet packet decomposition. Neural Processing Letters **45**, 365–378 (2017)

Zhu, Y., Wang, R.: Research on phase synchronization with spike-lfp coherence analysis. Neurocomputing **168**, 655–660 (2015)

Zhang, L., Wu, X., Pan, Q., Zhang, H.: Multiresolution modeling and estimation of multisensor data. IEEE transactions on Signal Processing **52**(11), 3170–3182 (2004)
<figure><img src="image_11.png" /><figcaption>Fig. 9 Performance evaluation for settings 1 to 3, from top to bottom. In each figure, up left: specific parameters for the setting; down left: average errors over 100 simulations; right: true connectivity graph.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10** Performance evaluation for settings 4 (top) and 5 (bottom). In each figure, up left: specific parameters for the setting; down left: average errors over 100 simulations; right: true connectivity graph. Color codes of the nodes for setting 5 are given in Figure 11 on real data.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11** Fitted connectivity graphs in aggregated segments of slow wave sleep cerebral activity before (left) and after (right) the aversion task. The HPC neurons and LFPs are shown in light and dark green, respectively. The BLA neurons and LFPs are shown in yellow and brown. BLA population activity appears in orange, as HPC population activity is not estimated as a node of relevance. Positive interactions are shown as red arrows, and negative interactions as blue arrows. The width of the arrow edges represents the strength (in $ℓ_1$ norm) of the interaction. The style of the arrows (solid/dashed) represents the type of interaction. As a reminder, "Power" means that the power, $|W_t^{p,j}|$, has an influence on the target rhythm or neuron. "Phase locking" means that the phase of the brain oscillation has an influence on the activity of the target neuron(s) or rhythm. The black oval is there to highlight a major difference between the two graphs. Note that the position of each node is fixed between the two graphs but if a node is connected in preRUN and not connected in postRUN (or reciprocally), it will only appear in the preRUN graph (reciprocally in the postRUN graph).</figcaption></figure>