Olivier David
email: odavid@fil.ion.ucl.ac.uk

Lee Harrison

Karl J Friston

Modelling event-related responses in the brain

The aim of this work was to investigate the mechanisms that shape evoked electroencephalographic (EEG) and magneto-encephalographic (MEG) responses. We used a neuronally plausible model to characterise the dependency of response components on the models parameters. This generative model was a neural mass model of hierarchically arranged areas using three kinds of inter-area connections (forward, backward and lateral).

We investigated how responses, at each level of a cortical hierarchy, depended on the strength of connections or coupling. Our strategy was to systematically add connections and examine the responses of each successive architecture. We did this in the context of deterministic responses and then with stochastic spontaneous activity. Our aim was to show, in a simple way, how event-related dynamics depend on extrinsic connectivity. To emphasise the importance of nonlinear interactions we tried to disambiguate the components of event-related potentials (ERPs) or event-related fields (ERFs) that can be explained by a linear superposition of trial-specific responses and those engendered nonlinearly (e.g. by phase-resetting).

Our key conclusions were; (i) when forward connections, mediating bottom-up or extrinsic inputs, are sufficiently strong, nonlinear mechanisms cause a saturation of excitatory interneuron responses. This endows the system with an inherent stability that precludes nondissipative population dynamics. (ii) The duration of evoked transients increases with the hierarchical depth or level of processing. (iii) When backward connections are added, evoked transients become more protracted, exhibiting damped oscillations. These are formally identical to late or endogenous components seen empirically. This suggests late components are mediated by reentrant dynamics within cortical hierarchies. (iv) Bilateral connections produce similar effects to backward connections but can also mediate zero-lag phase-locking among areas. (v) Finally, with spontaneous activity, ERPs/ERFs can arise from two distinct mechanisms: For low levels of (stimulus related and ongoing) activity the systems response conforms to a quasi-linear superposition of separable responses to the fixed and stochastic inputs. This is consistent with classical assumptions that motivate trial averaging to suppress spontaneous activity and disclose the ERP/ERF. However, when activity is sufficiently high, there are nonlinear interactions between the fixed and stochastic inputs. This interaction is expressed as a phase-resetting and represents a qualitatively different explanation for the ERP/ERF.

Introduction

Classical event-related potentials (ERPs) and event-related fields (ERFs) have been used for decades as putative electrophysiological correlates of perceptual and cognitive operations.

However, the exact neurobiological mechanisms underlying their generation are largely unknown. Recently there has been a special interest in the distinction between evoked and induced responses. Evoked responses are disclosed by conventional averaging procedures, whereas the latter usually call for single trial analyses of induced oscillations. In this paper we used neuronal simulations to examine the mechanisms that underpin ERPs/ERFs. In a companion paper we will examine induced responses using time-frequency analyses and other transforms of single trial data.

Neural-mass models

