Mauro Franceschelli
email: mauro.franceschelli@diee.unica.it.

‡ M Franceschelli

† P Frasca

Proportional Dynamic Consensus in Open Multi-Agent Systems

come

Proportional Dynamic Consensus in Open Multi-Agent Systems

Mauro Franceschelli ‡ and Paolo Frasca † Abstract-In this paper we investigate a dynamic consensus problem for an open multi-agent system. Open multi-agent systems are characterized by a time-varying set of agents connected by a network: agents may leave and new agents may join the network at any time, thus the term "open". The dynamic consensus problem consists in achieving agreement about the time-varying average of a set of reference signals that are assumed to be the agents' inputs. Dynamic consensus has recently found application in the context of distributed estimation for electric demand-side management, where a large population of connected domestic appliances needs to estimate its future average power consumption. Since the considered network of devices changes as new appliances log in and out, there is a need to develop and characterize dynamic consensus algorithms for these open scenarios. In this paper we give several initial contributions both to a general theory of open multi-agent systems and to the specific problem of dynamic consensus within this context. On the theoretical side, we propose a formal definition of open multi-agent system, a suitable notion of stability, and some sufficient conditions to establish it. On the applied side, we design a novel dynamic consensus algorithm, the Open Proportional Dynamic Consensus algorithm. We characterize some of its convergence properties in the proposed open-multi-agent systems framework and we illustrate its evolution by numerical simulations.

I. INTRODUCTION

A multi-agent system is a dynamical model for the behavior of a possibly large group of agents, e.g., robots, devices, sensors, oscillators etc., whose pattern of interactions due to sensing, communication or physical coupling is modeled by a graph that represents the network structure of the system. Most literature on multi-agent systems considers networks of fixed size, i.e., number of agents, and then considers several kinds of scenarios such as time-varying network topologies. In this paper we explicitly consider a more radical scenario of open multi-agent system where the set of agents is timevarying, i.e., agents may join and leave the network at any time.

Similar notions of "open" systems can be found in the computer science literature [START_REF] Huynh | An integrated trust and reputation model for open multi-agent systems[END_REF], [START_REF] Pinyol | Computational trust and reputation models for open multi-agent systems: a review[END_REF] when referring to software agents and the problem of evaluating reputation in open environments. Similar open systems have also been considered in game theory [START_REF] Shah | Dynamics in congestion games[END_REF]. Despite the abundance of works in multi-agent systems from the systems and control community, openness is rarely explicitly included in a rigorous analysis, but rather explored by simulations as in [START_REF] Zhu | Discrete-time dynamic average consensus[END_REF]. Indeed, openness implies some conceptual difficulties in adapting control-theoretic notions such as state or stability. For this reason, some authors have recently proposed to circumvent the mathematical hurdles by embedding the timevarying agent set in a time-invariant finite superset [START_REF] Varma | Open multi-agent systems with discrete states and stochastic interactions[END_REF]. In a different perspective, others have aimed to describe the open system through significant statistical properties: encouraging results have been presented in [START_REF] Abdelrahim | Maxconsensus in open multi-agent systems with gossip interactions[END_REF], where the authors study a max-consensus problem in open multi-agent systems, and in [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF], [START_REF]Open multi-agent systems: Gossiping with random arrivals and departures[END_REF], where the authors study the problem of averageconsensus by gossiping.

Our attempt to formalize open dynamical systems is motivated by modern power networks, where smart appliances that are endowed with communication and computation capabilities can arrive or depart at any time. In these applications, the appliances aim to estimate in a distributed fashion the current and future power consumption of the network they belong to: this information is then used locally to optimize the power consumption of the appliances [START_REF] Franceschelli | A heuristic approach for online distributed optimization of multi-agent networks of smart sockets and thermostatically controlled loads based on dynamic average consensus[END_REF], [START_REF] Franceschelli | Coordination of electric thermal systems for distributed demand-side management: A gossipbased cooperative approach[END_REF]. The use of dynamic consensus algorithms for distributed estimation is thus a key enabler for the development of distributed control protocols among appliances that do not require centralized supervision nor centralized data aggregation that could violate the privacy of the users. Lead by this application, we concentrate on one specific kind of open multi-agent dynamics, which generalizes to the open context the wellknow Proportional Dynamic Consensus dynamics. In this dynamics, each of the nodes receives an input signal and is tasked to track the average of all inputs over the network. In the smart grid scenario, the object of the distributed estimation is the time-varying average power consumption by the network. Thus, by considering the planned power consumption of each device as an external reference signal for each agent, a dynamic consensus algorithm can be used to estimate the time-varying average value of this potentially large set of reference signals. Since devices login and logout from the network without notice, the set of reference signals is, in general, time-varying.

The dynamic consensus problem has received significant attention, as demonstrated by the forthcoming tutorial [START_REF] Kia | Tutorial on dynamic average consensus: the problem, its applications, and the algorithms[END_REF]. Since the early work in [START_REF] Spanos | Dynamic consensus on mobile networks[END_REF], a fundamental idea to render consensus protocols "dynamic" has been adding the derivative of each agents' own reference signal to a consensus filter that would thus track the time-varying average of the references. Several algorithms that exploit this mechanism have been proposed [START_REF] Zhu | Discrete-time dynamic average consensus[END_REF], [START_REF] Nosrati | Dynamic average consensus via nonlinear protocols[END_REF], [START_REF] Kia | Distributed event-triggered communication for dynamic average consensus in networked systems[END_REF]: their main advantages are convergence speed and accuracy (which can be perfect for constant reference signals), while their main drawback is their lack of robustness with respect to errors in their initialization and, consequently, with respect to changes in the network composition. If the number of agents changes, these algorithms accumulate estimation errors that can severely deteriorate the estimation performance. Some algorithms, for instance those in [START_REF] Montijano | Robust discrete time dynamic average consensus[END_REF] and [START_REF] Freeman | Stability and convergence properties of dynamic average consensus estimators[END_REF], [START_REF] Bai | Robust dynamic average consensus of time-varying inputs[END_REF], [START_REF] Scoy | A fast robust nonlinear dynamic average consensus estimator in discrete time[END_REF], have instead shown superior robustness properties that can be useful to allow for the addition or removal of agents, even though their analysis has been so far limited to networks of fixed size. We take note that the main drawback of proportional dynamic consensus, as illustrated in [START_REF] Freeman | Stability and convergence properties of dynamic average consensus estimators[END_REF], consists in a tradeoff between steady-state error for constant reference signals, tracking error and convergence rate. In the recent work [START_REF] Franceschelli | Multi-stage discrete time dynamic average consensus[END_REF], the authors propose and characterize the so-called multi-stage dynamic consensus algorithm, which consists in a cascade of proportional consensus filters and has been proved to guarantee small steady-state and tracking error for a given convergence rate, thanks to exchanging a larger quantity of local information between the agents. Furthermore, the strategy proposed in [START_REF] Franceschelli | Multi-stage discrete time dynamic average consensus[END_REF] has been shown to be implementable with asynchronous and randomized (gossip-based) local interactions.

a) Contribution: The main contributions of this paper are: i) new definitions of equilibrium and stability, suitable for open-multi-agent systems, together with a sufficient condition to assess such stability; ii) a modified proportional dynamic consensus algorithm, presented in discrete time, which has been augmented with a mechanism to deal with agents joining and leaving the network: we thus refer to it as the Open Proportional Dynamic Consensus algorithm. The OPDC is the first dynamic consensus algorithms that is designed for open multi-agent systems.

b) Paper structure: In Section II we introduce both our general model of open multi-agent system and the specific Open Proportional Dynamic Consensus algorithm: the latter is a modified version of the proportional dynamic consensus algorithm, suitable to be executed by an open multi-agent system. Section III proposes the new definition of stability for the considered open multi-agent system, together with some sufficient conditions. These tools are then applied to the Open Proportional Dynamic Consensus algorithm. In Section IV the case study of distributed estimation in electric demand side management is presented to show a real case of application and numerical simulations are shown to corroborate the theoretical analysis. Finally, in Section V concluding remarks are given.

II. OPEN DYNAMICAL SYSTEMS

For all t ∈ Z ≥0 , let G(t) = (V (t), E(t)) be a timevarying directed graph with time-varying set of agents V (t) and time-varying set of edges E(t) ⊆ (V (t) × V (t)). Set V (t) contains the labels corresponding to the agents that are active at time t: we may think of V (t) as a subset of U, an infinite universe set of possible labels. The cardinality of set V (t), that is, the number of agents that belong to the network at time t, is denoted as n(t) = |V (t)|. Two agents i and j are said to be neighbors at time t if they share an edge at time t, i.e., (i, j) ∈ E(t). Let N i (t) be the set of neighbors of node i at time t, i.e., N i (t) = {j ∈ V (t) : (i, j) ∈ E(t)}.