The complexity of neural networks generating MEG/EEG signals [START_REF] Defelipe | Microstructure of the neocortex: comparative aspects[END_REF][START_REF] Thomson | Synaptic interactions in neocortical local circuits: dual intracellular recordings in vitro[END_REF] is considerable. This means that MEG/EEG observation models rely upon simplifying assumptions and empirical priors [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF][START_REF] Freeman | Models of the dynamics of neural populations[END_REF][START_REF] Valdes | Nonlinear EEG analysis based on a neural mass model[END_REF][START_REF] Lopes Da Silva | Model of brain rhythmic activity. The alpha-rhythm of the thalamus[END_REF][START_REF] Robinson | Prediction of electroencephalographic spectra from neurophysiology[END_REF][START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity?[END_REF][START_REF] Van Rotterdam | A model of the spatial-temporal characteristics of the alpha rhythm[END_REF][START_REF] Wendling | Relevance of nonlinear lumpedparameter models in the analysis of depth-EEG epileptic signals[END_REF]. The primary aim of this paper is to describe a candidate forward model and establish its face validity. This model was designed to reproduce responses seen empirically and enable mechanistic enquiries into the generation of evoked and induced responses. This is the focus of the current paper. However, we will also use this model in a forthcoming paper as an observation model, allowing its parameters to be inferred from real data (David et al., 2004b). In this context, face validity is especially important.

Neural-mass models of MEG/EEG usually comprise cortical macro-columns, which can be treated as surrogates for cortical areas and, sometimes, thalamic nuclei. These models use a small number of state variables to represent a neuronal population mean state. This approach, referred to loosely as a mean-field approximation, is efficient when determining the steadystate behaviour of neuronal systems but its accuracy in a dynamic or non-stationary context is less established [START_REF] Haskell | Population density methods for large-scale modelling of neuronal networks with realistic synaptic kinetics: cutting the dimension down to size[END_REF]. However, we will assume that the mean field approximation is sufficient for our purposes. The majority of neural mass models of MEG/EEG have been designed to generate alpha rhythms (Lopes da [START_REF] Lopes Da Silva | Model of brain rhythmic activity. The alpha-rhythm of the thalamus[END_REF][START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF][START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity?[END_REF][START_REF] Van Rotterdam | A model of the spatial-temporal characteristics of the alpha rhythm[END_REF]. Recent studies have shown that it is possible to reproduce the whole spectrum of MEG/EEG oscillations, using appropriate values of model parameters [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF][START_REF] Robinson | Prediction of electroencephalographic spectra from neurophysiology[END_REF]. In addition these models have been used to test specific hypotheses about brain function, e.g. focal attention [START_REF] Suffczynski | Computational model of thalamo-cortical networks: dynamical control of alpha rhythms in relation to focal attention[END_REF]. Pathological activity such as epilepsy can also be emulated. This means, in principle, generative models of the sort employed above could be used to characterise the pathophysiological mechanisms underlying seizure activity [START_REF] Robinson | Dynamics of large-scale brain activity in normal arousal states and epileptic seizures[END_REF][START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF].

To date, modelling event-related activity using neural mass models has received much less attention. An early attempt, in the context of visual ERPs, showed that it was possible to emulate ERP-like damped oscillations [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF]. A more sophisticated thalamocortical model has been used to simulate event-related synchronisation (ERS) and eventrelated desynchronisation (ERD), commonly found in the alpha band [START_REF] Suffczynski | Computational model of thalamo-cortical networks: dynamical control of alpha rhythms in relation to focal attention[END_REF]. Finally, it has been shown that model parameters can be adjusted to fit real ERPs [START_REF] Rennie | Unified neurophysical model of EEG spectra and evoked potentials[END_REF]. These studies [START_REF] Rennie | Unified neurophysical model of EEG spectra and evoked potentials[END_REF][START_REF] Suffczynski | Computational model of thalamo-cortical networks: dynamical control of alpha rhythms in relation to focal attention[END_REF], emphasise the role of the thalamo-cortical interactions by modelling the cortex as a single compartment.

Hierarchical models

It is well-known that the cortex has a hierarchical organisation [START_REF] Crick | Constraints on cortical and thalamic projections: the no-strong-loops hypothesis[END_REF][START_REF] Felleman | Distributed hierarchical processing in the primate cerebral cortex[END_REF], comprising bottom-up, top-down and lateral processes that can be understood from an anatomical and cognitive perspective [START_REF] Engel | Dynamic predictions: oscillations and synchrony in topdown processing[END_REF].

We have previously discussed the importance of hierarchical processes, in relation to perceptual inference in the brain, using the intimate relationship between hierarchical models and empirical Bayes [START_REF] Friston | Functional integration and inference in the brain[END_REF]. The current work was more physiologically motivated.

Using a hierarchical neural mass model, we were primarily interested in the effects, on eventrelated MEG/EEG activity, of connections strengths, and how these effects were expressed at different hierarchical levels. In addition, we were interested in how nonlinearities in these connections might be expressed in observed responses.

The neuronal model described below embodies many neuroanatomic and physiological constraints which lend it a neuronal plausibility. It has been designed to (i) explore emergent behaviours that may help understand empirical phenomena and, critically, (ii) as the basis of dynamic observation models. Although the model comprises coupled systems, the coupling is highly asymmetric and heterogeneous. This contrasts with homogenous and symmetrically coupled map lattices (CML) and globally coupled maps (GCM) encountered in more analytic treatments. Using the concepts of chaotic dynamical systems, GCMs have motivated a view of neuronal dynamics that is cast in terms of high-dimensional transitory dynamics among 'exotic' attractors [START_REF] Tsuda | Toward an interpretation of dynamic neural activity in terms of chaotic dynamical systems[END_REF]. Much of this work rests on uniform coupling, which induces a synchronisation manifold, around which the dynamics play. The ensuing chaotic itinerancy has many intriguing aspects that can be related to neuronal systems [START_REF] Breakspear | Modulation of excitatory synaptic coupling facilitates synchronization and complex dynamics in a biophysical model of neuronal dynamics[END_REF][START_REF] Kaneko | Chaotic itinerancy[END_REF]. However, the focus of this work is not chaotic itinerancy but chaotic transience (the transient dynamics evoked by perturbations to the systems state), in systems with asymmetric coupling. This focus precludes much of the analytic treatment available for GCMs (but see [START_REF] Jirsa | Spatiotemporal pattern formation in neural systems with heterogeneous connection topologies[END_REF] for an analytical description of coherent pattern formation in a spatially continuous neural system with a heterogeneous connection topology). However, as we hope to show, simply integrating the model, to simulate responses, can be a revealing exercise.

Mechanisms of ERP/ERF generation; linear or nonlinear?

It is generally held that an ERP/ERF is the result of averaging a set of discrete stimulusevoked brain transients [START_REF] Coles | Event-related brain potentials: An introduction[END_REF]. However, several groups [START_REF] Jansen | Phase synchronization of the ongoing EEG and auditory EP generation[END_REF][START_REF] Klimesch | Phaselocked alpha and theta oscillations generate the P1-N1 complex and are related to memory performance[END_REF][START_REF] Kolev | Analysis of phase-locking is informative for studying eventrelated EEG activity[END_REF][START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF] have suggested that some ERP/ERF components might be generated by stimulus-induced changes in ongoing brain dynamics. This is consistent with views emerging from several neuroscientific fields, suggesting that phase-synchronization, of ongoing rhythms, across different spatio-temporal scales mediates the functional integration necessary to perform higher cognitive tasks [START_REF] Penny | Event-related brain dynamics[END_REF][START_REF] Varela | The brainweb: phase synchronization and large-scale integration[END_REF]. In brief, a key issue is the distinction between processes that do and do not rely on phase-resetting of ongoing spontaneous activity. Both can lead to the expression of ERP/ERF components but their mechanisms are very different.

EEG and MEG signals are effectively ergodic and cancel when averaged over a sufficient number of randomly chosen epochs. The fact that ERPs/ERFs exhibit systematic waveforms, when the epochs are stimulus locked, suggests either a reproducible stimulus-dependent modulation of amplitude or phase-locking of ongoing MEG/EEG activity [START_REF] Tass | Stochastic phase resetting of stimulus-locked responses of two coupled oscillators: transient response clustering, synchronization, and desynchronization[END_REF]. The key distinction, between these two explanations, is whether the stimulus-related component interacts with ongoing or spontaneous activity. If there is no interaction the spontaneous component will be averaged out, because it has no consistent phase-relationship with stimulus onset. Conversely if there is an interaction, dominant frequencies of the spontaneous activity must experience a phase-change, so that they acquire a degree of phaselocking to the stimulus. Note that phase-resetting is a stronger-requirement than induced oscillations. It requires any induced dynamics to be phase-locked in peristimulus time. In short, phase-resetting is explicitly nonlinear and implies an interaction between stimulusrelated response and ongoing activity. Put simply, this means that the event-related response depends on ongoing activity. This dependency can be assessed with the difference between responses elicited with and without the stimulus (if we could reproduce exactly the same ongoing activity). In the absence of interactions there will be no difference. Any difference implies nonlinear interactions. Clearly this cannot be done empirically but it can be pursued using neuronal simulations.

The secondary aim of the current work was to use realistic neural-mass models of hierarchically organised cortical areas to see whether phase-resetting is an emergent phenomenon and a plausible candidate for causing ERPs/ERFs. Phase-resetting is used in this paper as an interesting example of nonlinear responses that have been observed empirically. We use it to show that nonlinear mechanisms can be usefully explored with neuronal models of the sort developed here. In particular, static nonlinearities, in neuronal mass models, are sufficient to explain phase-resetting. Phase-resetting represents nonlinear behaviour because, in the absence of amplitude changes, phase-changes can only be mediated in a nonlinear way. This is why phase-synchronization plays a central role in detecting nonlinear coupling among sources [START_REF] Breakspear | Nonlinear phase desynchronization in human electroencephalographic data[END_REF][START_REF] Tass | Stochastic phase resetting of stimulus-locked responses of two coupled oscillators: transient response clustering, synchronization, and desynchronization[END_REF].

Overview

This paper is structured as follows. In the first section we introduce the hierarchical neural mass model used in the remaining sections. It is based on previous neuroanatomic studies by [START_REF] Felleman | Distributed hierarchical processing in the primate cerebral cortex[END_REF] and work by Jansen and Rit on modelling MEG/EEG data [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF]. In the second section we demonstrate the basic behaviour of the model, by successive elaboration of a cortical hierarchy. We start with forward connections and then add backward and lateral connections. The goal of this approach was to provide an intuitive understanding of MEG/EEG like dynamics generated by coupled nonlinear systems. These simulations were performed in the absence of spontaneous activity. In the third section we examine the interaction between evoked and spontaneous activity, using a representative hierarchical architecture established in the previous section.

Finally, we discuss the potential benefits of this modelling approach, for the study of measured MEG/EEG activity.

Hierarchical models of event-related MEG/EEG activity

Cortico-cortical connections

Although neural mass models originated in the early 1970's [START_REF] Freeman | Models of the dynamics of neural populations[END_REF][START_REF] Lopes Da Silva | Model of brain rhythmic activity. The alpha-rhythm of the thalamus[END_REF][START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF], none have addressed explicitly the hierarchical nature of cortical organisation. The minimal model we propose, which accounts for directed extrinsic connections, uses the rules in [START_REF] Felleman | Distributed hierarchical processing in the primate cerebral cortex[END_REF]. Extrinsic connections are connections that traverse white matter and connect cortical regions (and subcortical structures). These rules, based upon a tri-partitioning of the cortical sheet (into supra-, infra-granular layers and granular layer 4), have been derived from experimental studies of cat visual cortex. We will assume that they can be generalised to the whole cortex.

The ensuing model is general, and can be used to model various cognitive paradigms (David et al., 2004b). However, variability among different cytoarchitectonic regions is restricted to differences in physiological parameters, under the same microcircuitry. Under this simplifying assumption, the connections can be defined as in Figure 1: (i) Bottom-up or forward connections originate in agranular layers and terminate in layer 4. (ii) Top-down or backward connections only connect agranular layers. (iii) Lateral connections originate in agranular layers and target all layers. All these long-range or extrinsic cortico-cortical connections are excitatory and are mediated through the axons of pyramidal cells.

Figure 1 about here Although the thalamocortical connections have been the focus of several modelling studies, they represent a minority of extrinsic connections: in contrast, it is thought that at least 99% of axons in white matter link cortical areas of the same hemisphere [START_REF] Abeles | Corticonics: neural circuits of the cerebral cortex[END_REF]. For this reason, and for simplicity, we do not include the thalamic nuclei in our model. However, they can be included for any application where the role of the thalamus (or other subcortical structure) is thought important.

Jansen model of a cortical area

The neocortex is commonly described as a 6-layered structure [START_REF] Defelipe | Microstructure of the neocortex: comparative aspects[END_REF]. Spiny neurons (pyramidal cells and spiny stellate cells) and smooth neurons are the two major groups of cortical neurons. The majority of cortical neurons are pyramidal cells that are found in layers 2 to 6. Most spiny stellate cells are interneurons that are located in the middle cortical layers. Smooth neurons are essentially GABAergic interneurons distributed in all layers. In general, cortical neurons are thought to be organized into multiple, small repeating microcircuits. In spite of cortical heterogeneity, a common basic microcircuit has emerged.

Its skeleton is formed by a pyramidal cell which receives excitatory inputs that originate from extrinsic afferent systems and spiny cells. Inhibitory inputs originate mostly from GABAergic interneurons. These micro-anatomical characteristics have been found in all cortical areas and species examined so far and, therefore, they can be considered as fundamental aspects of cortical organization [START_REF] Defelipe | Microstructure of the neocortex: comparative aspects[END_REF].

The Jansen model [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] follows the microcircuitry described above to emulate the MEG/EEG activity of a cortical area. Consequently, it is particularly well suited to embed in a hierarchical structure. A cortical area is modelled by three subpopulations. A population of excitatory pyramidal (output) cells receives inputs from inhibitory and excitatory populations of interneurons, via intrinsic connections (intrinsic connections are confined to the cortical sheet). Within this model, excitatory interneurons can be regarded as spiny stellate cells found predominantly in layer 4 and in receipt of forward connections [START_REF] Miller | Understanding layer 4 of the cortical circuit: a model based on cat V1[END_REF]. Excitatory pyramidal cells and inhibitory interneurons will be considered to occupy agranular layers and receive backward and lateral inputs.

We have described dynamics of these three subpopulations previously. We will review the model briefly but refer interested reader to [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF] for more details. The main difference (c.f. [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF] is that ERPs/ERFs are modelled as small perturbations around the resting potential. Therefore, all the variables below are zero-mean, centred on the resting state we assumed to be zero. The evolution of population dynamics rests on two operators. The first transforms p, the average density of pre-synaptic input arriving at the population, into v, the average postsynaptic membrane potential. This is modelled by the linear transformation

p h v  
, where  denotes convolution and h is the impulse response or kernel

           0 0 0 ) exp( t t t Ht t h   .
(1)

The excitatory (e) and inhibitory (i) kernels, h e and h i respectively, are parameterised by H e,i and i e,  modelling specific properties of inhibition and excitation. The parameters H e,i control the maximum post-synaptic potential and i e,  are lumped time-constants of passive membrane currents and other spatially distributed delays in the dendritic tree. The second operator transforms the average membrane potential of each subpopulation into an average firing rate. This is assumed to be instantaneous and is described by the sigmoid function

    0 0 exp 1 2 e rv e v S     (2)
where e 0 and r are parameters that determine its shape (e.g. voltage sensitivity). It is this function that endows the simulation with nonlinear behaviours that are critical for phenomena like phase-resetting.

Interactions, among the different subpopulations, depend on the constants i  , which control the strength of intrinsic connections and the total number of synapses expressed by each subpopulation. The relative values of these constants are fixed, using anatomical information from the literature, as described in [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF]:

1 4 3 1 2 25 . 0 , 8 . 0        
. The Jansen model is summarised in Figure 2.

Figure 2 about here

We assume MEG/EEG signals are a linear mixture of the average depolarisation of pyramidal cells. This mixture depends upon source lead-fields that model the spatial distribution of the ensuing electromagnetic fields [START_REF] Baillet | Electromagnetic brain mapping[END_REF]. We further assume the depolarisation of pyramidal cells is proportional to the cortical current source densities, which are estimated using inverse solutions [START_REF] Baillet | Electromagnetic brain mapping[END_REF]. Therefore, we consider the indirectly observed 'output' of the Jansen model to be the depolarisation of pyramidal cells.

For simplicity, we ignore the effects of instrumental amplifiers and assume that MEG/EEG cortical current densities can be estimated precisely and refer to them as "MEG/EEG signals". Thus, unless otherwise specified, signals are simply the depolarisation of pyramidal cells.

For given synaptic kernels h and sigmoid functions S, the Jansen model can reproduce a large variety of MEG/EEG-like waveforms [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF][START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF][START_REF] Wendling | Relevance of nonlinear lumpedparameter models in the analysis of depth-EEG epileptic signals[END_REF]. Although we have presented the dynamics in terms of a convolution operator and static nonlinearly, the integration of the model actually proceeds using the equivalent differential equations. These are provided in the Appendix. Being able to formulate the model in terms of differential equations is important for extending the neural mass formulation used here to a true mean field approximation using the Fokker-Plank formalism [START_REF] Groff | Stochastical aspects of neuronal dynamics: Fokker-Planck approach[END_REF][START_REF] Fourcaud | Dynamics of the firing probability of noisy integrate-and-fire neurons[END_REF]. This will be the subject of a future communication (Harrison et al in preparation).

Hierarchical extension of the Jansen model

Using the connection rules above, it is straightforward to construct hierarchical cortico- As discussed in section 2.1, we consider three types of extrinsic connections (Figure 3): forward, backward, and lateral. The strength of each type of connection is controlled by a coupling parameter a: a F for forward, a B for backward and a L for lateral. We model propagation delays for these connections (see Appendix).

cortical
Using these connections, hierarchical cortical models for MEG/EEG can be constructed to test various hypotheses, and represent examples of dynamic causal models [START_REF] Friston | Dynamic causal modelling[END_REF]. The causal model here is a multiple-input multiple-output system that comprises m inputs and l outputs with one output per region. The m inputs correspond to designed causes (e.g., stimulus functions encoding the occurrence of events) or stochastic processes modelling background activity. In principle, each input could have direct access to every region. However, in practice the effects of inputs are usually restricted to a single input region, usually the lowest in the hierarchy. Each of the l regions produces a measured output that corresponds to the MEG/EEG signal. Each region has five ( 1 , , , ,

  i e i e

H

) intrinsic parameters such as the membrane time constants described above. These play a crucial role in generating regional responses. However, in the present study, we will consider them fixed and focus on the extrinsic coupling parameters or effective connectivity. These are the matrices C, A F , A B and A L that contain the coupling parameters c, a F , a B and a L . The values of these parameters, used in the following simulations, are provided in the Appendix and in

Figure captions

.

Input-output behaviour

In this section we characterise the input-output behaviour of a series of canonical networks in terms of their impulse response functions. This is effectively the response (mean depolarisation of pyramidal subpopulations) to a delta-function-input or impulse. The simulations of this section can be regarded modelling event-related responses to events of short duration, in the absence of spontaneous activity or stochastic input. In the next section we will use more realistic inputs that comprise both stimulus-related and stochastic components.

The effects of inputs

Inputs u act directly on the spiny stellate neurons of layer 4. Their influence is mediated by the forward connections parameterised by the matrix C. When these connections are sufficiently strong, the output of the spiny stellate subpopulation saturates, due to the nonlinear sigmoid function in (2). This nonlinearly has important consequences for eventrelated responses and the ensuing dynamics. In brief, the form of the impulse response function changes qualitatively with input strength. To illustrate this point, we modelled a single area, which received an impulse at time zero and calculated the corresponding response for different values of c (Figure 4). With weak inputs, the response is linear, leading to a linear relationship between c and peak MEG/EEG responses. However, with large values of c neuronal activity leaves the linear domain of the sigmoid function in (2), the spiking saturates and the shape of the evoked response changes.

Figure 4 about here This behaviour is not surprising and simply reflects the nonlinear relationship between firing rates and postsynaptic depolarisation modelled by the nonlinearity. This nonlinearity causes saturation in the responses of units to intrinsic and extrinsic inputs. For example, when the input is strong enough to saturate spiny stellate spiking, the pyramidal response exhibits a short plateau (right panel in Figure 4). This saturation persists until the membrane potential of spiny stellate cells returns to its resting state. The sigmoid function models phenomenon; (i) at the single unit level, like refractoriness and spike rate adaptation and (ii) aspects of neuronal ensembles at the population level, like the distribution of thresholds involved in the generation of action potentials. The ensuing behaviour confers an inherent stability on dynamics because it is recapitulated in the response to all bottom-up influences, as shown next.

Bottom-up effects

The targets of the forward connections and extrinsic inputs are identical. Therefore, the effects of c and a F on event-related responses, are exactly the same. Figure 5 shows the simplest case of two areas (area 1 drives area 2). The difference, in relation to the previous configuration, is that area 1 has a gating effect. This is basically a low-pass filter, which leads to greater variation of the response in area 2, relative to responses elicited by direct input to area 2 (c.f. Figure 4) . For instance, the small negative response component in area 1, which follows the first positive deflection, is dramatically enhanced in area 2 for strong forward couplings. Again, this reflects the nonlinear behaviour of subpopulations responding to synaptic inputs.

Figure 5 about here

Level in hierarchy

As mentioned above, activity is subject to low-pass filtering, by synaptic processes, each time encounters a cortical region. A simple and intuitive consequence of this is that the form of event-related responses changes with each successive convolution in the hierarchy. To illustrate this point, we consider a feedforward configuration composed of five regions (Figure 6). We see in Figure 6 that, in addition to the propagation lag that delays the waveform at each level, the event-related response is more enduring and dispersed in higherlevel areas. A useful heuristic here is that late components of evoked responses may reflect hierarchical processing at a deeper level. This effect is independent of synaptic time constants and connectivity parameters.

Figure 6 about here This simple delay and dispersion is not necessarily seen with more realistic configurations that involve top-down effects. In this context, late response components in higher cortical areas can reenter [START_REF] Edelman | Neural Darwinism: selection and reentrant signaling in higher brain function[END_REF] lower levels engendering complicated and realistic impulse response functions. In the reminder of this section we look at the effects of adding backward and then lateral connections to the forward architecture considered above.

Top-down effects

Top-down connections mediate influences from high to low-level regions. Incoming sensory information is promulgated through the hierarchy via forward, and possibly lateral, connections to high-level areas. To demonstrate the effect of backward connections on MEG/EEG, we will consider a minimal configuration composed of two areas (Figure 7). The fact that the forward and backward connections are different renders this functionally asymmetric architecture hierarchical. Although asymmetric, the presence of forward and backward connections creates loops. This induces stability issues as shown in Figure 7:

when backward connections are made stronger, damped oscillations (a B =1; a B =10) are transformed into oscillations which ultimately stabilise (a B =50) because of the saturation described in the previous subsection. Therefore, with a B =50, the stable attractor is a limit cycle and the resting state point attractor looses it dynamic stability. The dependence of oscillations on layers, loops and propagation delays has been the subject of much study in computational models [START_REF] Lumer | Neural dynamics in a model of the thalamocortical system. I. Layers, loops and the emergence of fast synchronous rhythms[END_REF].

Figure 7 about here From a neurobiological perceptive, the most interesting behaviours are shown just prior to this phase-transition 1 when damped oscillations are evident. Note that the peaks of the evoked response, in this domain, occur every 100 milliseconds or so. This emulates the expression of late components seen empirically, such as the N300 or P400. The key point here is that late components, in the EEG/MEG, may reflect reentrant effects mediated by backward connections in hierarchical architectures. This observation fits comfortably with the notion that late EEG/MEG components reflect endogenous processing and depend explicitly on top-down effects. In short, late components may depend on backward connections and reflect a reentry of dynamics to hierarchically lower processing areas. This 1 A phase-transition refers to the qualitative change in the systems attractor caused by changes in the systems parameters, here the coupling parameters. In the present context, increasing the backward coupling causes the point attractor to loose its dynamic stability (stability under small perturbations) and the emergence of a limitcycle attractor. The nature of the phase-transition is usually assessed in terms of Lyapunov exponents (eigenvalues of the systems Jacobian

x f  

). When the system has a point attractor the imaginary part of the principal or largest exponent is zero. A limit cycle has non-zero imaginary parts and chaotic attractors have at dependency can be seen clearly by comparing the two left-hand panels in Figure 7 that show the emergence of late components on increasing the backward connection from one to ten.

The phase transition from damped late components to oscillations is critical. Before the transition the system is controllable. This means that the response can be determined analytically given the input. As discussed in (Friston, 2000a) long impulse responses endow the brain with a 'memory' of past inputs that enables perceptual processing of temporally extended events. In (Friston, 2000b) this was demonstrated using a Volterra kernel formulation and the simulation of spatiotemporal receptive fields in the visual system. However, after the transition it is no longer possible to determine when the input occurred given the output. This violates the principle of maximum information transfer [START_REF] Linsker | Perceptual neural organization: some approaches based on network models and information theory[END_REF] and precludes this sort of response in the brain. In short, it is likely that reentrant dynamics prolong neuronal transients but will stop short of incurring a phase-transition to oscillations. If this phase transition occurs it is likely to be short-lived or pathological (e.g. photosensitive seizure activity).

It should be noted that the oscillations in the right hand panels of Figure 7 do not represent a mechanism for induced oscillations. The oscillations here are deterministic component of the systems impulse response function and are time-locked to the stimulus. Induced oscillations, by definition, are not time-locked to the stimulus and probably arise from a stimulus-related change in the system's control parameters (i.e. short-term changes in connectivity). We will return to this point later.

Lateral connections

Lateral connections link different regions at the same level in the hierarchy. They can be unidirectional or bidirectional as shown for the model in Figure 8 with two areas. The main difference between forward and unidirectional lateral connections is that the latter target pyramidal cells. This means that the MEG/EEG signal is not so constrained by nonlinear saturation in layer 4 units. Therefore, as shown in Figure 8a, the event-related response does not saturate for strong lateral connectivity values a L . On the other hand, bilateral connections and hierarchical connections are differentiated by the fact that bilateral connections are completely symmetric, which enables them to create a synchronisation manifold [START_REF] Breakspear | Nonlinear interdependence in neural systems: motivation, theory, and relevance[END_REF][START_REF] Breakspear | Nonlinear phase desynchronization in human electroencephalographic data[END_REF]. A comparison of Figure 8b and Figure 7 shows that a special aspect of bilateral connections is their ability to support dynamics that are in phase. least one real positive exponent. We do not present a stability analysis or the Lyapunov exponents in this work, because the phase-transitions are self-evident in the trajectories of the system. This sort of zero-lag phase-synchronisation is commonplace in the brain. Its mediation by lateral connections in this model concurs with previous modelling studies of zero-lag coupling in triplets of cortical areas that involve at least one set of bilateral or reciprocal connections [START_REF] Chawla | Zero-lag synchronous dynamics in triplets of interconnected cortical areas[END_REF]. For very large values of a L , architectures with bilateral connections are highly nonlinear and eventually undergo a second phase transition (see Figure 8b).

Ongoing and event-related activity

So far, we have considered noise-free systems. Event-related responses were modelled in terms of deterministic impulse responses that were unique to a given neuronal configuration.

In this context it is not necessary to evoke the notion of averaging. However, real MEG/EEG signals show a great variability from trial to trial [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF]. In this section we model this variability by adding a stochastic component (a zero-mean Gaussian process) to the input u. The output corresponding to a stochastic input is referred to ongoing activity, i.e. oscillations in the MEG/EEG signal that share no phase relationship with the stimulus. This does not mean that ongoing activity should be considered as noise, just that there is no obvious stimulus-related structure.

Ongoing activity

Ongoing activity is shaped by the same nonlinear convolution experienced by deterministic inputs. In the context of stationary inputs, the outputs can be characterised in terms of their spectral properties, which are determined by the generalised transfer functions of the Volterra kernels associated with any controllable analytic system. The impulse response function is the first-order kernel. As soon as the connectivity parameters of a hierarchical network change, the principal modes of this network, defined by the principal frequencies of oscillations, are modulated [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF]. As an illustration, let us consider the simple hierarchical model of two cortical areas established in the previous section (Figure 9) with two configurations, which differ in the strength of backward connections (a B =1 and a B =10). The corresponding frequency spectra, of pyramidal cell depolarisation of the two areas, show that the change in connectivity induces a profound modulation of their spectral profile. As one might intuit, increasing backward connections induce a peak at the same frequency of the damped oscillations in the impulse response function. This is an important aspect of ongoing activity in the sense that its spectral behaviour may be very close to that of evoked transients as shown in [START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF].

Induced vs. evoked dynamics

This modulation of oscillatory dynamics, by the systems coupling parameters, provides a natural model for event-related changes in rhythmic activity. This phenomenon is known as event-related synchronisation (ERS) in frequency bands showing an evoked increase in power, or conversely, event-related desynchronisation (ERD) for decreases [START_REF] Basar | EEG-brain dynamics: Relation between EEG and brain evoked potentials[END_REF][START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. In light of the above connectivity-dependent changes in power, ERD and ERS may reflect the dynamics induced by evoked changes in short-term plasticity. The key difference between evoked and induced transients relates to the presence or absence of changes in the systems control parameters, here connectivity or synaptic efficacy. Evoked changes are not necessarily associated with parameter changes and any complicated response can be ascribed to transients that arise as the systems trajectory returns to its attractor. Conversely, induced responses arise from perturbation of the attractor manifold itself, by changes in the parameters and ensuing changes in the dynamics. This distinction was discussed in [START_REF] Friston | Another neural code?[END_REF] in relation to MEG dynamics and modelled using asymmetric connections between two areas in (Friston, 2000a).

Empirically, the ERS/ERD approach is used to look for MEG/EEG power changes of rhythmic activity induced by external events. This phenomenon has been modelled, in the case of alpha rhythms, by a computational model of thalamocortical networks [START_REF] Suffczynski | Computational model of thalamo-cortical networks: dynamical control of alpha rhythms in relation to focal attention[END_REF]. It has been shown that a key mechanism is the modulation of functional interaction between populations of thalamocortical cells and the reticular nucleus. This and related issues will be addressed in a subsequent paper on induced responses. Here we focus on evoked changes.

The present analysis concludes by examining the sensitivity of evoked transients to state changes caused by ongoing activity, when the parameters are held constant. In particular we see if this sensitivity can be expressed as phase-resetting. The key issue here is the presence of nonlinear interactions between the current state of the system and its perturbation by a stimulus.

Phase-resetting and nonlinear interactions.

In the following, we investigate the effect of ongoing activity, on stimulus-dependent responses, to reconcile apparently contradictory conclusions from studies of event-related potentials. On one hand, classical studies have shown that event-related potentials are associated with amplitude changes in the MEG/EEG signal that represent a linear summation of an impulse response and ongoing activity [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF][START_REF] Shah | Neural dynamics and the fundamental mechanisms of event-related brain potentials[END_REF]. In this scheme, the variability at the single-trial level is due to, and only to, ongoing activity, which is removed after averaging to estimate the impulse response. On the other hand, it has been hypothesised that event-related waveforms, obtained after averaging, could be due to a phase-resetting of ongoing activity with no necessary change in the amplitude (i.e. power) of any stimulus locked transient [START_REF] Jansen | Phase synchronization of the ongoing EEG and auditory EP generation[END_REF][START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF]. Although mathematically well defined, the neural mechanisms that could instantiate phase-resetting of ongoing activity are unknown.

We will take phase-resetting to imply a nonlinear interaction between ongoing activity and stimulus-related input that results in phase-locking to stimulus onset. Although phaselocking can be produced by evoking oscillatory transients (i.e. amplitude modulation) this mechanism involves no change or resetting of the ongoing dynamics. To assess the contribution of phase-resetting in our simulations we therefore need to look for interactions between ongoing and stimulus-related inputs that produced phase-locking in the outputs. As mentioned in the introduction this was addressed, in a simple way, by subtracting the response to ongoing activity alone from the response to a mixture of ongoing activity and stimulus input. In the absence of interactions this difference (the evoked response) should be the same. On the other hand, if interactions are prevalent, the difference should change with each realisation of ongoing activity. We performed these analyses with different levels of input and assessed the degree of phase-locking in the outputs with the phase-locking value (PLV) [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF][START_REF] Tallon-Baudry | Stimulus specificity of phaselocked and non-phase-locked 40 Hz visual responses in human[END_REF]:

trials t j t PLV )) ( exp( ) (  
where the instantaneous phase ) (t  was obtained from the Hilbert transform [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF].

Testing for interactions

To evaluate the effect of background activity on single-trial event-related responses, we used the two area hierarchical model above, with a B =1 (Figure 9). The first area was driven by an impulse function (stimulus) and Gaussian random noise (background activity) of standard deviation =0.05. The output of this region can be considered a mixture of evoked response and ongoing activity. We considered two conditions: one with low levels of mixed input (c=10 2 ) and another with high levels (c=2.10 4 ). These values were chosen to emphasise the systems nonlinear properties; with the smaller value of c, neuronal responses remain largely to the linear regime of the nonlinear function. The larger value of c was chosen so that excursions of the states encroached on the nonlinear regime, to produce neuronal saturation in some trials. In both cases, the stimulus was a delta-function. The simulated responses, for 100 trials, are shown in Figure 10.

When input levels are low (left hand side of Figure 10), event-related activity at the singletrial level shows a relatively reproducible waveform after stimulus onset (Figure 10B). This transient is reflected in the ERP/ERF after averaging (Figure 10C). To confirm the experimental results of [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF], we decomposed each event-related response into two components. First, the stochastic component (the response to ongoing activity alone -Figure 10D) and second, an extra component elicited by adding the stimulus (Figure 10E). This is the difference between the response elicited by stochastic component alone (Figure 10D) and the response to the mixed input (Figure 10B). If the system was linear, these differences should not exhibit any variability over trials, and thus define the "reproducible response" [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF]. Effectively, the stimulus-dependent component shows no variability and we can conclude that the response components due to stimulus and ongoing activity are linearly separable. In other words, there are no interactions that could mediate phase-resetting. Despite this, there is ample evidence for phase-locking. This is shown in Figure 10F, using the PLV index.

However, the situation is very different when we repeat the simulations with high input levels (right hand side of Figure 10). In this context, the event-related responses do not show any obvious increase in amplitude after the stimulus (Figure 10B). However, the averaged event-related activity (Figure 10C) is very similar to that above (left hand side of Figure 10C). The fact that one obtains an ERP by averaging in this way suggests that the stimulus input induced phase-resetting of the ongoing oscillations. This is confirmed by the large variation in stimulus-dependent components from trial to trial. This variation reflects nonlinear interactions between the stimulus and ongoing activity (Figure 10E). These interactions are associated with phase-locking as shown in Figure 10F.

Figure 10 about here In summary, the fact that the difference in evoked responses with and without background noise (panel E, Figure 10) shows so much variability, suggests that background activity interacts with the stimulus: When ongoing activity is high, stellate cells outputs saturate and the stimulus-related response is attenuated. Conversely, when ongoing activity is low the evoked-response is expressed fully. This dependency on ongoing activity is revealed by variation in the evoked responses with high input levels. In conclusion, the apparently contradictory results presented in [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF][START_REF] Jansen | Phase synchronization of the ongoing EEG and auditory EP generation[END_REF][START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF][START_REF] Shah | Neural dynamics and the fundamental mechanisms of event-related brain potentials[END_REF]) can be reproduced in most part and reconciled within the same framework. With high activity levels, the ongoing and stimulus-dependent components interact, through nonlinearities in the population dynamics, to produce phase-resetting and a classical ERP on averaging. When activity is lower, the stimulus and endogenous dynamics do not interact and the ERP simply reflects the transient evoked by stimuli that is linearly separable from ongoing dynamics.

Discussion

We have shown that it is possible to construct hierarchical models for MEG/EEG signals. To that end, we have assumed an architecture for cortical regions and their connections. In particular, we have used the Jansen model [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] for each source, and a simplified version of the connection rules of [START_REF] Felleman | Distributed hierarchical processing in the primate cerebral cortex[END_REF] to couple these sources. Here, we have fixed the parameters intrinsic to each source (synaptic time constants, output function, and intrinsic connections) and have focused on manipulating extrinsic connections among areas. We have used this model to address some issues in the genesis of evoked responses as observed with MEG/EEG.

Inferring neural mass models and their parameters, on the basis of the EEG and MEG alone, is a dynamical inverse problem that does not have a unique solution. The model we have described is one possible model among many others. We have used this model to afford a balance between simple neuronal models that are unrealistic but sufficient to emulate dynamics seen empirically, and more complicated models that have an explicit relation to electrophysiology but are difficult to manipulate. In future work we will compute the relative likelihood of, or evidence for, different models (with Bayesian model selection) using real MEG/EEG data (David et al., 2004b). At this stage we focus on establishing the face validity of a representative model by showing it can reproduce a range of biological behaviours.

Deterministic simulations

In section 3, we ignored the stochastic component of MEG/EEG signals. This allowed us to study various configurations in terms of their deterministic input-output behaviour. Each type of extrinsic connection (forward, backward and lateral) had specific effects on eventrelated dynamics. The key conclusions from these analyses were; (i) When forward connections, mediating bottom-up or extrinsic inputs, are sufficiently strong, nonlinear mechanisms cause a saturation of neuronal responses. This endows the system with an inherent stability that precludes non-dissipative dynamics. (ii) The duration of evoked transients increases with the hierarchical depth or level of processing. (iii) When backward or bilateral connections are added, evoked transients become more protracted, exhibiting damped oscillations. These are formally identical to late or endogenous components. These simulations suggest that late components are mediated by reentrant dynamics within cortical hierarchies. Increasing the strength of extrinsic backward or bilateral connections causes the number of late components to increase until a phase transition into a quasioscillatory regime, whose dynamic stability is assured by the nonlinear saturation above. We discussed, briefly, the importance of this phase transition in relation to information theoretic constraints on neuronal computations. In short, it is likely that selective pressure will extend the duration of transients so that the mutual information between current activity and past perceptual brain states is maximised (see Friston, 2000a for a fuller discussion). However, phase-transitions to oscillation are not adaptive because the system is no longer controllable and the informational link with the past is lost. A simple mechanism for increasing the memory of neuronal systems is to re-enter past states from hierarchical levels via backward connections. This is precisely what we simulated. In summary, late components may represent dynamic "echoes" that allow current sensory information to interact with perpetual constructs from higher areas, based on previous sensory input. We have discussed one form of this hierarchical reentry, using an empirical Bayesian perspective on predictive coding [START_REF] Friston | Functional integration and inference in the brain[END_REF].

Stochastic simulations

In the second section, we introduced stochastic components to the inputs, to study eventrelated responses in the context of background activity. This work was motivated by apparently contradictory views of ERP/ERF generation. Some authors [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF][START_REF] Shah | Neural dynamics and the fundamental mechanisms of event-related brain potentials[END_REF] have argued that ERPs/ERFs can be considered as the impulse response function of a linear system, despite the huge variability at the single-trial level. In contrast, it has been proposed [START_REF] Duzel | A multivariate, spatiotemporal analysis of electromagnetic time-frequency data of recognition memory[END_REF][START_REF] Jansen | Phase synchronization of the ongoing EEG and auditory EP generation[END_REF][START_REF] Klimesch | Phaselocked alpha and theta oscillations generate the P1-N1 complex and are related to memory performance[END_REF][START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF] that ERPs/ERFs could be due to a phase-resetting process that does not require the notion of a linear impulse response. This induces some hypothetical neural mechanism that can implement a phase-resetting of ongoing activity. These two views are supported by analyses of real MEG/EEG data. We showed that our model reproduces data that support both views. In our framework, the key factor is the level of activity. When activity is low, the model operates in a quasi-linear regime and reproduces the linear behaviours reported in [START_REF] Arieli | Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses[END_REF]. Conversely, when activity is high, single-trial responses do not show any amplitude modulation and yet still produce an ERP on averaging. This is consistent with the observations of [START_REF] Jansen | Phase synchronization of the ongoing EEG and auditory EP generation[END_REF] and [START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF].

Our analyses suggest that phase-resetting involves a nonlinear interaction between stimulus-related responses and ongoing activity. This observation is very reminiscent of a similar finding, pertaining to fast-oscillatory dynamics induced in simulated populations.

These simulations [START_REF] Chawla | Relating macroscopic measures of brain activity to fast, dynamic neuronal interactions[END_REF] showed that functional connectivity or dynamic integration between two populations increases with mean background activity and with stimulus-related rate modulation. Furthermore, as the background activity increases, the populations become increasingly sensitive to the intensity of the stimulus in terms of a predisposition to transient phase-locking. This reflected an interaction between background activity and stimulus-intensity in producing dynamic correlations. The mechanism of these interactions was modelled at the level of membrane time constants and temporal integration over milliseconds, using single-compartment units and the Hodgkin-Huxley formalism. Our simulations were at a much less detailed level, employing mean field-like approximations. Furthermore, they addressed evoked transients as opposed to induced fast dynamics.

However, they both speak to the central role of nonlinear interactions between ongoing activity and stimuli in causing phase-locking. In our case, of evoked responses, this phaselocking was to stimulus onset and the mechanism can be attributed directly to the only nonlinearity in our model: namely the sigmoid response function transforming depolarisation to firing rate. This nonlinearly renders ongoing activity sensitive to stimulus perturbations (this is the definition of an interaction). This sensitivity is due to the saturating nature of the sigmoid function. When ongoing activity is high the system states are close to the nonlinear regime of the sigmoid function and perturbations due to the stimulus cause neuronal saturation. This saturation causes the ongoing activity to "forget" its history and experience a resetting of its phase.

Phase resetting vs. Phase-locking

In this paper we have made a clear distinction between phase-locking and phase-resetting.

Phase-locking, as measured by the PLV implies a statistical dependency among phases, over ERPs (i.e. knowing the phase from one ERP allows one to predict the phase of another).

Phase-locking can be mediated in a number of ways. As shown in Figure 10F, there is no real difference between the phase-locking measures obtained from linear and nonlinear single-trial responses. Conversely, phase-resetting refers to a change in the phase, within ERPs. This can be mediated by linear or nonlinear mechanisms (i.e. a second or high-order interaction between ongoing dynamics and the stimulus). The key debate here is not about phase-resetting itself, but whether it is caused by nonlinear mechanisms (i.e. without changes in amplitude). We therefore focused on the distinction between linear or nonlinear, rather than phase-resetting per se. In a future paper we will look at phase-resetting, explicitly, using background activity that is periodic and of known phase (see below).

In this paper, phase-resetting is taken to imply a nonlinear interaction between ongoing dynamics and a stimulus that induces a phase-locking over trials. This generalizes the notion of alpha phase-resetting that is usually characterised empirically. Our simulations and analysis were not limited to a single frequency. The ongoing activity we used was broadspectrum and the phase-locking measure was not frequency-specific. In empirical studies it is not possible to analyse the interaction between stimulus and ongoing activity, because no single trial can be replicated in the absence of the stimulus. This means that empirical studies have to use epochs that show oscillations at a particular frequency to demonstrate phaseresetting. In a future communication we will apply this phase-resetting analysis (c.f. [START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF]) to stimulations in which our stochastic ongoing activity is replaced with a deterministic sinusoidal forcing term, in the alpha range. This will allow us to partition phase-locking into components that are mediated by linear superposition and phase-resetting respectively.

Evoked vs. induced oscillations

We have also noted that event-related activity is not restricted to ERPs/ERFs. It encompasses event-related changes in oscillations (power and synchronisation). We have suggested that this phenomenon can be explained by changes of connectivity due to shortterm neuronal plasticity (see (Friston, 2000a) for a related discussion on nonlinear connections). We will be pursuing this in a subsequent paper on induced responses. The same changes in connectivity are also responsible for the modulation of long-range synchrony (or any type of interdependence measure) between distant MEG/EEG oscillators [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF]David et al., 2004a).

Modelling and estimation

Generally, it was striking how different connections engendered such diverse event-related waveforms. This raises the issue of how models, of the sort described here, can be applied to real data. The long-term agenda of our modelling programme is to establish the validity of neuronal network models so that they can be used as forward models to explain real data.

The key advantage of this approach is that the parameters of the model, which are estimated, have a direct physiological interpretation. Several studies have shown that it is possible to estimate the parameters of observation models using real MEG/EEG data, in the case of both rhythmic activity [START_REF] Valdes | Nonlinear EEG analysis based on a neural mass model[END_REF] and ERPs [START_REF] Rennie | Unified neurophysical model of EEG spectra and evoked potentials[END_REF]. An important issue, in this context, is the ability to constrain, or regularise, the inverse problem that dynamic forward models like this pose. These issues can be addressed using Bayesian inference and dynamic causal modelling [START_REF] Friston | Dynamic causal modelling[END_REF]. A great advantage of using neuronally plausible forward models is that constraints on the solution can be applied, in an informed way, using priors on the model parameters. At present, our team has published Bayesian inference procedures for dynamic causal modelling of functional MRI [START_REF] Friston | Dynamic causal modelling[END_REF][START_REF] Penny | Comparing dynamic causal models[END_REF]. In forthcoming reports, we will extend this approach to a variety of models, including the neural mass model for MEG/EEG described here. A preliminary report has already been published in which we estimate the parameters of this neural mass model but without localising the sources (David et al., 2004b). In our next communication we will combine neural mass models and forward MEG/EEG modelling to study cognitive functions, and in particular to infer changes in connectivity among experimental conditions.

Conclusion

We have shown that neural mass models [START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF][START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF][START_REF] Lopes Da Silva | Alpha rhythms: noise, dynamics and models[END_REF][START_REF] Nunez | The brain wave equation: A model for the EEG[END_REF][START_REF] Rennie | Unified neurophysical model of EEG spectra and evoked potentials[END_REF][START_REF] Robinson | Prediction of electroencephalographic spectra from neurophysiology[END_REF][START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity?[END_REF][START_REF] Suffczynski | Computational model of thalamo-cortical networks: dynamical control of alpha rhythms in relation to focal attention[END_REF][START_REF] Valdes | Nonlinear EEG analysis based on a neural mass model[END_REF][START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF] can reproduce a large variety of MEG/EEG signal characteristics. The potential advantage they afford, in comparison to standard data analysis, is their ability to pinpoint specific neuronal mechanisms underlying normal or pathological activity. Effort is needed to incorporate them, more systematically, in MEG/EEG analyses to enable enquiry into mechanistic questions about macroscopic neuronal processes. In forthcoming studies, we will describe the estimation of such models in a Bayesian framework [START_REF] Friston | Classical and Bayesian inference in neuroimaging: theory[END_REF] and will apply the present model to the analysis of real ERPs.

                                                  2 ) ( 7 ) ( 8 3 1 ) ( ) ( 1 ) ( ) ( ) ( 8 ) ( 8 ) ( 7 2 ) ( 3 ) ( 6 ) ( 7 4 ) ( 6 ) ( 6 ) ( 3 2 ) ( 2 ) ( 5 ) ( 1 2 1 ) ( ) ( 1 ) ( ) ( ) ( 5 ) ( 5 ) ( 2 2 ) ( 1 ) ( 4 1 1 ) ( ) ( ) ( 1 ) ( ) ( ) ( 4 ) ( 4 ) ( 1 2 , , 2
j i ij j L l i j j ij j F e e i i i x x y S t y S j i A t y S j i A H x x x x x x S H x x x x x x S t y S j i A t y S j i A H x x x x x y S u i C t y S j i A t y S j i A H x x x                                                                                                            (A.2) where ) ( 3 ) ( 2 ) ( i i i x x y  
is the pyramidal cells membrane potential (depolarisation) of region i. This is assumed to be proportional to the cortically reconstructed current source densities obtained from MEG/EEG scalp signals. This equation was integrated using standard Runge-Kutta techniques [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] with an integration time-step of 1 ms. Apart from the extrinsic connectivity parameters, we used the same parameters for every simulation: H e =3.25, H i =29.3,  e =10ms,  i =15ms, = (ij) =10ms,  1 =50,  2 =40,  3 = 4 =12, e 0 =2.5, v 0 =0, r=0.56.

Figure Captions

Figure 1: Connection rules adopted for the construction of hierarchical models for MEG/EEG signals. These rules are a simplified version of those proposed by [START_REF] Felleman | Distributed hierarchical processing in the primate cerebral cortex[END_REF]. The cortical sheet is divided into two components: the granular layer (layer 4) and the agranular layers (supra-and infra-granular layers). Bottom-up connections originate in agranular layers and terminate in layer 4. Top-down connections only engage agranular layers. Lateral connections originate in agranular layers and target all layers. differ from those elicited by forward connections (Figure 4). In particular, the saturation of layer 4 is not so important and the signal exhibits less saturation for large a L . b) Increasing bi-directional lateral connections has a similar effect to increasing backward connections.

The main difference is the relative phase of evoked oscillations, which are synchronised at zero-lag. For very large values of a L , the model is highly nonlinear and eventually exits the oscillatory domain of parameter space.

  Figure 3 about here

Figure 8

 8 Figure 8 about here

Figure

  Figure 9 about here

Figure 2 :

 2 Figure 2: Jansen's model of a cortical area. Three neuronal subpopulations are considered to model a cortical area. Pyramidal cells interact with both excitatory and inhibitory interneurons with the connectivity constants 1 4 3 1 2 25 . 0 , 8 . 0       . The parameters H e,i

Figure 3 :

 3 Figure 3: Hierarchical connections among Jansen units (Figure 2) based on simplified Felleman and van Essen rules (Figure 1). Long range connectivity is mediated by pyramidal cells axons. Their targets depend upon the type of connections. Coupling or connectivity parameters control the strength of each type of connection: a F for forward, a B for backward, and a L for lateral.

Figure 4 :

 4 Figure 4:The strength of input modulates the shape of MEG/EEG signal. The output of one area (variable y, see Appendix) has been calculated for different values of c, the strength of forward connections mediating input u (delta function). When c is small (c=1, c=1000), the output is not saturated and the MEG/EEG signal (y) amplitude is linearly related to c. For large values of c (c=10 6 , c=10 9 ), spiny stellate cells saturate and the shape of event-related MEG/EEG response changes substantially.

Figure 5 :

 5 Figure 5: The MEG/EEG signal of area 1 (black) and area 2 (grey) is plotted as a function of the forward connectivity a F . Bottom-up connectivity has the same effect as input connectivity

Figure 6 :

 6 Figure 6: A feedforward system composed of five areas. The MEG/EEG signal (y) of each area elicited by a single pulse on area 1 is plotted in successive panels from left to right. Event-related activity lasts longer in high-level cortical areas of feedforward architectures. At each level in the hierarchy, the event-related response of pyramidal cells experiences successive low-pass filters, embodied by synaptic processes that transform the input signals to output.

Figure 7 :

 7 Figure 7: Backward connections have a key influence on the stability of MEG/EEG eventrelated activity as demonstrated by this simple model composed of two areas (area coded in black and area 2 coded in grey). The forward connectivity a F has been fixed to 40 and backward connectivity a B varies between 1 and 50 from left to right. When top-down effects are small, their reentry leads to longer lasting event-related responses characterised bydamped oscillations (a B =1; a B =10). However, over a critical threshold of a B (which depends upon a F ), the system undergoes a phase transition, loses its point attractor and expresses oscillatory dynamics (a B =25; a B =50).

Figure 8 :

 8 Figure 8: The effects of lateral connections are shown with a simple model composed of two areas (black: area 1, grey: area 2). The depolarisation of pyramidal cells (y) is plotted for several values of lateral connectivity a L . a) Unilateral connections support transients that

Figure 9 :

 9 Figure 9:The modulation of backward connectivity (a B =1 or a B =10) has huge effect on the power spectrum of ongoing MEG/EEG dynamics (variable y is plotted in black for area 1 and in grey for area 2). When a B increases from 1 to 10, there is loss of power below 3 Hz, and an excess between 3 and 7 Hz. The amplitude spectra in the right panel were obtained by

Figure 10 :

 10 Figure 10: Event-related responses in the context of ongoing activity (100 trials). Two hierarchically connected regions are considered. Input (ongoing and stimulus-related) enters into the system through region 1. Two levels of input are considered: weak on the left hand side (c=10 2 ), strong on the right hand side (c=2.10 4 ). The successive horizontal panels show different types of activity. The time scale is identical for each panel and shown at the bottom. A) Inputs, comprising a delta function and Gaussian noise of standard deviation 0.05 (stimulus onset at time 0). B) Event-related activity (y) at the single-trial level. The time series over trials is shown (area 1 is above area 2). C) Averaged event-related response estimated by averaging over epochs shown in B (area 1 in black, area 2 in grey). D) Responses (y) to the noisy input without the delta function, shown in the same format as in B. E). Stimulus-dependent component obtained from subtracting D from B. F). Phaselocking value computed form time series in B, which exhibits a transient phasesynchronisation to peristimulus time (Black: area 1, grey: area 2).