For each time t and each agent v ∈ V (t), we associate a "state" variable x v (t) ∈ R and an "input" variable u v (t) ∈ R. Note that these variables are defined only at time instants such that v ∈ V (t). With these ingredients we can define laws that describe how the sequence [x(t)] t∈Z ≥0 evolves. Note that we will not be able to write x(t + 1) as a function of x(t), so that the evolution of x(t) does not constitute a "closed" dynamical system. Instead, we shall take as given the sequences [V (t)] t∈Z ≥0 , [E(t)] t∈Z ≥0 , and [u(t)] t∈Z ≥0 (with the consistency conditions that E(t) ⊆ (V (t) × V (t)) and u(t) ∈ R V (t) for all t) and define the evolution law of x by

x(t+1) = f (x(t), V (t), V (t+1), E(t), E(t+1), u(t)). (1)
Such update rule f should distinguish three kinds of nodes v, respectively belonging to the sets:

• R(t) = V (t) ∩ V (t + 1), i.e., remaining nodes that that belong both to V (t) and to V (t + 1); • D(t) = V (t)\V (t+1), i.e., departing nodes that belong to V (t) but not to V (t + 1); • J(t) = V (t + 1) \ V (t), i.e., joining nodes that belong to V (t + 1) but not to V (t).

Since all x(s) takes values in R V (s) , the components corresponding to D(t) are simply left out from x(t + 1). Instead, components in J(t) need to be "initialized" according to some rule. Finally, for all v ∈ R(t) there shall be a causal evolution law in the form

x v (t + 1) = f v (x(t), V (t), E(t), u(t)).
For concreteness, we now describe one such map f which we call the Open Proportional Dynamic Consensus (OPDC) algorithm.

Definition 2.1 (Open Proportional Dynamic Consensus)

Let ε > 0 and α ∈ (0, 1). At each time t ∈ Z ≥0 , the generic agent v ∈ V (t) updates its state according to the rules:

If v ∈ R(t), then x v (t + 1) =x v (t) -α(x v (t) -u v (t)) -ε w∈Nv(t) (x v (t) -x w (t)). ( 2a
)
If v ∈ J(t), then

x v (t + 1) = 1 |N v (t + 1)| w∈Nv(t+1)
x w (t).

(2b)

Under the assumption that V (t+1) = V (t), i.e., the set of agents does not change, the OPDC reduces to what is called Proportional Dynamic Consensus. Namely, it reduces to the update (2a), which can be written in matrix form as

x(t + 1) = x(t) -α(x(t) -u(t)) -L(t)x(t) = (1 -α)I -εL(t) x + αu(t) = P (t)x(t) + αu(t) (3)
where matrix

P (t) = (1 -α)I -εL(t) in R V (t)×V (t)
is a non-negative substochastic matrix with row sum and column sum equal to (1 -α) and spectral radius ρ = max {|1 -α -ελ| : λ is eigenvalue of L(t)}. In what follows, we shall assume that for all t the directed graph G(t) is strongly connected and that the parameter ε is small enough to ensure that the spectral radius of P (t) is 1 -α.

III. STABILITY ANALYSIS

We now define the concept of trajectory of points of interest which will be useful in the considered scenario of open multi-agent system. Definition 3.1 (Trajectory of points of interest) Consider an open system (1). Assume that for every t, the equation

y = f (y, V (t), V (t), E(t), E(t), u(t))
has a unique solution and denote that solution as x e (t). Then, the sequence [x e (t)] t∈Z ≥0 is called trajectory of points of interest of the open multi-agent system.

Observe that x e (t 0 ) ∈ R V (t0) represents the hypothetical equilibrium of the dynamics followed by x(t) if the three given sequences V (t), E(t) and u(t) would be kept constant for all t ≥ t 0 . Consequently, x e (t 0 ) is determined only by information at time t 0 : the time-variance of V (t), E(t) and u(t) does not imply any ambiguity in the definition of the sequence x e (t).

As an example, consider system (2). Under our standing assumptions on P (t), the solution x e (t) is unique for every t and can be computed as

x e (t) =(I -P (t)) -1 αu(t) = I + ε α L(t) -1 u(t).
Next, we define a so-called "open" distance function which is used to evaluate the distance between two points with labeled elements that belong to Euclidean spaces of different dimensions. In the particular case in which the two points have elements with the same labels, i.e., the same agents, the open distance function reduces to the Euclidean distance.

Definition 3.2 (Open distance function) Let V 1 and V 2 be two finite sets of node indices. Let d : R V1 × R V2 → R ≥0 be defined as d(x, y) = v∈V1∩V2 (x v -y v ) 2 + v∈V1\V2 x 2 v + v∈V2\V1 y 2 v (4) for any x ∈ R V1 and y ∈ R V2 .
The open distance (4) satisfies several natural properties, which we summarize in the next statement. The next definition introduces a notion akin to a weak form of Lyapunov stability for open multi-agent systems. Definition 3.5 (Stability of a trajectory of points of interest) Let x(t) be the evolution of an open system. A trajectory of points of interest x e (t) with bounded variation is said to be stable if there exists R ≥ 0 such that for every ε > R, there exists δ > 0 such that if d(x(0), x e (0)) < δ then d(x(t), x e (t)) < ε for every t ≥ 0.

The next result is a very simple sufficient condition to ensure stability in the above sense.

Theorem 3.6 (Stability of Open Multi-agent Systems)

Consider an open multi-agent system with state trajectory x(t) and having x e (t) as trajectory of points of interest of bounded variation B. If there exist γ ∈ (0, 1) and Γ ≥ 0 such that d(x(t + 1), x e (t + 1)) ≤ γd(x(t), x e (t)) + Γ, then the open multi-agent system is stable according to Definition 3.5 with R ≥ B+Γ 1-γ . Proof: By hypothesis, the trajectory of points of interest has bounded variations as per Definition 3.4, thus it holds d(x e (t + 1), x e (t)) < B for any t ≥ 0. By the assumptions of the theorem, it holds

d(x(t + 1), x e (t + 1)) ≤ d(x(t + 1), x e (t)) + d(x e (t + 1), x e (t)) ≤ γd(x(t), x e (t)) + Γ + B. This last quantity is smaller than d(x(t), x e (t)) if d(x(t), x e (t)) > R ≥ B + Γ 1 -γ . (5)
We then conclude that, as long as d(x(t), x e (t)) > R, d(x(t + 1), x e (t + 1)) < d(x(t), x e (t)).

Now, for any ε > R, we can simply choose δ < , so that by (5) it holds d(x(t), x e (t)) < ε for all t ≥ t 0 . This proves the stability of the considered open system according to Definition 3.5.

In order to apply this general result to OPDC algorithm, we need to make some assumptions about the impact of leaving and joining nodes. A simple such assumption would be the existence of Γ > 0 such that for all t

Γ 2 ≥ v∈D(t)
x 2 e,v (t) + v∈J(t)

x 2 v (t + 1).

We believe that this assumption is not overly restrictive. For the OPDC algorithm, if both the input signals and the number of agents that join and leave at each instant of time are bounded, then it is easy to find such a bound Γ. For instance, consider the case in which initial conditions and input signals are bounded in the interval [0, 1] and the number of agents that may join or leave is limited to one at each instant of time. Then, one can choose Γ = √ 2.

Proposition 3.7 (Stability of Open Proportional Dynamic

Consensus) The Open Proportional Dynamic Consensus is stable with R = B+Γ α if its trajectory of points of interest has bounded variation B and there exists Γ > 0 that satisfies [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF].

Proof: Let us define for all t the vector y(t + 1) = P (t)x(t): note that y(t + 1) ∈ R V (t) and that the restrictions to R(t) of x(t + 1) and y(t + 1) coincide. We then observe that d x(t + 1),x e (t + 1) where we used [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF] and γ = 1 -α. By assumption, we have d(x e (t + 1), x e (t)) ≤ B for all t ≥ 0 and the conditions in Theorem 3.6 are satisfied.

≤ d x(t + 1), x e (t) + d x e (t + 1), x e (t) ≤ v∈R(t) x v (t + 1) -x e,v (t) 2 + v∈D(t) x 2 e,v (t)

IV. NUMERICAL SIMULATIONS

In this section we simulate numerically the open proportional dynamic consensus protocol with tuning parameter α = 0.01. We consider the problem of estimating the timevarying average of the power consumption by a network of smart appliances at a given time of the day. We consider a network composed at time t = 0 by 200 agents, each with constant input u i (t) chosen uniformly at random in the interval [0, 1], representing the power consumption of device i at a given hour of the day. We numerically simulate two scenarios. In the first, called Scenario A, at each instant of time t one agent leaves the network with probability p l = 0.05 and one joins the network with probability p j = 0.1, thus resulting in a network which, in expectation, increases in size. In the second, called Scenario B, at each instant of time t one agent leaves the network with probability p l = 0.1 and one joins the network with probability p j = 0.1, thus resulting in a network which, in expectation, does not change size. In both scenarios at each instant of time the network topology is that of a random graph with probability of edge existence between any two nodes equal to p = 0.1.

We now discuss the numerical simulation of Scenario A where an agent is more likely to join than to leave. In Figure 1 it is shown the time-varying number of agents. In Figure 2 it is shown the average of the inputs ū(t) = 1 n(t)