Acknowledgments

This work was supported by the Wellcome Trust.

Appendix

In this appendix, we provide differential equations for the i th of l cortical areas of the hierarchical Jansen model as described in section 2.3. The activity of the i th cortical area is described by eight state variables x (i) , the rate of change of which is described by an equation of the following form:

u (i) is the extrinsic input to region i not modeled by other areas. S(y (j) ) is the firing rate of pyramidal cells of region j. A F , A B , A L are l x l connectivity matrices for forward, backward and lateral connections, respectively. C is the connectivity vector for extrinsic inputs.  (ij) is the propagation delay between region i and region j. Equations ( 1) and (2) in the main text describe neuronal input and output operators that correspond to the following state equation which is the specific version of A.1 for our model:


domain: math info scco

text:
Roberto Rinaldi

Marco Ripà
email: marco.ripa@mensa.it

Optimal cycles enclosing all the nodes of a k-dimensional hypercube

Keywords: Euclidean space, Optimization, Link distance, Combinatorics. 2020 Mathematics Subject Classification: 05C38 (Primary); 05C12, 91A43 (Secondary)

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Given k ∈ N -{0, 1}, let the finite subset of 2 k points belonging to the Euclidean space E ⊂ R k be defined as H(2, k) := {{0, 1} × {0, 1} × • • • × {0, 1}} ⊂ R k , where the symbol "×" denotes the well-known cartesian product. Then, the problem of joining all the nodes of H(2, k) with a connected set of segments having minimum cardinality is equivalent to asking to ourselves which is the minimum-link covering tree embedding all the nodes of a k-dimensional hypercube. This question is not an open problem, since Dumitrescu and Tóth (see [START_REF] Dumitrescu | Covering grids by trees[END_REF], Figure 2), in 2014, easily showed that it is sufficient to connect all the 2 k-1 pairs of opposite nodes with as many segments so that all of them (i.e., exactly 2 k-1 line segments) meet in the center, C ≡ 1 2 , 1 2 , . . . , 1 2 . Now, it is obvious to understand why the number of (line) segments of the above mentioned minumum-link covering tree for H(2, k) cannot match the link-length of any polygonal chain covering the same set of 2 k nodes.

Thus, we are interested in solving a multidimensional thinking outside the box problem, quite similar to the k-dimensional generalization of the infamous nine dots puzzle [START_REF] Chein | Working memory and insight in the nine-dot problem[END_REF][START_REF] Kershaw | Training for insight: The case of the nine-dot problem[END_REF][START_REF] Loyd | Cyclopedia of Puzzles[END_REF] which was definitely solved in 2020 by Ripà (see [START_REF] Ripà | Solving the 106 years old 3 k points problem with the clockwise-algorithm[END_REF]): the crucial special case of finding a minimum-link polygonal chain covering any given H(2, k) := {{0, 1}×{0, 1}ו • •×{0, 1}} ⊂ R k , a fascinating challenge belonging to the general problem of finding a minimum-link covering trail for every set H(n, k) (see [START_REF]Sequence A318165 in the On-line Encyclopedia of Integer Sequences[END_REF][START_REF] Ripà | General conjecture on the optimal covering trails for any k-dimensional cubic lattice N k[END_REF][START_REF] Ripà | Minimum-link covering trails for any hypercubic lattice[END_REF]).

In order to introduce the original results of the present paper (Section 2), let us give a few of definitions first. Definition 1.1. Let P(m) := (S 1 )-(S 2 )-. . . -(S m+1 ), a polygonal chain consisting of m links, be well-defined through the sequence of its m + 1 vertices, so