v∈V (t)
u v (t), which is time-varying due to joining and leaving agents. In Figure 3 it is shown the evolution of the open distance function d(x(t), x e (t)), i.e., the distance between the state of the network at time t and the point of interest x e at time t. Finally, in Figure 4 it is shown the evolution of d(x(t), ū(t)), which is a natural measure of the estimation error in dynamic consensus algorithms. From the simulations it can be seen that the distance functions are monotonically decreasing despite the time-varying number of agents up to threshold value which is determined by the join/leave dynamics of the agents. We now discuss the numerical simulation of Scenario B where agents are equally likely to join or leave. In Figure 5 it is shown the time-varying number of agents. In Figure 6 it is shown the evolution of the open distance function d(x(t), x e (t)), i.e., the distance between the state of the network at time t and the point of interest x e at time t. Finally, in Figure 7 it is shown the evolution of d(x(t), ū(t)). It can be seen that also in this scenario the distance functions are monotonically decreasing up to threshold value which is determined by join/leave dynamics of the agents. Open distance function value Open distance function value Open distance function value a distributed estimation problem in the context of electric demand-side management, where a large population of smart appliances interconnected via a peer-to-peer network over the internet needs to estimate the average of the future planned power consumption in an environment where devices may log-in and log-out without notice. Thus motivated, we proposed a novel definition of stability for open multi-agent systems and proposed a dynamic consensus algorithm, the Open Proportional Dynamic Consensus algorithm, which enables the agents to estimate in distributed fashion the timevarying average of a time-varying set of reference signals.

Preliminary results on the characterization of the stability of the proposed algorithm have been provided. Numerical simulations show that the proposed algorithm is robust to time-varying number of agents. This paper opens several research avenues. On the theoretical side, we have only scratched the surface of a possible theory of open dynamical systems. On the practical side, other dynamic consensus protocols can be envisioned and tested in concrete estimation problems.

Proposition 3 . 3 ( 3 )

 333 Open distance function) Function d(x, y) in (4) is such that for any vectors x, y, and z:1) d(x, y) ≥ 0; 2) d(x, y) = d(y, x); If x = y, then d(x, y) = 0; 4) d(x, z) ≤ d(x, y) + d(y, z)Note that the converse of the third implication does not hold. Indeed, consider x ∈ R {1,2} to be x = [1, 0] and y ∈ R {1} to be[START_REF] Huynh | An integrated trust and reputation model for open multi-agent systems[END_REF]. Then, d(x, y) = 0 despite the two vectors being different.Having this open distance available, we can naturally use it on the trajectories of points of interest. Definition 3.4 (Trajectory of points of interest of bounded variation) A trajectory of points of interest x e (t) is said to have bounded variation if there exists a constant B ≥ 0 such that d(x e (t + 1), x e (t)) ≤ B for all t ∈ Z ≥0 where d(•, •) is the open distance function in (4).

Fig. 1 .Fig. 2 .

 12 Fig. 1. Scenario A: Time-varying number of agents |V (t)|

Fig. 3 .

 3 Fig. 3. Scenario A: Evolution of open distance function d(x(t), xe(t)).

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Scenario A: Evolution of open distance function d(x(t), ū(t)) (estimation error).

Fig. 7 .

 7 Fig. 7. Scenario B: Evolution of open distance function d(x(t), ū(t)) (estimation error).

  +

		x 2 v (t + 1) + B
		v∈J(t)
	≤	y v (t + 1) -x e,v (t)	2
	v∈R(t) D(t)
	+	x 2 e,v (t) +	x 2 v (t + 1) + B
	v∈D(t)	v∈J(t)
	≤ P (t)x(t) + αu(t) -x e (t) 2 + Γ + B
	= P (t)x(t) -P (t)x e (t) 2 + Γ + B
	≤ P (t) 2 x(t) -x e (t) 2 + Γ + B
	≤ (1 -α) x(t) -x e (t) 2 + Γ + B
	= γd(x(t), x e (t)) + Γ + B,

‡ This work was supported by the Italian Ministry of Research and Education (MIUR) under call SIR 2014 with the grant "CoNetDomeSys", code RBSI14OF6H, and by Regione Autonoma della Sardegna with the 2016/2017 "Visiting Professor Program".


Laurens R Krol
email: lrkrol@gmail.com

Juliane Pawlitzki

Fabien Lotte

Klaus Gramann

Thorsten O Zander

SEREEGA: Simulating event-related EEG activity

Keywords: Electroencephalography, Simulation, Evaluation, Event-Related, Classification, Brain-Computer Interface, Ground Truth Highlights:

published or not. The documents may come

Introduction

Having seen almost a century of continuous research and development since its first application on humans in the 1920s [START_REF] Berger | Über das Elektrenkephalogramm des Menschen[END_REF], electroencephalography (EEG) is now widely used in, among others, clinical settings, neuroscience, cognitive science, psychophysiology, and brain-computer interfacing, while its use continues to expand in fields such as neuroergonomics [START_REF] Parasuraman | Neuroergonomics: the brain at work[END_REF][START_REF] Frey | Framework for electroencephalography-based evaluation of user experience[END_REF], neurogaming [START_REF] Krol | Meyendtris: A hands-free, multimodal Tetris clone using eye tracking and passive BCI for intuitive neuroadaptive gaming[END_REF], neuromarketing [START_REF] Vecchiato | On the use of EEG or MEG brain imaging tools in neuromarketing research[END_REF], neuroadaptive technology [START_REF] Zander | Neuroadaptive technology enables implicit cursor control based on medial prefrontal cortex activity[END_REF] and mobile brain/body imaging [START_REF] Gramann | Cognition in action: imaging brain/body dynamics in mobile humans[END_REF]. As of December 2017, PubMed reported over 140 000 publications related to EEG, with over 4 000 published in each of the past five years.

EEG reflects the electric fields that arise primarily due to the synchronous activity of postsynaptic potentials at apical dendrites in the cortical surface of the brain, as recorded by electrodes placed on the scalp. As such, it measures a specific subset of brain activity. This enables cognitive and affective correlates to be found in the EEG, allowing post-recording or even real-time evaluation of certain mental states exhibited by the recorded person, such as surprise [START_REF] Donchin | Surprise![END_REF], error perception [START_REF] Falkenstein | Effects of errors in choice reaction tasks on the ERP under focused and divided attention[END_REF][START_REF] Blankertz | Single Trial Detection of EEG Error Potentials: A Tool for Increasing BCI Transmission Rates[END_REF], task load [START_REF] Klimesch | EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis[END_REF][START_REF] Mühl | EEG-based workload estimation across affective contexts[END_REF][START_REF] Zander | Automated task load detection with electroencephalography: Towards passive brain-computer interfacing in robotic surgery[END_REF], or imagined movement [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF][START_REF] Blankertz | The non-invasive Berlin Brain-Computer Interface: Fast acquisition of effective performance in untrained subjects[END_REF]. Compared to other brain monitoring and imaging methods, EEG is relatively inexpensive and provides a high temporal resolution. It is also becoming increasingly portable and ready for personal use [START_REF] Zander | Evaluation of a dry EEG system for application of passive braincomputer interfaces in autonomous driving[END_REF] in various forms of neuroadaptive technology [START_REF] Krol | Passive Brain-Computer Interfaces: A Perspective on Increased Interactivity[END_REF]. This explains EEG's relative popularity.

However, these advantages present a trade-off, with costs incurred primarily in spatial resolution. A single electrode records the average activity of up to a billion neurons, and probably never less than 10 million neurons [START_REF] Nunez | Electric fields of the brain: the neurophysics of EEG[END_REF]. This and other issues including volume conduction, the placement and distance of the electrodes relative to the cortical generators of the activity they measure, and the complex relation between cortical functions and features of scalp potentials, require that great care is taken when analysing and interpreting raw EEG recordings.

A host of methods have been developed over the last decades to extract robust features from the recorded EEG that correlate to cortical functions and can be understood in a neurophysiological or statistical sense. Examples are the event-related potential technique [START_REF] Luck | An Introduction to the Event-Related Potential Technique[END_REF], independent component analysis [START_REF] Makeig | Independent component analysis of electroencephalographic data[END_REF], common spatial patterns [START_REF] Guger | Real-time EEG analysis with subject-specific spatial patterns for a brain-computer interface (BCI)[END_REF], hierarchical linear modelling [START_REF] Pernet | LIMO EEG: A toolbox for hierarchical linear modeling of electroencephalographic data[END_REF], and a many signal processing and machine learning algorithms [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF].

One major difficulty in developing techniques for EEG analyses is that there is no ground truth that describes the exact brain activity: The recorded EEG data cannot be compared to the actual neuro-electric activity of the brain because no technique exists to provide such a reference measurement. Thus, developers must use other ways to examine the validity of their EEG analysis approaches. The exact nature and detail of the required ground truth of course differs depending on the analysis.