P(m) ≡ {S 1 S 2 ∪ S 2 S 3 ∪ • • • ∪ S m S m+1 }.
In particular, for any d ∈ Z + : d ≤ m+1, let the d-th vertex of P(m) ⊂ R k be univocally identified by the k-tuple (x 1 , x 2 , . . . , x k ) (e.g., given P(3) := (0, 0)-(1, 0)-(1, 1)-(0, 1) and d : d = 2, we have

x 1 (S d ) = 1 ∧ x 2 (S d ) = 0, since (x 1 (S 2 ), x 2 (S 2 )) ≡ (1, 0) for the aforementioned minimum-link polygonal chain covering H(2, 2) := {{0, 1} × {0, 1}}). Definition 1.2. Accordingly to Definition 1.1, let h(2, k) denote the link-length of the minimum-link polygonal chain P(h(2, k)) := (S 1 )-(S 2 )-. . . -(S h(2,k)+1 ) visiting all the nodes H(2, k) := {{0, 1} × {0, 1} × • • • × {0, 1}} of the k-dimensional hypercube {[0, 1] × [0, 1] × • • • × [0, 1]}.
Definition 1.3. Let P(h(2, k)) be a (possibly self-intersecting) path if there is not any element of the set H(2, k) which belongs to more than one link of P(h(2, k)). Then, let P(h(2, k)) be a cycle if it is a path such that (S 1 ) ≡ (S m+1 ) (i.e., we call a cycle any closed path). Furthermore, we define as "perfect covering cycle", C(h(2, k)), any closed path such that no element of the set {S j | j = 1, 2, . . . , m + 1} belongs to more than two links of the given covering path for H(2, k).

Lastly, for clarity sake, let us specify that we will use vertices and links when we are referring to the turning points (usually we consider Steiner points which do not belong to

{[0, 1] × [0, 1] × • • • × [0, 1]} ⊂ R k ) of
the polygonal chains that we are taking into account, whereas we will prefer nodes and edges for the respective subsets of points entirely belonging to given k-dimensional hypercube.

Since the aim of this paper is to find polygonal chains, embedding H(2, k), which are optimal with respect to the number of line segments, we immediately point out that any covering cycle for H(2, k) is a covering path for the same set of points, and it is also a polygonal chain covering all the 2 k nodes of the hypercube

{[0, 1] × [0, 1] × • • • × [0, 1]} ⊂ R k .
Now, a constructive proof of the existence of covering cycles for H(2, 2) and H(2, 3) having link-length of only h(2, 2) = 3 and h(2, 3) = 6, respectively, has already been shown in Reference [START_REF] Ripà | General uncrossing covering paths inside the axis-aligned bounding box[END_REF] (e.g., the aforementioned paper, see pages 163-164 and Figures 7 to 9, provides an optimal upper bound for any k ∈ {2, 3} that is also achievable by taking into account only covering cycles, instead of generic polygonal chains, since H(2, 2)

⊂ P(3) = 1 2 , 3 2 -(2, 0)-(-1, 0)-1 2 , 3 2 and H(2, 3) ⊂ P(6) = 1 2 , 1 2 , 3 2 -(2, 2, 0)-(-1, -1, 0)-1 2 , 1 2 , 3 2 -(2, -1, 0)-(-1, 2, 0)-1 2 , 1 2 , 3 2
). Thus, from here on, let us assume that k ≥ 2 is given. The goal of the present research paper is to show that h l (2, k) = 3 • 2 k-2 is a valid lower bound for any polygonal chain visiting all the elements H(2, k) and to constructively prove that the aforementioned lower bound is equal to the upper bound h u (2, k) = 3 • 2 k-2 [START_REF]Sequence A007283 in the On-line Encyclopedia of Integer Sequences[END_REF] which returns the minimum link-lenght of any perfect covering cycle for the same set of nodes.