To that end, simulated EEG data (or "toy data") has often been used to test and validate methods, as for example with blind source separation [START_REF] Makeig | Independent component analysis of simulated ERP data[END_REF][START_REF] Potter | Separation performance of ICA on simulated EEG and ECG signals contaminated by noise[END_REF], connectivity measures [START_REF] Silfverhuth | Experimental comparison of connectivity measures with simulated eeg signals[END_REF][START_REF] Stam | Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF], artefact removal [START_REF] Romo-Vazquez | EEG ocular artefacts and noise removal[END_REF][START_REF] He | Removal of ocular artifacts from the EEG: a comparison between time-domain regression method and adaptive filtering method using simulated data[END_REF], functional brain imaging [START_REF] Gramfort | Functional brain imaging with m/eeg using structured sparsity in time-frequency dictionaries[END_REF] and neurophysiological weight vector interpretation [START_REF] Haufe | On the interpretation of weight vectors of linear models in multivariate neuroimaging[END_REF]. The authors of these examples all implemented simulation approaches from scratch, usually by linearly mixing a number of independent signals. This linear mixing was done using random weights, i.e., no realistic spatial information was taken into account. Other approaches use head models in order to provide more realistic linear mixing and add spatial dependencies to the simulation (e.g., [START_REF] Giraldo | Estimation of dynamic neural activity using a Kalman filter approach based on physiological models[END_REF][START_REF] Haufe | Modeling sparse connectivity between underlying brain sources for EEG/MEG[END_REF][START_REF] Haufe | A critical assessment of connectivity measures for EEG data: A simulation study[END_REF][START_REF] Huiskamp | Interindividual variability of skull conductivity: an EEG-MEG analysis[END_REF].

Such custom-made approaches are often difficult to reproduce, as they have been implemented using different software packages, are reported at different levels of abstraction, and/or may be using head models that are not publicly available. Some authors have not implemented their own methods, but have instead relied on commercially available packages (e.g., [START_REF] Yao | Evaluation of different cortical source localization methods using simulated and experimental EEG data[END_REF][START_REF] Lansbergen | ADHD and EEG-neurofeedback: a double-blind randomized placebo-controlled feasibility study[END_REF]. These packages, however, are not specifically designed for the purpose of simulation, nor are they freely available to the general public or open to scrutiny.

One open-source package known to the authors that provides simulation functionality is the Source Information Flow Toolbox (SIFT; [START_REF] Delorme | ERICA: New Tools for Advanced EEG Processing[END_REF]. Its main purpose is to investigate structural, functional, and effective connectivity between brain regions and networks, but the toolbox can also be used to simulate scalp EEG using a system of coupled oscillators. Such simulations can serve as effective ground truth for the same connectivity measures that the toolbox investigates, and have been used to evaluate blind source separation methods as well [START_REF] Hsu | Online recursive independent component analysis for real-time source separation of high-density EEG[END_REF]. The SIFT toolbox, however, is not intended for wider-purpose simulation and as such, the methods that can be meaningfully applied to the data it can generate are limited.

Another work allows a number of predefined responses (e.g. a P300, an eye blink, a frequency spike) to be configured along a given timeline of events [START_REF] Lindgren | SimBCI -tool to simulate EEG and BCI[END_REF]. It focuses on effects specific to brain-computer interface (BCI) applications and is aimed at the batch analysis of BCI classification algorithms.

In a more normative approach, [START_REF] Haufe | A simulation framework for benchmarking eeg-based brain connectivity estimation methodologies[END_REF] proposed a simulation and evaluation framework and made their code publicly available. The framework includes forward modelling using a realistic head model, both source and sensor noise with controlled signal-to-noise ratios, and signal generation based on autoregressive models. As with SIFT, this approach focuses solely on connectivity measures on continuous data, and thus provides no further signal generation options.

A final EEG simulation approach of note is the phantom head created by [START_REF] Oliveira | Induction and separation of motion artifacts in EEG data using a mobile phantom head device[END_REF]. They constructed a physical head model, filling a mannequin head with a conductive plaster and inserting antennae to simulate current dipoles. This allows unique effects to be investigated: For example, the authors investigated the effects of head motion by placing the model on a motion platform. Constructing such hardware, however, is not a viable approach for most questions where simulation can provide an answer. Importantly, software is easier to share, maintain, adapt, and extend.

Clearly, EEG data simulation is widely used as a tool to assess and validate the methods that are in use and in development. However, to the best of the authors' knowledge, there currently exists no software package whose sole or primary purpose it is to simulate different, configurable, and/or custom types of EEG data, i.e. a dedicated, general-purpose EEG data simulation tool. We therefore present SEREEGA, short for Simulating Event-Related EEG Activity, making EEG data simulation more accessible to researchers.

SEREEGA is a free and open-source MATLAB-based toolbox to generate simulated, eventrelated EEG data. Starting with a forward model obtained from a head model or pre-generated lead field, dipolar brain components can be defined. Each component has a specified position and orientation in the head model. Different activation patterns or signals can be assigned to these components. Scalp EEG data is simulated by projecting all signals from all components onto the scalp and summing these projections together.

SEREEGA is modular in that different head models and lead fields can be supported, as well as different activation signals. Five lead fields are currently supported, four pre-generated, the fifth customisable according to the user's needs from a standard head model. Five types of activation signals are provided, allowing the simulation of different types of systematic (eventrelated) activity in both the time and the frequency domain, as well as the inclusion of any already existing time series as an activation signal.

This toolbox is intended to be a tool to generate data with a known ground truth in order to evaluate neuroscientific and signal processing methods, such as blind source separation, source localisation, connectivity measures, brain-computer interface classifier accuracy, derivative EEG measures, et cetera.

In the following, we first introduce the architecture of the toolbox and provide a brief introduction to its basic functionality and workflow. We then provide an analysis using established neuroscientific and brain-computer interfacing methods of a sample data set created with the toolbox.

2 SEREEGA Architecture and Functionality

Principles of EEG Simulation

To simulate EEG data, the toolbox solves the forward problem of EEG, prescribing how activation signals from specific sources in the brain are projected onto an array of electrodes on the scalp. In matrix notation, this can be written as x = As + , with x denoting the vector of the recorded or simulated scalp signal, s the source activation signal, A the projection matrix used to project signals from the source to the scalp electrodes, and denoting a vector of noise.

In SEREEGA, the user defines s, A, and , allowing EEG data x to be simulated. Source activations are defined on the basis of so-called components: for each component, any number of different signal classes can be defined, which prescribe how corresponding activation signals are to be generated. For each signal class, a type of activation signal (for example, an event-related potential; see Section 2.3.3) is defined, as well as all corresponding parameters. Furthermore, for each component, a source from the lead field is specified, which is modelled as a dipole at a specific location in the brain. The component also contains this dipole's (i.e. this source's) orientation. A component thus prescribes what signal is to be simulated, and how it is to be projected onto the scalp.

Notation Description

 := A •,h,• ∈ R n chs ×3 Projection matrix of one source A ∈ R n chs ×nsrcs×3
SEREEGA simulates one segment of EEG data, or one epoch, at a time, and repeats this n eps times to obtain a larger data set. In the following, we first consider a single epoch. The activation signal of each class, g k c ∈ R ntps , is considered as a vector which consists of the activation signal (i.e. the amplitude time series) for all n tps time points in this epoch.

Thus, for each component k, its activation s k consists of the summed activations of all signal classes g k c assigned to that component, which are n k cls many:

s k = n k cls c=1 g k c .
For most signal classes, its exact activation signal g k c is determined procedurally at runtime based on the specified parameters.

We now describe how an activation signal is projected to the channels. For this, we consider the signal at a single time point t, by taking the t-th entry of s k , denoted by s k t . We project this signal from a source denoted by h.

The activation signal is projected onto the electrodes on the scalp through a projection matrix  which we obtain using lead field theory [START_REF] Ferree | Lead field theory and the spatial sensitivity of scalp EEG[END_REF]. A lead field contains projection parameters that indicate for each electrode and source how an activation at that source is scaled when it is recorded at that electrode on the scalp. To that end, the activation signal is split in three directions (represented by the three base vectors of Euclidean space), resulting in a third-order tensor consisting of n chs layers, n srcs rows and 3 columns:

A ∈ R n chs ×nsrcs×3 .
Each row describes the projection matrix for a single source. For a specific source h corresponding to a component k, Â := A •,h,• ∈ R n chs ×3 thus describes the projection matrix of that component. The orientation o h of source h can be expressed as a linear combination of the different base vectors contained in the three columns of Â:

(1, 0, 0) T the x-direction, pointing to the left ear, (0, 1, 0) T the y-direction, pointing to the nose, and (0, 0, 1) T the z-direction, pointing to the top of the head.

This can be expressed by

o h = o h 1   1 0 0   + o h 2   0 1 0   + o h 3   0 0 1   =   o h 1 o h 2 o h 3   , for o h 1 , o h 2 , o h 3 ∈ R.
To project a component's activation onto the scalp, this activation s k t is first oriented by scaling it in the corresponding directions, yielding the oriented activation vector ŝk t ∈ R 3 :

ŝk t := s k t • o h =   s k t • o h 1 s k t • o h 2 s k t • o h 3   .
Then, the oriented activation vector ŝk t is projected through the lead field which corresponds to the source of the activation, by multiplication with the projection matrix  = [a ij ] i=1,...,n chs j=1,...,3 :

xk t := Â • ŝk t =   3 j=1 a ij ŝk t   i=1,...,n chs .
This yields a vector xk t ∈ R n chs in which every element corresponds to the simulated signal amplitude at time point t, projected from source h to one electrode. To obtain the corresponding matrix for all time points for the component k, the vectors of all time points are concatenated:

Xk = xk 1 , xk 2 , xk 3 , • • • , xk ntps .
X e = ncmps k=1
Xk .

The parameters in the signal classes describe how each component's activation signals vary between epochs. The projection matrix  can also be made to vary from epoch to epoch, to simulate non-stationarities in the components' projections, for example due to shifts in electrode positions. Finally, all scalp signal matrices for all epochs are concatenated in the third dimension, yielding a third-order tensor in R n chs ×ntps×neps . It is at this point that sensor noise E is optionally added (source-level noise is defined as an element of g k c ). This yields the tensor of simulated EEG data:

X = X 1 , X 2 , X 3 , • • • , X neps + E.
This is the data format used by most software packages to represent epoched EEG data.

Platform and License

MATLAB R2014b or higher is recommended for SEREEGA. Some optional functions depend on the Digital Signal Processing (DSP) toolbox version 8.6 (R2014a) or higher. EEGLAB (Delorme & Makeig, 2004) is required as it is used for a number of functions. Lead field generation either requires additional head model files which can be downloaded from their respective websites, or the FieldTrip toolbox [START_REF] Oostenveld | FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data[END_REF]. Since SEREEGA is modular, future functions may have further dependencies.

SEREEGA is licensed under the GNU General Public License, version 3, and the code is publicly available on GitHub1 . This manuscript refers to version 1.0.5-beta.

Terminology and Workflow

SEREEGA is available as an EEGLAB plug-in including a graphical user interface (GUI) that allows the core steps of designing and running a simulation to be performed; see figure 1. For more advanced use, SEREEGA is based on written commands and assignments. A general configuration variable holds information about the number of epochs to simulate, their length, and their sampling rate. This configuration, as well as many other variables, are contained within MATLAB structure arrays.

SEREEGA's forward model is contained in a lead field structure array. The lead field structure contains all possible sources within a virtual brain, modelled as dipoles at specific locations along with their projection patterns. The word source thus corresponds to an index in this lead field. For each such index, the lead field contains a position (x, y, z), and projection patterns along three axes (x, y, z). As mentioned in Section 2.1, a linear combination of the three projection patterns can effectively be used to virtually define the source dipole's orientation in 3D space. A source together with its orientation can be saved as a component structure array, which additionally needs to be assigned a signal : the simulated neuro-electrical activity that will be projected from it, such as an event-related potential or an oscillation at a specific frequency. These are defined separately and added to the components.

The main workflow consists of defining any number of such components, each containing at least one source, orientation, and signal. The simulation of scalp data finally generates each component's signals and projects them onto the scalp.

The following sections go through the main steps in more detail. Note that these steps must not necessarily be followed in this order.

Lead Field Generation

The lead field determines both the maximum possible number of sources as well as the number of channels that will be simulated. Currently, SEREEGA supports two processes to obtain a lead field: it can be obtained from an existing, pre-generated lead field, or a custom lead field can be generated from a given head model. Currently, support for four pre-generated lead fields is included.

The New York Head (ICBM-NY) pre-generated lead field includes almost 75 000 source locations and their projections onto to 228 channels [START_REF] Huang | The New York Head-a precise standardized volume conductor model for EEG source localization and tES targeting[END_REF]. The electrode positions follow the international 10-05 system [START_REF] Oostenveld | The five percent electrode system for highresolution EEG and ERP measurements[END_REF] but also include two rows of channels below the ears. Each source comes with a default orientation that orients it perpendicular to the cortical surface.

The Pediatric Head Atlases comprise three different head models with pre-generated lead fields for up to three different electrode layouts each, ranging from 128 to more than 2500 electrodes [START_REF] Song | Anatomically accurate head models and their derivatives for dense array EEG source localization[END_REF]. The models cover heads from three pediatric age clusters, 0-2, 4-8, and 8-18 years old. Lead field sources are spaced in an approximately 1 × 1 × 1 mm grid and range from 3188 to 4837 in number. These lead fields do not contain default dipole orientations. For inclusion in SEREEGA, the given electrode and dipole coordinates are transformed upon initialisation to be centred around (0, 0, 0) and aligned to the axes.

FieldTrip [START_REF] Oostenveld | FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data[END_REF] can be used to generate a lead field as needed. Using a given head volume, it can generate any number of sources with a given resolution, projecting to any number of channels. By default, a standard head model and a 339-channel definition file following the international 10-05 system are included. FieldTrip does not provide default source orientations.

When obtaining a lead field, it is possible to only use a subset of the available channels by indicating their labels. Convenience functions are available to simulate standard EEG cap layouts with e.g. 64 channels. Figure 2 shows the available channel locations for three head models.

Source Selection and Orientation

The lead field structure contains all possible sources. From these, one or more sources can be selected using different methods to be included in the simulation. It is possible to simply obtain a random source, or to obtain multiple random sources that are at least a given distance apart from each other. It is also possible to select the source nearest to a specific location in the brain, or all sources within a certain radius from a specific location or another source. Sources are referenced using their index in the lead field. Figure 3, left, shows the location of one source using EEGLAB's standard head model and corresponding plotting function. (In this case, EEGLAB's standard head model fits the lead field's head model. SEREEGA's default plotting functions do not depend on this being the case.)

Sources represent dipoles at the indicated location. Having selected a source location, this dipole can be oriented to face different directions, which determines how it projects onto the scalp. The lead field makes this possible by containing the dipole's projection pattern onto all selected electrodes into three mutually perpendicular directions x, y, and z-the base vectors mentioned in Section 2.1. A linear combination of these three projection patterns can then be used to effectively orient the dipole in space. Thus, a dipole's orientation is indicated using three values, with (1, 0, 0) representing a perfect x-positive orientation, (0, 1, 0) y-positive, ( 1 √ 2 , 1 √ 2 , 0) a 45-degree angle with the same amplitude, et cetera. Figure 3, right, shows the three separate patterns as well as a combined projection pattern for one source.

Signal Definition

Having selected a source location and orientation, i.e. a projection pattern, one or more signals can be defined to be projected onto the scalp as if they originate from that location. These are defined using classes. A class contains a signal's type and any parameters corresponding to that type in a structure array. For each simulated epoch, a signal will be procedurally generated using the parameters of that class.

SEREEGA currently includes four types of signals: systematic deflections in the time domain (i.e., event-related potentials), systematic modulations of oscillatory activity (i.e., event-related spectral perturbations), stochastic processes in the time domain (i.e., noise), and autoregressive signals (i.e., spatio-temporal signals with systematic interactions between them). Figure 4 provides an overview of three of these and their main parameters. A final data type enables the inclusion of any already existing time series as an activation signal. The following sections mention the primary free parameters for each type.

Event-Related Potentials An event-related potential (ERP; [START_REF] Luck | An Introduction to the Event-Related Potential Technique[END_REF] class defines one or more positive and/or negative 'peaks' or deviations from a baseline in sequence. Each peak is determined by its latency in ms, its width in ms, and its amplitude in µV. For example, the ERP in figure 4 is a single positive peak at 500 ms, 200 ms wide, of 1 µV.

A peak is generated by centring a normal probability density function around the indicated latency with the given width covering 6 standard deviations. This is then scaled to the indicated amplitude. For multiple peaks, each peak is generated individually and then summed together. From top left to bottom right: projection in the x-positive, y-positive, and z-positive direction, and a linear combination of the three (arbitrary units). This effectively orients the dipole at the source's location in 3D space.

Event-Related Spectral Perturbations Where ERPs represent systematic activity in the time domain, event-related spectral perturbation (ERSP) classes are defined primarily using spectral features. The term ERSP refers to event-related spectral power changes [START_REF] Makeig | Mining eventrelated brain dynamics[END_REF]; in SEREEGA, it refers to the class of signals that generate specific spectral powers and/or changes therein.

First, a base frequency is defined, either as a pure sine wave of a single frequency, or as a frequency band. In this latter case, uniform white noise is band-pass filtered in the indicated band using a window-based finite impulse response filter with a Kaiser window and an automatically decided filter order [START_REF] Oppenheim | Discrete-time signal processing[END_REF]. The left middle panel of figure 4 shows an oscillation of 10 Hz, the panel to the right of that shows filtered noise with power focused between 5 and 10 Hz.

The minimal definition of an ERSP class consists of a frequency (or frequency band) and its amplitude. In case of a single frequency, a phase can optionally be indicated as well.

In other cases, these parameters merely define a base oscillation that is further modulated. The lower row of figure 4 presents two types of modulations that can be defined for ERSP classes: (inverse) burst modulation and amplitude modulation. The former attenuates or amplifies the signal to a given degree using a Tukey window of given width, latency, and tapering. A Tukey window can be adjusted ranging from a square to a cosine shape. This mimics event-related synchronisation and desynchronisation. The latter modulates the amplitude of the base signal to a given degree according to the phase of a sine wave, whose frequency and phase can be defined in the class. Additionally, a pre-stimulus period can be indicated such that the modulation only starts at a given latency.

Noise A class of a noise signal simply produces coloured noise, from either a normal Gaussian or a uniform distribution, of a given amplitude. A noise class is defined by at least a colour and an amplitude. The upper right panel of figure 4 shows white noise with an amplitude 1 µV.

The Gaussian function uses MATLAB's DSP toolbox to generate noise with a spectral characteristic of 1/f n , with n either -2, -1, 0, 1, or 2. The uniform function first samples the signal from a uniform distribution and then applies a Fourier transformation to achieve the same spectral characteristic.

Data For the inclusion of any given pre-existing time series, for example data generated elsewhere or taken from existing EEG recordings, a data class can be defined. It must be given a matrix of data, epoch-dependent indices, and an absolute or relative amplitude. It will project the indicated data scaled to the given amplitude for each simulated epoch.

Autoregressive model A final type of class generates signals based on a linear autoregressive model (ARM), following the approach by [START_REF] Haufe | A simulation framework for benchmarking eeg-based brain connectivity estimation methodologies[END_REF]. For a single source activation, this generates a time series where each sample depends linearly on its preceding samples plus a stochastic term. At least the order of the model, i.e. the number of non-zero weights on the preceding samples, and the signal amplitude must be configured. The exact weights of the model are determined automatically.

It is also possible to simulate connectivity between multiple sources. Uni-and bidirectional interactions can be indicated such that a time-delayed influence of one source on another is included in the signal. Either the number of interactions is configured and the toolbox will randomly select the interactions, or a manual configuration of the exact directionalities can be indicated. The model order is the same for all interactions.

Interacting signals are not simulated at runtime, but generated beforehand and included in the simulation as data classes.

Variability

It defeats the purpose of simulating multiple epochs if every epoch were to be exactly the same. To take an ERP as example, it may not be desirable to have a peak defined at 500 ms appear at exactly 500 ms in every epoch. In fact, perhaps the ERP should sometimes not appear at all. To that end, both classes and components contain parameters to add variability to the generated signal.

For all parameters displayed in figure 4 controlling the shape of the various signals, two more parameters exist that control those base parameters: each parameter has a corresponding deviation parameter, and a slope parameter.

The deviation indicates a range of allowed values. A base value of 2 with a deviation of 0.3, for example, will vary from epoch to epoch between 1.7 (2 -0.3) and 2.3 (2 + 0.3) following a Gaussian function. Specifically, the range of deviation covers 6 standard deviations. This can be used to model general Gaussian variability within the system. A deviation is calculated for each parameter value separately for each epoch. For some parameters, one additional shift parameter is available to make these values deviate with equal magnitude within one epoch. This can be useful to e.g. have the latency of multiple peaks in an ERP vary together, maintaining the overall shape of the ERP.

The slope indicates a systematic change over time. A base value of 2 with a slope of -0.3, and no deviation, will be 2 on the first epoch, and 1.7 (2 -0.3) on the final epoch. In between, it will scale linearly with the number of the current epoch. Slopes can be used to model e.g. effects of fatigue and habituation, or to cover a range of intended values.

Each class also has a probability parameter that determines, for each epoch, that signal's probability of appearance. A probability of 0.5 indicates a 50% chance of the defined signal being generated for any simulated epoch. In case it is not generated, a flatline is returned instead, meaning it does not interfere with any other activations at that time.

Deviation and slope parameters are also available for components, to control the dipole orientation. Components can have one additional source of variability, with respect to the source location itself. Instead of one source index, multiple source indices can be indicated as that component's source. In that case, one of them will be randomly picked for each epoch. As such, together with the orientation variability, spatial variability can be added as well.

Simulating Scalp EEG

A component is a structure array that contains, for each element, at least one source index, a corresponding orientation, and at least one signal. Thus, by combining the results of the previous two steps into components, it is specified which of the defined signals are projected from which source at which orientation. It is possible to assign more than one signal to a component. To add noise to an ERP, for example, both an ERP and a noise class can be added to the same component. All signals assigned to a component are summed together before being projected. When the same signal class is given to multiple components, a new instance of that signal will be generated for each component at runtime.

There is no limit to the number of components that can be defined, nor to the number of signal classes that can be added to each component.

A simulation function takes the defined components, the lead field, and the general configuration as input, and outputs the simulated scalp data in a channels × samples × epochs matrix, as well as components × samples × epochs source data. For each epoch, this function follows the forward model described in Section 2.1 to produce the scalp data. That is, it generates and sums each component's signals, projects them onto the scalp using the given lead field, source location, and orientation, and then sums all the projected activations. Each signal is generated independently of the others. For interacting signals, as with ARM classes, the signals are generated beforehand and included as data classes.

At this stage, a final layer of noise can optionally be added to the simulated scalp data to simulate sensor noise. This is uniform, temporally and spatially uncorrelated white noise.

Manual and Random Definition of Signal Classes

As a final point it is worth mentioning that not all of the available parameters must necessarily be set by the user. Only a small number are required, and a validation function automatically sets all others to their default values. All deviation-related values default to zero. Non-zero default values are chosen to be most useful; the phase of an oscillation, for example, defaults to a new random value for each epoch. A convenience function exists to add deviations and/or slopes of a certain percentage of each base value to signal classes. When the lead field contains default source orientations, these will be used automatically when no other orientation has been explicitly indicated.

When multiple components are to be simulated, it is also not necessary to manually define all signals. For one, the same signal can be added to multiple components. Also, convenience functions exist to generate any number of 'random' ERP or ERSP classes, based on a range of allowed base values.

Sample Code

The following lines of code reflect the workflow explained above, generating 100, 1-second epochs of brown noise from 64 sources spaced at least 25 mm apart, projected onto 64 channels.

epochs

= struct('n', 100, 'length', 1000, 'srate', 1000); leadfield = lf_generate_fromnyhead('montage', 'S64'); sourcelocs = lf_get_source_spaced(leadfield, 64, 25); signal = struct('type', 'noise', 'color', 'brown', ... 'amplitude', 1, 'amplitudeDv', .5); components = utl_create_component(sourcelocs, signal, leadfield); data = generate_scalpdata(components, leadfield, epochs); EEG = utl_create_eeglabdataset(data, epochs, leadfield); EEG = utl_add_icaweights_toeeglabdataset(EEG, components, leadfield);

The epochs structure array contains three fields configuring the number of epochs, their length in ms, and the sampling rate in Hz. On the next line, a lead field is obtained from the New York Head using a preconfigured montage of 64 channels. Next, 64 source locations are randomly selected from the lead field, such that each source is at least 25 mm apart from all the others. As a signal, brown noise is used with an amplitude varying between 0.5 and 1.5 for each epoch. 64 components are then created by combining that same signal with the 64 previouslydefined source locations. Because the New York Head is used, no orientations must be indicated explicitly: the toolbox will use the default values included in the lead field. Finally, scalp data is generated: generate scalpdata simulates the components' signals, projects them according to their respective orientations through the lead field, and generates the data set as configured in epochs. In the last two lines of code, the data set is transformed into an EEGLAB-compatible format, and a ground-truth ICA decomposition is added.

When using the GUI, these steps would be followed in more or less the same order.

Sample Data Set

To demonstrate some of the toolbox's basic functionality, we now present a sample data set that was generated using SEREEGA and simulates a number of known cortical effects during a hypothetical experiment. We then apply established analysis methods to this data and present the results.

Simulated Experiment and Cortical Effects

The hypothetical experiment represents a visual go/no-go experiment. The visual stimulus in question is a pattern change, eliciting a N70-P100-N135 complex over occipital sites [START_REF] Pratt | Sensory ERP components[END_REF]. In the 'go' condition, the pattern change represents a target stimulus eliciting a strong P3a-P3b complex, most pronounced over the parietal-central-frontal midline [START_REF] Polich | Updating P300: An integrative theory of P3a and P3b[END_REF]. Upon perception of the target stimulus, the participants initiate a manual response with their right hand, resulting in mu-and beta band desynchronisation over sites covering the contralateral motor cortex.

Specifically, the visually evoked potential complex was modelled after findings by Gomez Gonzalez, Clark, Fan, [START_REF] Gomez Gonzalez | Sources of attention-sensitive visual event-related potentials[END_REF] and [START_REF] Pratt | Sensory ERP components[END_REF]. The N70 was a single dipole projecting centrally from the visual cortex. The P100 and N135 projected bilaterally from Brodmann areas BA 17 and BA 18, respectively. The P3a-P3b was modelled after [START_REF] Polich | Updating P300: An integrative theory of P3a and P3b[END_REF] and [START_REF] Debener | What is novel in the novelty oddball paradigm? functional significance of the novelty P3 event-related potential as revealed by independent component analysis[END_REF] with dipole locations taken from the latter. The manual response, finally, mimics mu-and beta-band findings of manually responded visual targets by Makeig, Delorme, et al. (2004).

The visually evoked potentials were modelled as single-peak ERP classes centred around 70, 100, and 135 ms, respectively, with widths of 60, 60, and 100 ms, and amplitudes of 7.5, 7.5, and 10 µV, as well as amplitude slopes of -2, -2, and -3 to simulate a level of habituation [START_REF] Pratt | Sensory ERP components[END_REF].

The P3a is an ERP class with a single 400 ms wide peak at 300 ms of amplitude 10 µV, and an amplitude slope of -5, simulating fatigue [START_REF] Polich | Updating P300: An integrative theory of P3a and P3b[END_REF]. The P3b ERP consists of two peaks. One at 400 ms (width 600 ms, amplitude 7 µV, amplitude slope -5), and one at 500 ms (width 1000 ms, amplitude 2 µV, amplitude slope -1). This latter peak represents the longer-lasting effects of the P3. The stronger slope of the P3b as compared to the P3a is a known finding reflecting habituation to the stimulus [START_REF] Polich | Updating P300: An integrative theory of P3a and P3b[END_REF]. All parameters of the ERP classes were given a deviation of 20% of their value relative to stimulus onset.

The motor cortex is simulated using four components for each class: both left and right motor cortex components are included, each with mu and beta band effects. In the non-target condition all components simply emit resting state mu (here as 8-12 Hz) and beta (19-26 Hz) activity implemented using the ERSP class with no modulation, and an amplitude of 2 µV. In the target condition, the left motor cortex classes are given an inverse burst modulation centred around 650 ms (mu; width 600 ms, taper .5, relative amplitude .5) and 600 ms (beta; width 500 ms, taper .8, relative amplitude .5) respectively. ERSP parameters had a deviation of 10% relative to stimulus onset, with the exception of latency deviation, which was fixed at 100 ms.

To simulate background processes, brown noise was added to each component (as per [START_REF] Freeman | Combining fMRI with EEG and MEG in order to relate patterns of brain activity to cognition[END_REF], with an amplitude of 5 µV. Furthermore, brown noise was also projected from a random selection of source locations, such that the total number of components in the data set was 64, all at least 25 mm apart from each other.

A single 'base' participant was defined and simulated using the specified component locations and orientations. For 10 subsequent participants, each component's location was randomly chosen within 20 mm of the original locations. Dipole orientations were randomly varied to be within ±25% of the original value for each axis.

This data set was simulated using the New York Head Model and projected onto the S64 montage of 64 channels. For each condition (target versus non-target), 100 epochs were simulated of 2000 ms each with a 500 ms pre-stimulus period, at a sampling rate of 1000 Hz.

Analyses and Results

The generated data set was subjected to a number of standard analyses to demonstrate their applicability and to illustrate the data itself.

First, figure 5 shows a single-subject ERP image analysis at three midline electrodes of the 'base participant'. This clearly shows the N70-P100-N135 complex, pronounced most strongly at parietal-occipital sites, as well as the P3, with the P3a most prominent over fronto-central sites and the P3b more pronounced further centrally and parietally. Also visible is the epoch-toepoch variability of the ERP signals in time and in amplitude, as well as their spatial differences with peaks pronounced differently at different sites. The effects over time can also be seen, with amplitudes diminishing at different rates.

Second, a cluster-based group-level analysis was performed. Using DIPFIT [START_REF] Oostenveld | DIPFIT 2.x[END_REF], equivalent dipole locations were reconstructed from the data set's ICA decomposition. These locations as well as the independent components' calculated ERSPs were used for k-means clustering. Due to the variability of the data as well as the clustering algorithm, manual adjustment was still necessary to create one cluster containing all left motor cortex components. The calculated dipole locations are given in figure 6. As was defined, they vary around a given mean location. Figure 7 shows the mean calculated event-related spectral perturbation for this cluster for both the target (left) and the non-target (right) conditions. The target condition shows clear mu and beta effects similar to those found by Makeig, Delorme, et al. (2004, figure 8).

Finally, the data was subjected to an analysis of classifiability for a hypothetical BCI system to be able to detect the elevated P3 response in the target condition compared to the non-target condition, and, separately, the manual response.

The P3 classifier followed a windowed means approach [START_REF] Blankertz | Single-trial analysis and classification of ERP components -a tutorial[END_REF], using the mean amplitude of six consecutive 50 ms time windows starting 200 ms after stimulus onset as features. Before, the signal was bandpass-filtered between 0.1 and 15 Hz.

The manual response was classified using a logBP approach [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF], using as features the power between 7 and 27 Hz, 400 to 1000 ms after stimulus onset, focused on sites covering the motor cortex: FC1, FC2, FC3, FC4, FC5, FC6, C1, C2, C3, C4, C5, and C6.

Both classifiers were implemented using BCILAB [START_REF] Kothe | BCILAB: a platform for brain-computer interface development[END_REF]. A regularised linear discriminant analysis classifier was trained to separate classes [START_REF] Duda | Pattern Classification[END_REF]. A 5 × 5 nested cross-validation with margins of 5 was used to select the shrinkage regularisation This is an authors' version of: Krol, L. R., Pawlitzki, J., Lotte, F., Gramann, K., & Zander, T. O. (2018). SEREEGA: Simulating Event-Related EEG Activity. Journal of Neuroscience Methods, 309, 13-24. doi: 10.1016Methods, 309, 13-24. doi: 10. /j.jneumeth.2018.08.001.08.001

Discussion

In this paper we presented the core functionality of SEREEGA, a free and open-source MATLABbased toolbox to simulate event-related EEG activity. SEREEGA allows any number of components to be defined in a virtual brain model, each with a specific location in the brain, a freely oriented scalp projection, and any number of signals. Five different types of activity are currently supported, allowing systematic effects in both the time and frequency domain to be simulated, mimicking known EEG features. Simulated EEG data can serve as ground truth for the development of new EEG analysis methods or the validation of existing ones.

The presented data set serves as a demonstration of some of SEREEGA's functions: a set of known cortical effects was simulated, mimicking a hypothetical experiment. It furthermore demonstrates how a basic initial configuration can be procedurally varied to simulate any number of additional participants with different anatomical layouts and variations in their cortical responses, leading to different results for example with respect to classification accuracies.

The classification accuracies presented here were only around 7 to 12 percentage points above significance on average. This is primarily due to the relative amplitude of the added noise, i.e., the signal-to-noise ratio of the generated data set. This can easily be varied to increase or decrease the classification accuracies as needed, or to investigate the accuracy as a function of this ratio. SEREEGA allows signal and noise components to be simulated separately, and then mixed to achieve a specific signal-to-noise ratio.

The ability to generate data as needed should be met with some caution. When validating signal detection mechanisms, for example, insufficiently considered parameters may lead to trivial results. Similarly, there is a potential danger of making circular arguments. For example, when SEREEGA is used to model a linear effect, and a linear method is subsequently used to again identify this effect, the obtained results could be essentially tautological. In such a case, the results would be applicable only to the extent that real EEG data contains a comparable linear effect. More generally, properties of signal simulation methods are based on assumptions that must not necessarily, or not fully, be met in real EEG data, and this limits the transfer of findings from simulated data to real recordings. This is a fundamental property of all simulation approaches, and users of simulations should be aware of the generator algorithms underlying the data, how they relate to the methods they are investigating, and how they may or may not relate to real EEG recordings.

Indeed, the open architecture of any general-purpose toolbox places some responsibility with the user to make use of the tools in a way that is compatible with the tested hypothesis. Different analysis methods have different assumptions as to e.g. the number of active sources or the statistical properties of the source signals. While SEREEGA is capable of adhering to a large variety of assumptions, not all data that is produced by it will adhere to the same ones. This depends on the design of the data.

It is thus highly important to carefully consider a simulation's parameters. With SEREEGA, different levels of complexity and realism can be generated. For some purposes, there may be no explicit need to simulate complex data sets or otherwise precisely model specific known cortical effects: The results of the methods may hold regardless of the complexity of the simulation. For example, to validate, in theory, a method such as logBP, a pure sine wave among white noise in a low number of randomly mixed channels may suffice. However, other methods do require a higher level of realism with e.g. a higher number of channels and realistic spatial dependencies, a higher number of noise sources, et cetera. For the sake of demonstration, we opted for a more realistic hypothetical experiment to demonstrate in this paper. We do not mean to say however that this level of realism is the explicit goal of SEREEGA, nor that it is methodically preferable. We did intend to illustrates how SEREEGA can be used for this purpose.

It should also be mentioned that realism here refers to the ultimately simulated source and scalp signals, not their underlying neurophysiology. As such, for example, ERPs are simulated directly as ERP-like time series, not as e.g. lower-level phase-resetting oscillations. Similarly, other known EEG effects such as phase-amplitude coupling, travelling waves, EEG microstates, and certain cross-trial dependencies are also not currently implemented. Additional signal classes or component properties would have to be written before testing methods that require such types of signal generators. The online documentation provides a tutorial.

Another important factor in a simulation's level of realism is how realistic the noise is. The included stochastic noise processes follow current simulation standards, where in particular 1/f and 1/f 2 frequency scales have been identified to characterise EEG [START_REF] Bédard | Does the 1/f frequency scaling of brain signals reflect self-organized critical states?[END_REF][START_REF] Freeman | Combining fMRI with EEG and MEG in order to relate patterns of brain activity to cognition[END_REF]. This is most likely a simplification. Empirical EEG recordings may have different spectral properties, especially non-stationarities over time, and are influenced by deeper and smaller processes whose activity is not directly measurable using EEG. However, there is currently no consensus regarding an alternative noise simulation method. In SEREEGA, any time series can be included in a simulation using the data class. This also allows noise extracted from empirical EEG recordings to be used. Where available, this may be a better solution. In the future, we may consider a public database of such recordings for all users to draw from. In fact, such a database could contain empirical recordings reflecting known cortical responses as well, allowing simulations to be pieced together by combining different parts of real EEG data. For this to be feasible, a large number of empirical data sets must be made available using a single, standardised format to describe them. Hierarchical Event Descriptor (HED) tags [START_REF] Bigdely-Shamlo | Hierarchical event descriptors (HED): Semi-structured tagging for real-world events in large-scale eeg[END_REF] offer a way to do that, and the Human Electrophysiology, Anatomic Data and Integrated Tools (HeadIT) Resource is one such database that could serve as precursor to the database needed for this purpose (Swartz Center for Computational Neuroscience (SCCN), n.d.).

Notable sources of artefacts are the eyes and muscles. The head models currently included in SEREEGA are cortical models and do not contain sources in the eyes or muscles. To the best of the authors' knowledge, no head model is currently available that does. A project is currently underway to construct a head model that contains a larger part of the human anatomy. For the moment, SEREEGA allows custom sources to be added to a simulation, by manually providing their coordinates and projection patterns. These can be taken from empirical recordings.

When taking data from empirical recordings, potential anatomical differences should be taken into account. Since SEREEGA supports different head models, it is not necessarily the case that they use the same coordinate format (e.g. Talairach, MNI, or custom measurements). Sources found at the same [x, y, z] location in different empirical studies or SEREEGA scripts may correspond to different actual locations depending on the head model used. Although authors who publish head models and lead fields often take great care to standardise their formats, they are not all uniform. Sometimes this is necessary: the dimensions of the Pediatric Head Atlas covering zero-to two-year-olds are of course smaller than those of an adult human's head model. SEREEGA includes plotting functions to investigate the size, shape, and location of a lead field's head model and the sources contained within, as well as to compare them against a standard adult head.

To the extent that their details were sufficiently reported, all (software) simulations that were performed by the authors mentioned in the introduction of this paper can be done using SEREEGA. The five signal types and other functions of SEREEGA cover and extend the processes used in those methods. This refers to the core steps of each simulation-further implementation details may of course differ. SIFT, for example, being focused on connectivity, has more options to control the parameters of the autoregressive models than SEREEGA currently does.

A limitation in the current architecture is that the defined components are necessarily independent at runtime: the procedurally generated activity of one component cannot presently depend on the activity in another. For such interactions to be modelled, the signals must be generated beforehand and included as data classes, as is the case for the autoregressive models.

In its current version, SEREEGA makes five accurate, high-density head models accessible in one toolbox, and allows five types of activation signals, enabling coloured noise, ERP-like deflections in the time domain, oscillatory activity with specific spectral characteristics and/or temporal modulations, and autoregressive signals to be simulated, as well as the inclusion of any other pre-generated or pre-recorded time series. This covers and extends the vast majority of past and present-day EEG simulation approaches. Additionally, FieldTrip-generated lead fields based on standard or custom head models can be used, and the toolbox's architecture allows it to be readily extended with additional head models and signals. For more specific needs, SEREEGA has been designed to be modular with extendability in mind.

SEREEGA thus represents, for the first time, an extensive general-purpose toolbox dedicated to simulating event-related EEG activity, making simulation methods more accessible, standardised, and reproducible.

With SEREEGA, we hope to help developers of EEG-based methods by making data simulation more accessible. With the consistent and increasing popularity of EEG, there is an accompanying need to further develop and validate EEG analysis methods. Simulated data can help with that by providing a ground truth to verify these methods' results. Seeing the current scope of application of simulated EEG studies, we believe that SEREEGA can help with the great majority of these, and can be easily extended for most of the remaining needs.

Figure 1 :

 1 Figure 1: Impression of the SEREEGA GUI.

Figure 2 :

 2 Figure 2: Channel locations and estimated brain boundaries for three head models. Left: The New York Head with its 228 channels. Middle: a 339-channel layout around a standard average head model available by default in FieldTrip. Right: The Pediatric Head Atlas (8 to 18, version 2) showing 2562 channels.

Figure 3 :

 3 Figure 3: Left: The location of source number 7479 in the New York Head as plotted onto EEGLAB's standard head model. Right: Four projection patterns of that source onto the scalp.From top left to bottom right: projection in the x-positive, y-positive, and z-positive direction, and a linear combination of the three (arbitrary units). This effectively orients the dipole at the source's location in 3D space.

Figure 4 :

 4 Figure 4: Sample signals annotated with the base variables that control their behaviour. From top left to bottom right: event-related potential (ERP); coloured noise; single-frequency unmodulated event-related spectral perturbation (ERSP); frequency band unmodulated ERSP; inverse burst-modulated ERSP; amplitude-modulated ERSP. Each base variable can additionally be given a deviation and a slope to control the epoch-to-epoch and systematic progressive variability of the signal. Not illustrated: the autoregressive model (ARM) and data classes.

Figure 5 :

 5 Figure 5: ERP image (Makeig, Debener, et al., 2004) at channel Fz, Cz, and POz showing the epoch-to-epoch variability of the simulated ERP signals as well as their spatial differences and effects over time.

Figure 6 :

 6 Figure 6: Locations of the left motor cortex equivalent dipoles in each simulated participant, from DIPFIT (Oostenveld & Delorme, 2003) calculations based on the independent components' scalp maps.

Figure 7 :

 7 Figure 7: Mean event-related spectral perturbation calculated for the left motor cortex cluster. Left: target condition (right-manual response). Right: non-target condition (no response).

Table 1 :

 1 Notation used in Section 2.1.

	Projection tensor of all sources

Table 2 :

 2 . As per the Journal's requirements, this version is licensed under CC BY-NC-ND 4.0.parameter, and to generate the estimates of the classifiers' accuracy listed in table 3.2. With chance level at 50% and significance reached at 57%[START_REF] Müller-Putz | Better than random? A closer look on BCI results[END_REF], we see that both classifiers reach significance for all 'participants', with the exception of participant 5, where no significant difference could be detected between the elicited P3 responses of the two conditions.This shows the data's variability on another level and demonstrates that levels of significance can be achieved with neither a ceiling nor a floor effect. Classification accuracies for two classifiers, one using the windowed means approach (WM) and one based on band power (BP). TP = true positive rate, TN = true negative rate, Acc = overall accuracy. Significance is reached at an accuracy of 0.57. The WM classifier for participant 5 does not reach significance.

	Participant WM TP WM TN WM Acc BP TP BP TN BP Acc
	0	0.78	0.75	0.76	0.70	0.64	0.67
	1	0.76	0.69	0.72	0.67	0.69	0.68
	2	0.74	0.74	0.74	0.57	0.58	0.57
	3	0.77	0.72	0.74	0.61	0.62	0.61
	4	0.72	0.57	0.64	0.65	0.64	0.64
	5	0.55	0.58	0.56	0.66	0.69	0.67
	6	0.68	0.57	0.62	0.61	0.70	0.65
	7	0.71	0.64	0.67	0.70	0.68	0.69
	8	0.71	0.61	0.66	0.63	0.55	0.59
	9	0.76	0.67	0.71	0.56	0.59	0.57
	10	0.67	0.57	0.62	0.69	0.64	0.66
	Mean	0.72	0.65	0.69	0.64	0.64	0.64

https://github.com/lrkrol/SEREEGA

This is an authors' version of: Krol, L. R., Pawlitzki, J.,Lotte, F., Gramann, K., & Zander, T. O. (2018). SEREEGA: Simulating Event-Related EEG Activity. Journal of NeuroscienceMethods, 309, 13-24. doi: 10.1016Methods, 309, 13-24. doi: 10. /j.jneumeth.2018.08.001.08.001. As per the Journal's requirements, this version is licensed under CC BY-NC-ND 4.0.

Acknowledgements

Part of this work was supported by the Deutsche Forschungsgemeinschaft (grant number ZA 821/3-1), Idex Bordeaux and LabEX CPU/SysNum, and the European Research Council with the BrainConquest project (grant number ERC-2016-STG-714567). The authors thank Mahta Mousavi for the band-pass filter design, Stefan Haufe for the autoregressive signal generation code and accompanying support, and Ramón Martínez-Cancino for his assistance with the GUI development.

Conflict of Interest

The authors declare no conflict of interest.