It follows that we are going to prove that h u (2, k) = h l (2, k) by providing optimal covering cycles consisting of 3 • 2 k-2 links for any k ∈ N -{0, 1}, and this result will be shown in the next section.

Main Result

In order to prove that h l (2, k) = 3 • 2 k-2 holds for any k ∈ N -{0, 1}, we will introduce the following lemma.

Lemma 2.1. For any k ∈ N -{0, 1}, it is not possible to visit more than 4 distinct elements of the set

H(2, k) := {{0, 1} × {0, 1} × • • • × {0, 1}} ⊂ R k
by using a polygonal chain with 3 links.

Proof. We prove Lemma 2.1 by studying the generic trail P(3) ≡ {S 1 S 2 ∪ S 2 S 3 ∪ S 3 S 4 } that passes through 4 (distinct) nodes of H(2, k). Then, we will show that there is not a choice for these four nodes (and also for the considered Steiner points) that implies the existence of a fifth node belonging to P(3). We will start by demonstrating that there is not a trail P(2) that, passing through at least 3 nodes, visits a fourth node. Considering that we need to pass through at least 3 nodes, and being able to visit at most 2 nodes with the first segment, we can impose that S 1 S 2 passes through 2 nodes. Although, for convenience, we will impose S 1 ≡ V 1 ≡ O, the obtained result can be extended to any choice of S 1 .

Let S j be the origin of a given half-line q j . Given a parameter t j ∈ R : t j ≥ 0, the corresponding parametric equation is of the form q j = S j + t j • # » S j V j+1 , where V j indicates the i-th node of H(2, k) visited by P(m). In this way, for t j = 1, each of the Cartesian coordinates of the nodes must assume the value 0 or 1. Our goal is to show that this happens only for t j = 1, and if this occurs for other values of t j , then we want to show that the visited node is a node already visited previously.

Since the Steiner point S j+1 belongs to the considered half-line, let us denote by tj the value of the parameter t j such that S j+1 = S j + tj • # » S j V j+1 (i.e., S j+1 := S j+1 ( tj )).

The generic point S 2 is obtained as the last endpoint of a segment passing through V 2 ,

S 2 = S 1 + t1 • # » S 1 V 2 . (1)
Now, from here on, we will indicate the i-th coordinate of the generic point P, belonging to the Euclidean space R k , as x i (P) (e.g., P ≡ (x 1 (P), x 2 (P), . . . , x k (P))).

Consequently, from Equation (1), it follows that

x i (S 2 ) = t1 • x i (V 2 ) , (2)
where t1 ≥ 1.

Similarly, we will have that

S 3 = S 2 + t2 • # » S 2 V 3 . (3)
Hence,

x i (S 3 ) = t1 • x i (V 2 ) + t2 • (x i (V 3 ) -t1 • x i • (V 2 )) = t1 • x i (V 2 ) • (1 -t2 ) + t2 • x i (V 3 ) , (4)
where t2 ≥ 1.

Let us consider the segment S 2 S 3 . By disregarding the node V 3 , obtained by imposing t 2 = 1, we verify that it does not exist any node V j of H(2, k), belonging to S 2 S 3 , which has not been previously visited.

Thus, for i < k, we need to study all the x i (V j ) equations, showing that there are no solutions such that 0 < t 2 < 1. It is not necessary to continue beyond point V 3 , studying solutions for t 2 > 1.

If we encounter a node of the set H(2, k) after point V 3 , then we will necessarily have to study also the case in which point V 3 represents the furthest node and, with t 2 < 1, we will find the point studied previously.

As a result, we have to study the above mentioned x i (V j ) equations,

       t1 • x i (V 2 ) + t 2 • (x i (V 3 ) -t1 • x i (V 2 )) = x i (V j ) t1 > 1 . 0 < t 2 < 1 (5)
Hence,

x i (V 2 ) = x i (V 3 ) = 0 ⇒ t 2 ∈ R : 0 < t 2 < 1, (6)
x i (V 3 ) = 0 ⇒ t 2 = t1 -1 t1 . (7)
Consequently, all the solutions imply

x i (V 3 ) = 0. If x i (V 3 ) = 0 holds for all i : i < k, then (V 3 ) = (V 1 ).
Now, we are finally ready to study the generic trail

P(3) ≡ {S 1 S 2 ∪ S 2 S 3 ∪ S 3 S 4 }.
Thanks to the results discussed above, we know that if S 2 S 3 visits two nodes, then S 1 S 2 and S 3 S 4 visit one node each, so we can impose (from the beginning) that S 1 S 2 visits two nodes of H(2, k).

Such trail, P(3), can be built starting from the just described trail P(2), by simply adding a fourth generic Steiner point whose coordinates satisfy

S 4 = S 3 + t3 • # » S 3 V 4 , (8)
so we have

x i (S 4 ) = t1 • x i (V 2 ) • (1 -t2 ) + t2 • x i (V 3 ) + t3 • (x i (V 4 ) -( t1 • x i (V 2 ) • (1 -t2 ) + t2 • x i (V 3 ))) =( t1 • x i (V 2 ) • (1 -t2 ) + t2 • x i (V 3 )) • (1 -t3 ) + t3 • x i (V 4 ), (9)
with t3 = 1.

Before moving on the segment S 3 S 4 , let us make some considerations for a better understanding of the nature of the next step of the present proof.

We have that t1 ≥ 1 and t2 > 1, since t2 = 1 would imply S 3 = V 3 . It follows that S 3 S 4 would visit V 2 and V 3 , whereas it could not visit other nodes, since the set H(2, k) has not more than 2 collinear nodes.

Under these constraints, the following results are obtained so that we can use them to find the solutions of Equation [START_REF] Ripà | General conjecture on the optimal covering trails for any k-dimensional cubic lattice N k[END_REF].

1. t1 • (1 -t2 ) < 0, since t1 > 0 and (1 -t2 ) < 0. 2. t1 • (1 -t2 ) + t2 < 1, since t1 = 1 implies t1 • (1 -t2 ) + t2 = 1 and ∂ ∂ t1 ( t1 • (1 -t2 ) + t2 ) < 0 ∀ t1 > 1, t2 > 1.
Now, we consider the segment S 3 S 4 . By disregarding the node V 4 , obtained by imposing t 3 = 1, we verify that it does not exist any unvisited node V j of H(2, k), belonging to S 3 S 4 .

Thus,

             t1 • x i (V 2 ) • (1 -t2 ) + t2 • x i (V 3 ) + t 3 • (x i (V 4 ) -( t1 • x i (V 2 ) • (1 -t2 ) + t2 • x i (V 3 ))) = x i (V j ) t1 > 1 t2 > 1 0 < t 3 < 1 .
(10) Hence,

x i (V 2 ) = x i (V 3 ) = x i (V 4 ) = x i (V j ) = 0 ⇒ t 3 ∈ R : 0 < t 3 < 1, (11)
x i (V 2 ) = x i (V 3 ) = 1, x i (V 4 ) = x i (V j ) = 0 ∧ t1 = t2 t2 -1 ⇒ t 3 ∈ R : 0 < t 3 < 1, (12)
x i (V 2 ) = x i (V 4 ) = 0, x i (V 3 ) = x i (V j ) = 1 ⇒ t 3 = t2 -1 t2 , (13)
x i (V 2 ) = x i (V 3 ) = x i (V 4 ) = 1, x i (V j ) = 0 ⇒ t 3 = t1 (1 -t2 ) + t2 t1 (1 -t2 ) + t2 -1 . (14)
There cannot be two indices i, i ′ such that 12)&( 14)), since by imposing t1 = t2 t2 -1 we obtain t 3 = 0.

(x i (V 2 ) = x i (V 3 ) = 1, x i (V 4 ) = 0) ∧ (x ′ i (V 2 ) = x ′ i (V 3 ) = x ′ i (V 4 ) = 1, x ′ i (V j ) = 0) (see (
The uniqueness of the indices that simultaneously verify (11),( 12)&(13) implies that V 1 is visited twice, while the uniqueness of the indices that simultaneously verify (11),( 14)&(13) implies that V 2 is visited twice. Therefore, it is not possible to join more than 4 nodes of the given set with a polygonal chain consisting of only 3 links, and this concludes the proof of Lemma 2.1.

By invoking Lemma 2.1, we can easily prove the following theorem.

Theorem 2.2. Let k ∈ N -{0, 1} be given. The link-length of the covering trail

P(h(2, k)) for the set H(2, k) satisfies h(2, k) ≥ 3 • 2 k-2 .
Proof. There is a total of 2 k nodes to be visited. By Lemma 2.1, we can join a maximum of 4 nodes with a polygonal chain of link-length 3. Since H(2, k) has 2 k 4 groups of 4 nodes, the 4 nodes of each of these groups require a minimum of 3 segments to be visited. Therefore, h(2, k) ≥ 3 • 2 k-2 . Now, we need to find the shortest possible covering path, P(h(2, k)). At this purpose, it is possible to prove the following result.

Theorem 2.3. Given H(2, k) with k ∈ N -{0, 1}, it is always possible to construct a covering cycle P(h(2, k)) of link-length h(2, k) = 3 • 2 k-2 .
Proof. It is possible to create an algorithm that generates a covering circuit for H(2, k) whose link-length exactly coincides with the lower bound stated by Theorem 2.2. The algorithm is valid for any finite number of dimensions.

First of all, we notice that it is always possible to join the four nodes of a rectangle with a covering circuit of link-length 3. In fact, given the set {{0, 1} × {0, b}}, we can have the covering circuit (S 1 )-(S 2 )-(S 3 )-(S 1 ), with the elements of the set {S 1 , S 2 , S 3 } given by

S 1 ≡ b 2 , 3 2 ;
S 2 ≡ (-b, 0) S 1 + 1 3 • # » S 1 S 2 = (0, 1); S 3 ≡ (0, 2 • b) S 2 + 1 3 • # » S 2 S 3 = (0, 0), S 2 + 2 3 • # » S 2 S 3 = (b, 0); S 4 ≡ S 1 S 3 + 2 3 • # » S 3 S 4 = (b, 1).
We consider the sheaf of planes that has in common the line r := C + t • #» e k , where C ≡ 1 2 , 1 2 , . . . , 1 2 , 0 and #» e k := (0, 0, . . . , 0, 1) is the second vector of the canonical basis. These planes have parametric equation

C + t • #» e k + u • #» s with #» e k and #» s being linearly independent vectors. Let #» s l ∈ {{-1 2 } × {-1 2 , 1 2 } × {-1 2 , 1 2 } × • • • × {-1 2 , 1 2 } × {0}} ⊆ R k-2 be the vector such that l ∈ N : l < 2 k-2 . If t = u = 1, then we obtain the point V 4l+1 ∈ H(2, k) such that V 4l+1 ≡ (0, x 2 (V 4l+1 ), x 3 (V 4l+1 ), . . . , x k-1 (V 4l+1 ), 1), where    x i (V 4l+1 ) = 1 if x i ( #» s l ) = + 1 2 x i (V 4l+1 ) = 0 if x i ( #» s l ) = -1 2 . ( 15
) If t = 0 ∧ u = 1, then we obtain the point V 4l+2 ∈ H(2, k) such that V 4l+2 ≡ (0, x 2 (V 4l+2 ), x 3 (V 4l+2 ), . . . , x k-1 (V 4l+2 ), 1), where    x i (V 4l+2 ) = 1 if x i ( #» s l ) = + 1 2 x i (V 4l+2 ) = 0 if x i ( #» s l ) = -1 2 . ( 16
) If t = 0 ∧ u = -1, then we obtain the point V 4l+3 ∈ H(2, k) such that V 4l+3 ≡ (0, x 2 (V 4l+3 ), x 3 (V 4l+3 ), . . . , x k-1 (V 4l+3 ), 1), where    x i (V 4l+3 ) = 0 if x i ( #» s l ) = + 1 2 x i (V 4l+3 ) = 1 if x i ( #» s l ) = -1 2 . ( 17
) If t = 1 ∧ u = -1, then we obtain the point V 4l+4 ∈ H(2, k) such that V 4l+4 ≡ (0, x 2 (V 4l+4 ), x 3 (V 4l+4 ), . . . , x k-1 (V 4l+4 ), 1), where    x i (V 4l+4 ) = 0 if x i ( #» s l ) = + 1 2 x i (V 4l+4 ) = 1 if x i ( #» s l ) = -1 2 . ( 18
)
Being {σ} l := {σ l : l = 0, 1, 2, . . . , 2 k-2 -2, 2 k-2 -1} the set of the planes containing r and C + #» s l , in total we will have exactly 4 nodes of H(2, k) lying on each one of the aforementioned planes (i.e., there are 2 k-2 planes σ l such that 0 ≤ l < 2 k-2 , being 2 k-2 the multisubsets of size k -2 from the set { 1 2 , -1 2 }). Since each of the 2 k-2 multisubsets is different from all the others, it follows that it does not exists any pair of positive integers (j, j ′ ) : j ̸ = j ′ such that V j ≡ V j ′ . Consequently, each plane contains exactly 4 nodes that does not belong to any other plane σ l .

Thus, there are 2 k-2 planes that include 4 different nodes each, for a total of 2 k distinct nodes. Since a k-dimensional hypercube has exactly 2 k nodes, we conclude that each point lies on one and only one plane σ l .

Let l ∈ N 0 : 0 ≤ l ≤ 2 k-2 -1 be given. Then, the four points V 4l+1 ,V 4l+2 ,V 4l+3 , and V 4l+4 identify the nodes of a rectangle with base 1 and height 2 • (k -1). In fact, #» e 2 = (0, 1, 0, 0, . . . , 0) forms an orthogonal basis with vector #» s that has coordinate s 2 = 0.

We have already proven that the nodes in this orthogonal basis are in position (0, -1); (0, 1); (1, 1); (1, -1) so that, in this basis, they are vertices of a rectangle of base 1 and height 2.

Consequently, being 1 2 • (k -1) the magnitude of vector #» s , we get a height of 2• 1 2 • (k -1) = 2 • (k -1) in the canonical basis of the space E.

Finally, we have 2 k-2 rectangles, whose vertices can be covered by 2 k-2 covering circuits of link-length 3 (see Figures 1&2).

= 1 2 , 1 2 , 2 --1 2 , -1 2 , 0 -( 3 2 , 3 2 , 0)-1 2 , 1 2 , 2 -(-1 2 , 3 2 , 0)-( 3 2 , -1 2 , 0)-1 2 , 1 2
, 2 visits all the nodes of H(2, 3) once and only once (picture realized with GeoGebra [START_REF] Hohenwarter | [END_REF]). Lastly, using the described covering circuits, we get a circuit that starts and ends at a point that lies on the generating line of the sheaf of planes that contains all the planes belonging to {σ} l . Thus,

{V 4l+1 , V 4l+2 , V 4l+3 , V 4l+4 } ⊂ P m (4) ⇒ ∃! l ∈ {0, 1, . . . , 2 k-2 -1} : (S 3l+1 )-(S 3l+2 )-(S 3l+3 )-(S 3l+4 ) ⊂ {σ} l . (19
) As shown in Figure 3, we obtain a covering cycle for H(2, k) by the repetition, for every l ∈ {0, 1, 2, . . . , 2 k-2 -1}, of the covering circuit described by

S 3l+1 ≡ C + 3 2 • #» e 1 ; S 3l+2 ≡ C + 3 • #» s l S 3l+1 + 1 3 • # » S 3l+1 S 3l+2 ≡ V 4l+1 ; S 3l+3 ≡ C -3 • #» s l S 3l+2 + 1 3 • # » S 3l+2 S 3l+3 ≡ V 4l+2 , S 3l+2 + 2 3 • # » S 3l+2 S 3l+3 ≡ V 4l+3 ; S 3l+4 ≡ S 3l+1 S 3l+3 + 2 3 • # » S 3l+3 S 3l+4 ≡ V 4l+4 .
= 1 2 , 1 2 , 1 2 , 3 2 -(-1, -1, -1, 0)-(2, 2, 2, 0)-1 2 , 1 2 , 1 2 , 3 2 -(-1, -1, 2, 0)-(2, 2, -1, 0)- 1 2 , 1 2 , 1 2 , 3 2 -(-1, 2, -1, 0)-(2, -1, 2, 0)-1 2 , 1 2 , 1 2 , 3 2 -(-1, 2, 2, 0)-(2, -1, -1, 0)-1 2 , 1 2 , 1 2 , 3 2
joins all the nodes of H(2, 4) (picture realized with GeoGebra [START_REF] Hohenwarter | [END_REF]).

Therefore, we have constructively proven that h(2, k) ≤ 3 • 2 k-2 , for any k ∈ N -{0, 1}.

Lastly, we note that it is also possible to generate a covering cycle that does not have coincident Steiner points, except for the first and the last one, following a variation of the previous algorithm, as shown by Corollary 2.4 (see also Figures 4&5).

We can choose x k (S 3l+1 ) := l+2 2 to obtain

S 3l+1 ≡ C + l + 2 2 • #» e 2 ;
(20)

S 3l+2 ≡ C + 1 + 1 l • #» s l S 3l+1 + l + 1 2l + 3 • # » S 3l+1 S 3l+2 ≡ V 4l+1 ; S 3l+3 ≡ C -1 + 1 l • #» s l S 3l+2 + l + 1 2l + 3 • # » S 3l+2 S 3l+3 ≡ V 4l+2 ; S 3l+2 + l + 2 2l + 3 • # » S 3l+2 S 3l+3 ≡ V 4l+3 ; S 3l+4 ≡ C + l + 3 2 • #» e 2 S 3l+3 + l + 2 2l + 3 • # » S 3l+3 S 3l+4 ≡ V 4l+4 .
Therefore, for any given k ∈ N -{0, 1} we have provided a perfect covering cycle, for H(2, k), which is characterized by a link-length of 3 • 2 k-2 and such that no Steiner point is visited more than once, with the only exception of the starting/ending point,

S 1 ≡ S 1+3•2 k-2 .
This concludes the proof of Corollary 2.4.

Conclusion

Although for any H(2, k) we have constructively shown the existence of perfect covering cycles whose link-length is equal to 3 • 2 k-2 , the problem of finding an analogous formula concerning optimal covering paths for any given set H(n, k) := {1, 2, . . . , n} × {1, 2, . . . , n} × • • • × {1, 2, . . . , n} ⊂ R k , such that n ≥ 4 ∧ k ≥ 3, remains completely open [START_REF] Ripà | Minimum-link covering trails for any hypercubic lattice[END_REF] (e.g., we can only say that h(4, 3) ∈ {21, 22, 23} [START_REF]Sequence A318165 in the On-line Encyclopedia of Integer Sequences[END_REF][START_REF] Ripà | General conjecture on the optimal covering trails for any k-dimensional cubic lattice N k[END_REF]).

Figure 1 : 3 2

 13 Figure 1: The minimum-link perfect covering cycle C(h(2, 2)) := 1 2 , 3 2 -(-1, 0)-(2, 0)-1 2 , 3 2

Figure 2 :

 2 Figure 2: The minimum-link closed polygonal chain P(6) := 1 2 , 1 2 , 2 --1 2 , -1 2 , 0 -( 3 2 , 3 2 , 0)-1 2 , 1 2 , 2 -(-1 2 , 3 2 , 0)-( 3 2 , -1 2 , 0)-1 2 , 1 2, 2 visits all the nodes of H(2, 3) once and only once (picture realized with GeoGebra[START_REF] Hohenwarter | [END_REF]).

Figure 3 :

 3 Figure 3: The minimum-link closed polygonal chain P(12) := 1 2 , 1 2 , 1 2 , 3 2 -(-1, -1, -1, 0)-(2, 2, 2, 0)-1 2 , 1 2 , 1 2 , 3 2 -(-1, -1, 2, 0)-(2, 2, -1, 0)-

Acknowledgments

We sincerely thank Luca Onnis for his kind assistance on the initial phase of the present preprint.

, 2 joins all the nodes of H(2, 3) (picture realized with GeoGebra [START_REF] Hohenwarter | [END_REF]).

visits all the nodes of H(2, 4) once and only once (picture realized with GeoGebra [START_REF] Hohenwarter | [END_REF]).

Proof. We constructively prove the corollary by following the same approach that has been introduced in the proof of Theorem 2.3, taking also into account that the Steiner points of the type S 3l+1 , that lie on the straight line r, must have the coordinates x k (S 3l+1 ) different from each other.