Johann Francou
email: yvon.francou@univ-reunion.fr

Cédric Abbezzot

Paulisimone Rasoavonjy

Didier Calogine

An Islanded Microgrid Design : A Case Study

Keywords: Hybrid power plant, solar power, conventional energy sources, batteries

Mafate is a natural protected cirque in Reunion Island, where about 300 dwellings are settled. The locals suffer from a massive energy insecurity because Mafate is not connected to the main grid. Although stand-alone solar microgrid installations are the preferred solution, they cause significant costs. "Micro Réseau Mafate" is a research project funded mainly by the European Regional Development Fund aiming at proposing strategies for microgrid development. In this paper, an actual islanded microgrid in Mafate is under scrutiny. The objective of the study is to rightsize the actual microgrid design considering the actual energy demand and solar resource measured in-situ. A Levelized Cost of Energy analysis is carried out and energy demand increasing scenarios are considered.

I. INTRODUCTION

Global electricity consumption has drastically increased over the years. However, some areas are still highly affected by energy insecurity for various reasons such as topography, local authority regulations or remoteness. Each of these 3 common causes generates additional costs for electrification. Solar based microgrids are suitable for rural electrification, as primary energy transportation is costless. However, rural microgrid development requests a major capital investment, especially for a reliable energy storage system. To secure power supply, particularly during bad weather periods, microgrid designers tend to oversize installations. Although microgrid oversizing is a security guarantee, it leads to important extra costs. The case study is described in Calogine et al. [START_REF] Calogine | A fractional derivative approach to modelling a smart grid-off cluster of houses in an isolated area[END_REF], it is a 7 kWp solar microgrid supplying 3 actual residential houses in Mafate of Reunion Island (France). The dwellings seen in figure 1 are located at about 50 meters away from the power plant. A 140 kWh lead acid battery capacity is installed in-situ. The microgrid, designed by the SIDELEC [2], is operational since january 2019 and tested in the scope of the European project entitled "Micro Réseau Mafate", aiming at proposing strategies for microgrids development. The installation is fully monitored to obtain weather, solar production, energy storage and energy consumption data. A microgrid model is used to simulate various configurations.

To foresee an increase of the actual consumption, dummy loads are generated by adding individual appliance profiles to the measured loads. The sub-metered profiles are drawn from public databases and aggregated according to probabilistic models. Hybridizing the solar microgrid with a diesel generator (DG) is also under examination to secure the power supply and to lower initial costs. As we have a 3 years feedback on the microgrid operation through the collected data, a rightsizing is proposed in this work, to determine a more profitable power plant configuration. An analytical method [START_REF] Yang | Battery Energy Storage System Size Determination in Renewable Energy Systems: A Review[END_REF] is applied to optimize the sizing against Levelized Cost of Energy (LCOE) and fuel consumption criteria. The process is carried out considering energy demand increase scenarios, and realistic configurations of the case study. The paper is organized as follows, section II describes the actual microgrid in Mafate, section III presents the microgrid model and section IV the load increasing scenarios, finally section V contains the sizing results and discussions.

II. CASE STUDY MICROGRID

In this section, a feedback from the case study microgrid operation is given through the measured data. Actual aggregated energy demand and solar resources in-situ are described.

A. Load evaluation

A majority of Mafate inhabitants has a very rural lifestyle, with a low access to energy. On the case study site, before microgrids were installed, they used weak solar installation and DG requiring a regular fuel supply. Storing a high fuel quantity in Mafate is not possible to avoid fire incident, as it is an UNESCO protected site. Mafate is not accessible on foot, each significant supplying has to be made by helicopter, leading to extra costs. Furthermore, wood fire cooking is generally favoured by the locals. Consequently, the inhabitants are naturally accustomed to saving electrical energy. This behaviour is demonstrated by the energy demand data, as seen in figure 2. Indeed, the 3 dwellings consume far less energy, with an average difference of 7.80 kWh between sizing demand hypothesis and actual energy demand. Thus, the microgrid is at least 58.00% too oversized. However, oversizing is necessary to mitigate the uncertainties on demand increases. To evaluate the load increase, the energy consumption is determined on a quarterly basis. Figure 3 demonstrates a 24% load increase between the fourth quarter of 2020 and 2021.

B. Solar resource evaluation

Another source of uncertainty is the solar resource. Daily insolation is represented in figure 4, the seasonal cycle can be clearly seen. June to July is the peak winter period while January to December is the peak austral summer. During cloudy days and rainy days, insolation values decrease sharply reaching 2 kWh/m²/day. Even though low insolation days is at a first sight dispersed over time, in the beginning of February 2022 we had 3 very cloudy days (from the 3rd to the 5th of February) in a row during Batsirai cyclone. At the end of February, a second cyclone named Emnati struck Reunion Island, responsible for the three consecutive very low insolation days (from the 20th to the 22nd of February). Those two extreme climatic events are important to be considered for the microgrid rightsizing process.

III. A MICROGRID MODEL

In this section, a microgrid model is presented. The model inputs are the weather data measured in-situ (irradiance and ambient temperature), and the real/dummy energy demand.

A. Microgrid architecture

The microgrid architecture is presented in figure 5. The energy demand is met by prioritizing the energy from the photovoltaic (PV) field. If there is a PV production surplus the battery is charged, while if there is not enough solar production, the battery is drawn. The DG works as a backup, whenever the battery is depleted.

Designers tend to limit the use of the DG because Mafate is a protected site where it is substantial to limit emitted pollution (noise pollution and greenhouse gas (GHG) emissions).

B. PV model

The PV model used is presented in Guezgouz et al. [START_REF] Guezgouz | Optimal Hybrid Pumped Hydro-Battery Storage Scheme for off-Grid Renewable Energy Systems[END_REF], the paramaters are adjusted according to the parameters of the PV panel used. PV output power is defined as :

E pv (t) = ∆t.N pv .P pv .f pv . G(t) 1000 .(1 + K t (T c (t) -T ref )), (1)
where ∆t, f pv , G, P pv , N pv , K t are respectively the time step, the losses factor, the solar radiation, the PV panel power peak, the number of panels composing the solar field, and the temperature coefficient of the PV panel. T c is the cell temperature, calculated as follow :

T c (t) = T em (t) + N OCT -20 800 .G(t), (2)
T em being the ambient temperature and N OCT the nominal operating cell temperature.

C. Battery model

In this work an energy balanced battery model is used. It has been demonstrated in Zhang et al. [START_REF] Zhang | The Influence of Photovoltaic Models and Battery Models in System Simulation and Optimization[END_REF], simplified battery model accuracy is not too far from a more detailed model. The battery model output are the State of Charge (SOC) and the depth of discharge (DOD) defined as :

SOC(t) = E b (t) E bnom × 100, (3)
DOD(t) = 100 -SOC(t), (4)
E b (t), E bnom , respectively the energy capacity of the battery at time t and its nominal energy capacity. E b (t) is calculated according to the previous energy level with

E b (t) = E b (t -1) + (E pv (t) - E load (t) η v ).η, (5)
η = η c in charging mode 1 η d in discharge , (6)
η v being the inverter efficiency and E load (t) the energy demand. For standalone microgrid design, battery lifetime is important to consider as it is generally the most expensive component. For our study, the simulated DOD for an average operating day is determined to estimate the battery lifetime with respect to the figure 6.

IV. LOAD INCREASE SCENARIOS

Each Mafate dwelling connected to the case study microgrid has a far lower energy demand than the average dwelling in Reunion Island due to their rural lifestyle. However, the figure 3 hints an eventual load augmentation for the future. There are several possible causes for energy demand increasing, two reasons are considered, load augmentation due to the use of more electrical appliances and load increasing due to connection of more dwellings to the microgrid. We note L k j Fig. 6. Cycle number as function of the DOD at 25°C. The figure is given by the OPzS battery supplier in technical datasheets [START_REF]Charging of HOPPECKE OPzS solar.power battery in Solar Applications[END_REF].

the load profile with j additional appliances and k connected dwellings. The applied load increasing process is defined as :

L m n = n j=0 a j + m k=0 d k , (m, n) ∈ N 2 , ( 7
)
a j is the j th appliance sub-meter power profile, a 0 is the always-off signal, i.e a constant function equal to 0 W. d k is the individual profile of the k th dwelling, d 0 is the alwaysoff signal. We note n and m respectively the number of added appliances and the number of dwellings connected to the power plant.

A. Additional appliances

With the microgrid, dwellers are more confident to consume energy, thus a likely scenario would be that the inhabitants will purchase and use more appliance types to ease their daily life. Between 2019 and 2021, dwellers used mainly fridge, freezer, television set, phone charger and lighting. For simulations, sub-meter load profiles are added incrementally to the aggregated load to evaluate the microgrid capabilities to meet higher loads due to the use of more electrical appliances. To construct the a j time-series in equation ( 7), an operating profile of the appliance to add is extracted from the Dataport public dataset [START_REF] Parson | Dataport and NILMTK: A Building Data Set Designed for Non-Intrusive Load Monitoring[END_REF]. Chunk extracted are repeatedly added to an "always-off" signal, a 0 , according to a realistic normal distribution of the hours of use for each specific appliance.

B. Additional dwellers

For profitability purpose, the microgrid operators intend to connect more dwellings to the power plant. Up to m = 3 in equation 7, the measured data are considered. To create a fictive load for an additional dweller, i.e when m > 3, an actual measured individual load, randomly chosen between the 3 monitored dwellers, is shuffled by day. The day shuffling allows to avoid biased power peak after aggregation. Choosing a load base between the 3 monitored dwellings allows to retain local consumption behaviours.

V. EXPERIMENTATION AND RESULTS

A. Microgrid configurations

For the present work, 4 configurations are under examination :

• configuration 1 (C1) : 7 kWp solar field, and a 140 kWh energy storage, • configuration 2 (C2) : 7 kWp solar field, and a 70 kWh energy storage, • configuration 3 (C3) : 3.5 kWp solar field, and a 140 kWh energy storage, • configuration 4 (C4) : 3.5 kWp solar field, and a 70 kWh energy storage. For each configuration, a 3 kVA DG is used if backup is required. The solar field for the initial configuration, C1, is composed of two solar panel arrays, each array has its own MPPT charge regulator. Therefore, halving the total solar field is equivalent to stopping one of the arrays. A same logic is applied to the batteries, as in C1 there are two battery banks connected in parallel where each battery bank is composed of 24 serially mounted cells of 2 V each. Consequently, cutting one of the battery fleet would keep the DC bus voltage at 48 V. Adding more solar panels and more batteries to the C1 configuration are not viable options as it would require more lands to be purchased.

B. Simulations

LCOE has been used and calculated in several works on microgrid sizing, such in [START_REF] Pawel | The Cost of Storage -How to Calculate the Levelized Cost of Stored Energy (LCOE) and Applications to Renewable Energy Generation[END_REF] and [START_REF] Lotfi | Levelized Cost of Energy Calculations for Microgrids[END_REF]. This economic metric is useful to compare energy system configurations based on each subsystem features. We use the LCOE for configuration comparison. We suppose a 25 lifetime years for the PV panels, the battery lifetime depends on its DOD as indicated in figure 6. Project lifetime is supposed to be 25 years for the simulation. A DG is used only when the battery is depleted. When backup is required, we suppose a 30% higher maintenance cost. For the DG operation, the fuel price is supposed to be slightly higher than usual prices due to helicopter deliveries. Simulations are processed on MATLAB/SIMULINK. The 4 different microgrid configurations are simulated according to 11 load scenarios, see Table II.

C. Results

The actual configuration C1 with the actual load, L 3 0 , has a very high LCOE, between 0.12 and 0.16 C/kWh more than the actual energy price in Reunion Island, which is considerable. Simulation and real operation of the microgrid show C1 with L 3 0 load does not need DG backup, see figure 8. Figure 7 demonstrates the LCOE results. From scenario 0 to scenario 4, for which only appliance profiles are added to the aggregated load, we see a slight decrease for each configuration. From scenario 5, new dwellings are connected to the power plant making the decrease more substantial. Indeed, increasing the load demand, thus increasing the produced energy, reduces the LCOE value. For grid operators to aim at a profitability, they have to connect new end-users to the actual microgrid. From scenario 6, with load L 5 0 , the energy cost becomes interesting regardless of the configuration. However, increasing the energy demand weakens the power plant. As seen on figure 8, from load scenario 9, C1 needs an auxiliary DG. C1 and C2 keep reasonable DG energy needs, even with 10 dwellings, L 10 0 , where about respectively 60 and 110 L of fuel per year are required. For those last configurations, the DG is activated only for the lasting bad weather periods. Hence, C1 and C2 does not need a too frequent fuel deliveries. After scenario 6, fossil fuel requirements for C3 and C4 increase drastically, contrary to C2 and C1. The last two have a higher solar energy collection capability, which is very important for the case of Mafate. Indeed, the afternoons are often cloudy, so the most of the solar energy collection has to be made in a short period in the morning otherwise the battery is deeply discharged and slowly charged, shortening its lifespan. That makes C3 and C4 less resilient to load increasing. C3 is then the most unfavourable configuration as it has a higher LCOE than C4 with an equal conventional source need. C2 gives an interesting trade-off between LCOE and fuel requirements. Up to 6 end-users, the C2 power plant is capable to operate without a DG at 0.11 C/kWh cost, slightly under the local average energy price. C2 can reach an even more interesting LCOE value (0.07 C/kWh) at a low fuel consumption (110 L/year) for 10 dwellers connected (scenario 11).

VI. CONCLUSION

The present work broaches the case of Mafate. The high surrounding relief makes it an off-grid area. To meet the local energy demand, solar microgrids are designed. Our case study is on a microgrid in Mafate powering 3 actual residential buildings. In this paper, the case study configuration is rightsized to reduce the LCOE and the DG use. Simulation results demonstrate more dwellings have to be connected to the installation to reach a more profitable microgrid operation. Batteries are the more expensive among the subsystems, halving battery capacity allows an average 9.00% LCOE decreasing, with an interesting resilience to load increasing, it can support until 6 dwellings connected without any backup needed. Our results show that hybrid systems are more resilient, with the use of a backup generator more than 10 dwellings can be powered with an affordable energy cost. Although fossil fuels have high impact on the environment, an economical and limited usage of conventional sources decreases energy storage investment costs, thus contributes to limit GHG emissions from electrochemical battery manufacture and transportation. For future investigations, it is worth wondering to what extent this GHG reduction compensates for the environmental impact of the DG operation.

Fig. 1 .

 1 Fig. 1. Drone photography of the case study microgrid.

Fig. 2 .

 2 Fig. 2. Daily aggregated energy demand, in-situ measurement compared to sizing hypothesis.

Fig. 3 .

 3 Fig. 3. Energy demand evolution on a quarterly basis. The bottom dashed line is the energy demand for the fourth quarter of 2020, while the top dashed line represents the energy demand for the fourth quarter of 2021.

Fig. 4 .

 4 Fig. 4. Daily insolation received on a horizontal surface at the experimentation site in Mafate. The dashed line represents the threshold for considering a very low insolation level.

Fig. 5 .

 5 Fig. 5. The case study microgrid architecture.

6thFig. 7 .

 7 Fig. 7. LCOE values for each load scenario. Energy prices are obtained from the local grid operator [3].

Fig. 8 .

 8 Fig. 8. Energy to be supplied by the DG for the 25 years simulation.

TABLE I LIST

 I OF THE ADDED APPLIANCES WITH THEIR NORMAL DISTRIBUTION DEFINITIONS, N (µ, σ 2 ) IS THE NORMAL DISTRIBUTION DEFINITION OF THE APPLIANCE WHERE µ AND σ ARE RESPECTIVELY THE MEAN HOUR OF USE AND THE STANDARD DEVIATION. THE OCCURRENCE IS THE NUMBER TIMES OF USE OF EACH APPLIANCE BY THE 3 INITIAL DWELLINGS.

	Appliances	Mean power (W)	Occurrence	N (µ, σ 2 )
	Washing machine (wm)	136	370	(18, 4 2 )
	Electrical oven (ov)	1558	85	(18, 4 2 )
	Microwave (mw)	1007	300	(9, 3 2 )
	Dish washer (dw)		549	150	(21, 3 2 )
			TABLE II
		CONSIDERED LOAD SCENARIOS
	id	Load	Energy (kWh/year)	Additional appliances
	0 1 2 3 4 5 6 7 8 9 10 11	L 3 0 L 3 1 L 3 2 L 3 3 L 3 4 L 4 0 L 5 0 L 6 0 L 7 0 L 8 0 L 9 0 L 10 0	2067 2128 2221 2234 2359 2773 4045 5120 6253 7318 8484 9182		-wm wm -ov wm -ov -mw wm -ov -mw -dw -------

ACKNOWLEDGMENT

The authors would like to thank the European Regional Development Fund, the Reunion Council, and the French Government for the funding to undertake this work. The authors would like to thank the SIDELEC for their valuable contribution on the "Micro Réseau Mafate" project.


Jan Fousek

Giovanni Rabuffo

Kashyap Gudibanda

Hiba Sheheitli

Viktor Jirsa

Spase Petkoski

Symmetry breaking organizes the brain's resting state manifold

Spontaneously fluctuating brain activity patterns that emerge at rest have been linked to brain's health and cognition. Despite detailed descriptions of the spatio-temporal brain patterns, our understanding of their generative mechanism is still incomplete. Using a combination of computational modeling and dynamical systems analysis we provide a mechanistic description of the formation of a resting state manifold via the network connectivity. We demonstrate that the symmetry breaking by the connectivity creates a characteristic flow on the manifold, which produces the major data features across scales and imaging modalities. These include spontaneous high amplitude co-activations, neuronal cascades, spectral cortical gradients, multistability and characteristic functional connectivity dynamics. When aggregated across cortical hierarchies, these match the profiles from empirical data. The understanding of the brain's resting state manifold is fundamental for the construction of task-specific flows and manifolds used in theories of brain function such as predictive coding. In addition, it shifts the focus from the single recordings towards brain's capacity to generate certain dynamics characteristic of health and pathology.

Introduction

The human brain at rest exhibits remarkable richness of neural activity structured both in time and space. Early computational modeling studies explored how these spontaneous fluctuations are constrained and how their organisation is shaped by the anatomic connectivity [START_REF] Honey | Network structure of cerebral cortex shapes functional connectivity on multiple time scales[END_REF][START_REF] Ghosh | Noise during rest enables the exploration of the brain's dynamic repertoire[END_REF][START_REF] Deco | Resting brains never rest: computational insights into potential cognitive architectures[END_REF][START_REF] Cabral | Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: Models and mechanisms[END_REF] enabling to start disentangling the mechanisms of the resting state dynamics in silico. A substantial body of work has related the emergent activity patterns at rest to the brain functional networks involved in task conditions [START_REF] Gusnard | Searching for a baseline: functional imaging and the resting human brain[END_REF][START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF], and shown that the spatio-temporal variability of resting-state activity possesses functional significance [START_REF] Deco | Emerging concepts for the dynamical organization of resting-state activity in the brain[END_REF][START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF][START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF], relevance to cognitive task performance [START_REF] Shine | The Dynamics of Functional Brain Networks: Integrated Network States during Cognitive Task Performance[END_REF], consciousness levels [START_REF] Cavanna | Dynamic functional connectivity and brain metastability during altered states of consciousness[END_REF], changes during ageing [START_REF] Battaglia | Dynamic Functional Connectivity between order and randomness and its evolution across the human adult lifespan[END_REF][START_REF] Petkoski | White-matter degradation and dynamical compensation support age-related functional alterations in human brain[END_REF], mental disorders [START_REF] Braun | From Maps to Multi-dimensional Network Mechanisms of Mental Disorders[END_REF], and neurodegenerative diseases (e.g. Alzheimer's dementia; [START_REF] Jones | Non-stationarity in the "resting brain's" modular architecture[END_REF]). The structure of the resting state dynamics changes over time [START_REF] Zalesky | Timeresolved resting-state brain networks[END_REF] and is characterized by a range of properties such as metastability [START_REF] Baker | Fast transient networks in spontaneous human brain activity[END_REF][START_REF] Beim Graben | Metastable Resting State Brain Dynamics[END_REF], event-like coactivations [START_REF] Tagliazucchi | Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis[END_REF][START_REF] Liu | Co-activation patterns in resting-state fMRI signals[END_REF][START_REF] Esfahlani | High-amplitude cofluctuations in cortical activity drive functional connectivity[END_REF] and traveling waves [START_REF] Gu | Brain Activity Fluctuations Propagate as Waves Traversing the Cortical Hierarchy[END_REF]. However, our understanding of the mechanisms underlying these spatiotemporal patterns of the brain activity at rest is still incomplete [START_REF] Mišić | Network-Level Structure-Function Relationships in Human Neocortex[END_REF] and whole brain network models have a crucial role to play on that front [START_REF] Breakspear | Dynamic models of large-scale brain activity[END_REF].

There is general agreement that the resting brain operates near criticality [START_REF] O'byrne | How critical is brain criticality?[END_REF]. This is supported by a large range of analyses performed on simulated and empirical data using network based measures (functional connectivity, functional connectivity dynamics), information theoretical measures (entropy, ignition) and descriptions of spatiotemporal dynamics (avalanches, cascades).

Modeling efforts provide further evidence for the close relationship between the empirical data features and the properties of the structural network, local dynamics, coupling strength, neural gain [START_REF] Cabral | Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: Models and mechanisms[END_REF][START_REF] Petkoski | White-matter degradation and dynamical compensation support age-related functional alterations in human brain[END_REF][START_REF] Deco | Key role of coupling, delay, and noise in resting brain fluctuations[END_REF][START_REF] Deco | Dynamical consequences of regional heterogeneity in the brain's transcriptional landscape[END_REF][START_REF] Melozzi | Individual structural features constrain the mouse functional connectome[END_REF][START_REF] Shine | The dynamic basis of cognition: an integrative core under the control of the ascending neuromodulatory system en[END_REF][START_REF] Roberts | Metastable brain waves[END_REF][START_REF] Courtiol | Dynamical mechanisms of interictal resting-state functional connectivity in epilepsy[END_REF]. The resting state dynamics can then be understood as noise-driven fluctuations of brain activity, operating near criticality and constrained by the brain connectivity [START_REF] Ghosh | Noise during rest enables the exploration of the brain's dynamic repertoire[END_REF][START_REF] Hansen | Functional connectivity dynamics: modeling the switching behavior of the resting state[END_REF]. However, none of the above qualifies as a description of a mechanism. Descriptions of mechanisms require formulation in terms of causal activities of their constituent entities and render the end stage, in our context the resting state dynamics, intelligible by showing how it is produced [START_REF] Machamer | Thinking about mechanisms[END_REF]. To explain is thus not merely to redescribe one regularity (e.g. functional connectivity dynamics, or maximization of entropy) as a series of several (such as near-criticality, cascades, ignition). Rather, explanation involves revealing the productive relation between causal activities linked to their constituent entities.

In this paper we aim to remedy this situation and provide this explanation using Structured Flows on Manifolds (SFMs) [START_REF] Jirsa | Selbstorganisation-ein Paradigma für die Humanwissenschaften[END_REF][START_REF] Pillai | Symmetry Breaking in Space-Time Hierarchies Shapes Brain Dynamics and Behavior[END_REF][START_REF] Huys | Functional architectures and structured flows on manifolds: a dynamical framework for motor behavior[END_REF][START_REF] Mcintosh | The hidden repertoire of brain dynamics and dysfunction[END_REF][START_REF] Woodman | Emergent dynamics from spiking neuron networks through symmetry breaking of connectivity[END_REF]. SFMs is a mathematical framework explaining how low dimensional dynamics, reflecting generative sets of rules underlying behavior, emerges in high-dimensional nonlinear systems, specifically dynamical systems on networks modeling macroscale brain dynamics. When properly linked to the network's constituent entities (functional nodes and connectivity), we will demonstrate how their causal activities lead to the formation of brain's resting SFM, comprising all its dynamic signatures (see Figure 1). If we distill the previous reports of brain resting state data analysis from the dynamical systems point of view, we arrive at the following main empirical signatures that should be part of the end stage of a successful mechanistic description: bistability of single region activation [START_REF] Watanabe | A pairwise maximum entropy model accurately describes resting-state human brain networks[END_REF][START_REF] Watanabe | Energy landscapes of resting-state brain networks[END_REF][START_REF] Gu | The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure[END_REF], low-dimensionality of the global system dynamics in state space [START_REF] Deco | Emerging concepts for the dynamical organization of resting-state activity in the brain[END_REF][START_REF] Ashourvan | The energy landscape underpinning module dynamics in the human brain connectome[END_REF][START_REF] Vohryzek | Ghost Attractors in Spontaneous Brain Activity: Recurrent Excursions Into Functionally-Relevant BOLD Phase-Locking States[END_REF], cascade propagation [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF], multistability of recurrent coactivation spatial patterns [START_REF] Beim Graben | Metastable Resting State Brain Dynamics[END_REF][START_REF] Cornblath | Temporal sequences of brain activity at rest are constrained by white matter structure and modulated by cognitive demands[END_REF] and their nontrivial temporal dynamics or intermittency [START_REF] Esfahlani | High-amplitude cofluctuations in cortical activity drive functional connectivity[END_REF][START_REF] Hansen | Functional connectivity dynamics: modeling the switching behavior of the resting state[END_REF][START_REF] Wong | A recurrent network mechanism of time integration in perceptual decisions[END_REF]. These signatures will constitute the key features of what we will describe as structured flows on the low dimensional resting state manifold. Structured flows on manifolds as focus of resting state characterization. With respect to the structure of the connectivity of the dynamical system, we consider spectrum defined by the two symmetrical limit cases: fully connected and fully disconnected network. Driven by noise, the disconnected system exhibits fully statistical, high-dimensional dynamics -it explores the whole state space in a equidirectional manner. On the other hand, the dynamics of the fully connected system is fully constrained corresponding to a SO(n) hypersphere with zero flow. The dynamics on the sparsely connected system leads to an object in between -a low-dimensional attractive manifold with an associated flow (SFM). It is this object we wish to put in the center of interest and characterize. While the SFM object remains the same, connections are made to data of various modalities with the help of suitable data features.

Results

In what follows we employ whole-brain modeling to study the low dimensional manifold and the associated structured flows of the spontaneous resting state dynamics, and how these relate to the structural connectome. We constructed a brain network model (BNM) in the Virtual Brain [START_REF] Sanz-Leon | The Virtual Brain: a simulator of primate brain network dynamics[END_REF] using the two-dimensional mean-field model of an ensemble of quadratic integrate-and-fire neurons ( [START_REF] Montbrió | Macroscopic description for networks of spiking neurons[END_REF]; MPR) to govern the regional dynamics coupled with a connectome derived from a subject from the Human Connectome Project [START_REF] Van Essen | The WU-Minn human connectome project: an overview[END_REF]. We applied Balloon-Windkessel model [START_REF] Stephan | Comparing hemodynamic models with DCM[END_REF] to the simulated neuronal mass activity to generate realistic BOLD signals. From these, we computed the Dynamical Functional Connectivity (dFC) that captures the changes in the system's dynamics on the slow time scale, which we compared with empirical recordings. The fast neuronal activity is decomposed in a 2N-dimensional state space using Principal Component Analysis (PCA) to unveil the low-dimensional manifold on which the system evolves (see Materials and Methods for more details).

When driven by noise, the network of the bistable MPR nodes has the capacity to exhibit realistic dFC when the network input is scaled appropriately [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF]. The noise together with the network input drives the switching between up-and down-state of the individual nodes, while the network mediates the coordination reflected in the functional connectivity. In the following sections, we explore how the manifold of the resting state activity arises from the networked interactions, how it shapes the multistability of the functional connectivity in the simulated BOLD, and how it relates to empirical observations.

Symmetry breaking: working point for dFC

To assess the impact of the symmetry breaking by the connectome, we simulated 10 minutes of spontaneous activity for a range of values of the coupling scaling parameter G and noise variance σ, and applied PCA to the source signal Ψ(t) and dFC to the BOLD (Figure 2). We used the variance accounted for (VAF) of the first two PCA components as an estimate for the dimensionality of the system's dynamics in the state-space (Figure 2D), and the variance of the upper triangle of the dFC matrix as a measure of the fluidity of the system's dynamics-that is the propensity to dwell in specific brain states (defined by the functional connectivity) and shift and return between several such states (Fig- ure 2A). In addition, using Kolmogorov-Smirnov distance between the centered distributions of the values of the upper triangle of the dF C w in the empirical and simulated data, we have verified that the region of the parameter space where dFC is most similar to the one derived from empirical data overlaps with the region with the highest fluidity, Figure 2A. and noise variance σ to produce both time-series of the state space variables r(t), V(t), and the BOLD signal. For each combination of G and σ we compute the sliding window dF C w matrix from the simulated BOLD signal, and quantify the "switching index" of the dFC as the variance of the upper triangle (A). Kolmogorov-Smirnov distance of the centered (mean-subtracted) distributions of the values of the upper triangle of the dFC computed from empirical and simulated resting state BOLD time series. The region of parameter space where the distributions are closest overlaps with the region with high fluidity of dFC (B). For selected values of (G, σ) we show the sliding window dF C w (C), edge based dF C e (D) and the projection of r(t) time series in the first two PCA components (E) annotated with corresponding fractional variance accounted for (VAF). In the working point around G = 0.54 and intermediate values of σ the system exhibits recurrence in the large-scale dynamics as captured by nonzero switching index, and reduction of dimensionality as captured in the increase in explained variance by the first PCA components and the asymmetry in the respective projection. For values of G below or above the working point, the systems loses the fluidity property as reflected in the absence of the off-diagonal blocks on the dFC, and exhibits high-dimensional dynamics.

Network dynamics

Before we delve into the characterization of the low-dimensional manifold, let us first describe the network dynamics in detail. For the MPR model, the dynamical profile of an isolated node in the bistable parametrization consists of an unstable fixed point (saddle node) and two stable fixed points: down-state stable node and up-state focus (Figure 3A). Considering the uncoupled system, that is, the joint dynamics of the N populations (nodes) in the absence of any inter-population synaptic coupling, the phase flow is represented by 2 N stable fixed points that contain all possible combinations of the populations firing at either their low or high mean firing rates (down or up state, respectively). Starting from an initial condition, the system will settle into the nearest accessible such fixed point, a stable network state composed of a corresponding combination of regions in their up or down state. Thus, the dynamics of the uncoupled system in phase space can be thought of as being driven by a potential energy landscape with multiple stable local minima representing the stable attractor states of the network. In fact, the uncoupled system, as such, is invariant under permutation of the indexes of the populations, such that these latter attractor network states are distinguished only in terms of the respective number of nodes in up and down states. The global dynamics of the system, thus, collapses in finite time onto this stable attractor state composed of a finite set of stable equilibrium points that is invariant under shuffling of indexes of the nodes.

The associated global phase flow can be decomposed into N projections onto the identical 2D phase planes of individual populations, depicted in Figure 3A.

Viewed from this perspective, the structure of the basins of attraction of the 2 N stable system equilibrium points redundantly inherits, in higher dimensions, the relative structure of the basins of attraction of the two stable fixed points of an individual population. For an isolated node, varying the external input I i changes the size of the basins of attraction of the stable fixed points. This modulates the probability of switching between the two states when driven by noise as captured by the mean escape times (Figure 3A, see Methods for more detail). For a connected node, the external input I e depends on the state of the neighboring nodes (see Equation 4), fluctuating as they transition between the up-and down-state. On the network level, given right scaling of the network connections, this enables the cascades of up-and down-state switching at the fast time-scale, and the co-fluctuation of the BOLD signal (Figure 3B).

To understand better the dynamical underpinning of the increase of fluidity of the dFC we assess the characteristics of the co-fluctuations of the BOLD signal and the cascades in the source signal. For the co-activations, we start from the edge time series which is defined as pairwise dot product of z-scored BOLD signal (an average over the edge time series would correspond to the pearson correlation). The correlation across time-points yield the dF C e matrix capturing the recurrence of the edge configurations,and the root sum squared (RSS) over the edges at each time point captures the contribution of that particular time point to the overall functional connectivity (see Methods for more details). The time-points crossing the 95-th percentile threshold of the RSS are considered as strong co-activation events. The neuronal cascades [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF] are long lasting perturbations of the neuroelectric activity and are measured on a global level as a sum over regions of the binarized firing rate activity (at the threshold of 3 standard deviations). We compared these measures between the working point G w , the disconnected system G = 0, the strong network coupling regime G >> G w , and the empirical data (Figure 3C,D).

In the working point G w the co-activations include large number of edges (Figure 3D) and the RSS follows the number of cascades up to a short delay corresponding to the delay of the BOLD signal. Moreover, some of the strong co-activations re-occur partially in time as reflected in the non-zero elements of the dF C e matrix. The same profiles can be observed in the empirical data, namely in the simultaneous EEG and fMRI recordings. On the other hand, the characteristic spatial and temporal structure is lost outside of the working point, that is either for the weakly coupled system (G << G w ), or for too strong coupling (G >> G w ).

To quantify how the co-activation events contribute to the characteristic similarity across time, we compare the correlation of the edge vectors during the events, during the non-events, and between events and non-events. As a result we observe an increased similarity of the edge vectors during the events both in the empirical data and in the simulations in the working point G w . Again, this property is lost for too weak (G << G w ) or for too strong coupling (G >> G w ).

Together, these results show, that the system has a similar dynamical profile in the working point G w as observed in the empirical data with respect to the network-carried fluctuations on both the fast and slow timescales (as captured by dF C and cascades respectively). Example of a cascade-coordinated increase in activity translating to a delayed correlated peak in the BOLD signal. Below we compare the network dynamics in and outside the working point, and the empirical data. In both empirical data and the working point (G = 0.54), the BOLD co-activations follow the neuronal cascades with a latency (C), and show distinct spatial profiles which are recurrent in time (D): edge time series on the (top panel) captures the spatial profiles of the co-activations, the similarity across time is captured by the dF C e matrix (middle panel), and the distributions of correlation between co-activation events (CA) and non-events (nCA) is compared (bottom panel).

Manifold of the resting state and characteristic subspaces

Having characterized the dynamics of the system in the working point with an appropriate measure, we proceed with the description of the manifold on which it evolves-that is to show that this behavior is low-dimensional and constrained to a specific subspace. To relate cascades and co-activations to the trajectories of the system in the 2N-dimensional state space, we first select time intervals with similar functional connectivity. Starting from the edge time series for the magnitude of co-fluctuations, we clustered the time points using k-means (k=5). This separated the high-activity intervals (majority of the nodes in the up-state), low-activity intervals (majority of the nodes in the down-state), and the co-fluctuation events (Figure 4A).

Next, we identified the trajectories of the system underlying each cluster in the low dimensional projection of state space. For each cluster, we have selected the corresponding time points in the state space of the system, and projected them into the first two principal components of the PCA computed on the complete time series. We have observed that while the corresponding subspaces overlap partially in the projection (Figure 4B, colors correspond to the clusters), the activity within the clusters concentrates to different subspaces. This concentration in different subspaces is reflected in the distance between the centroids of the cluster time points in the PCA projection.

While the cluster activity overlaps in the projection in the components of the PCA computed from the whole time series, the co-activation trajectories become clearer by choosing different basis to span the low-dimensional space, that is, to compute the PCA from the time points corresponding to the individual co-activation events. To project the trajectories of the events observed at the slow time scale of the BOLD on the manifold, we have shifted the BOLD signal by the characteristic lag, and for each BOLD time point belonging to the cluster we selected the corresponding time points in r(t), and convolved the resulting data with a Gaussian kernel to smoothen out the noisy fluctuations (see Methods for details). We then spanned the subspace corresponding to the first two PCA components of the co-fluctuation trajectory and overlayed the smoothened trajectory over the density plot of the full r(t) time series. The density plots (shades of red, Figure 4C) of the example events show a separation of the event subspace marked by the peak in the RSS (shown in yellow in Figure 4C) on the smoothened trajectory from the rest of the manifold. This suggests that the event subspace is relatively stable, allowing the system to dwell in it long enough to cause the significant peaks in the slow BOLD signal, and that the intermediate states are less stable than the event subspace or the rest of the manifold and visited only transiently.

Although the linear embedding of the whole time-series does not separate the event trajectories well when applied to the r(t) time series, the event trajectories concentrate in the high-activity subspace spanned by the first two PCA components of the BOLD signal (Figure 4D) .

Together, these results chart the low-dimensional manifold of the system in the working point regime, associating the subspaces with specific flows. The fluid dynamics as characterized in the previous section then arise from the slow transitions between the low-and high-activity subspaces, where the latter supports the strong co-activation events which are reflected in the dF C.

Fixed point skeleton and structured flow

To understand how the resting state manifold arises, we start by considering the uncoupled system, that is, the joint dynamics of the N populations (nodes) in the absence of any inter-population synaptic coupling. This uncoupled system's phase flow is dominated by 2 N stable fixed points that represent all possible combinations of the populations firing at either their low or high mean firing rates (down or up state, respectively). Starting from an initial condition and in the absence of noise, the BNM will settle into the nearest accessible such fixed point, a stable network state composed of a corresponding combination of regions in their up or down state.

The dynamical effects of the symmetry breaking in the BNM are delineated by the topology of the connectome. The heterogeneity of the in-degree (total connectivity) of individual nodes of the network drives a variation in the relative positioning of the separatrices between the basins of attraction of the equilibrium points, mirrored in the variation of the corresponding projections onto the 2D phase planes of corresponding nodes (see Figure 3A). In conjunction, connectivity strength and topology give rise to gradients in the relative attractiveness of the system's equilibrium states. This attractiveness (or stability) can be quantified by the largest negative real eigenvalues obtained from the linearization of the system about the respective equilibrium state (linear stability analysis).

To map the complete manifold outside the simulated trajectories we sampled the stable fixed points for varying coupling scaling parameter G from the 2 N combinations of up-and down-states, and evaluated their stability (see Methods for more details). We found that the number of stable fixed points in the sample decreases with increasing G. This decrease is due to the loss of states with mixed composition of up-and down-state due to the bifurcation of the down state in nodes with high input (Figure 5A). Projecting the r component of the fixedpoints in the first two eigenvectors of the Laplacian confirms this thinning of the intermediate compositions biased towards those with higher number of nodes in the up-state (corresponds to the first Laplacian eigenvector λ 1 ). Additionally, the stability of the fixed points was inversely proportional to the number of nodes in the up-state, that is in the direction of E 1 the first eigenvector of the Laplacian (Figure 5B).

To put this into the context of the simulated trajectories, we have next identified the fixed points around which the simulated trajectory evolved by taking initial conditions from the simulated trajectory, and integrating the system without noise to the equilibrium. We have confirmed that in all instances the system reached a stable fixed point composed of combination of up-and down-states, and that the stability of these fixed points follows the same gradient in terms of the composition (Figure 5C).

Furthermore, the nodes of the network exhibit a frequency gradient of the oscillations in the up-state (Figure 5D). This gradient reflects variability of the characteristic frequency in the up-state across nodes in the network. In the fixed-point state, if the nodes are treated as isolated systems with an input current term based on the existing network state, then

r * i = r * + δ r i v * i = v * + δ v i (1)
where (r * , v * ) are the symmetric fixed-points of the network and (δ r i , δ v i ) are the excursions from the symmetric fixed-point and change according to the existing network state. These excursions depend directly on the in-strength of the ith node and the local states of its first neighbours.

Following linear stability analysis of the ith system around the fixed-point (see Methods), the eigenvalues of the Jacobian matrix are given by

λ 1,2 = 2v * i ± Jr * i -4π 2 r * 2 i = 2v * i ± 2Jr * i -4π2r * 2 i + 2Jδ r i -4π 2 δ r 2 i -8π 2 r * δ r i (2)
From the above equation, we see that the frequency of oscillations in the upstate of the ith node increases proportional to δ r i and therefore proportional to the in-strength of the node, which we also observe in the simulation (Figure 5D).

Furthermore, applying the PCA projection on the empirical BOLD timeseries, we have identified a similar separation of the event trajectories in the global embedding as observed in the simulations (Figure 5E). However in the case of the empirical data the system exhibits both the co-activations and codeactivations as seen on the separation through the first PCA component.

Lastly, symmetry breaking by the connectivity alone results in a spatial organization of the above described flow which is aligned with empirically observed trends (Figure 5F). In particular, per region the time spent in avalanches and the cumulative z-scored BOLD signal within events both decrease across the cortical hierarchy from the primary to paralimbic regions. The established principal functional gradient extracted from the empirical fMRI data is also aligned along this axis .

Discussion

Using a combination of computational modeling and dynamical systems analysis we have provided a complete mechanistic description in terms of constituent entities and their causal activities leading to spontaneous co-activations and neuronal cascades in the brain's resting state [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF]. We showed how the breaking of the symmetry of the BNM's connectivity gives rise to the structured lowdimensional dynamics in the phase space and recurrent fluctuations of the func-tional connectivity (Figure 2). These fluctuations arise from network-mediated cascades of up-and down-state switching and capture well the empirically found relationship between the strong co-activation events and the recurrence structure reflected by the functional connectivity dynamics (Figure 3). The subspace accessible to the brain in this regime was charted and partitioned according to the characteristic flow associated with each of the partition (Figure 4). Finally, this subspace and its associated flows arise from the rich fixed point structure of the system and the differential stability of the nodes in these fixed points is not only reflected in the propensity to state switching that reflects the cortical hierarchy, but also influences the dominant oscillation frequency (Figure 5).

In summary, these results support our hypothesis that the recurrent functional connectivity states of the resting state correspond to distinct subspaces on a low-dimensional manifold associated with distinct structured flows.

The central result of our work is that the symmetry breaking via the structural connectivity carves out an attractive subspace of all the possible states of the brain, and that the flow on this manifold governs the characteristic dynamics of the brain (that is discarding the transient towards the manifold from arbitrary initial conditions). In this regime, the model captures the multistability and noise-driven exploration of the dynamic repertoire explored previously in computational studies [START_REF] Ghosh | Noise during rest enables the exploration of the brain's dynamic repertoire[END_REF][START_REF] Courtiol | Dynamical mechanisms of interictal resting-state functional connectivity in epilepsy[END_REF][START_REF] Hansen | Functional connectivity dynamics: modeling the switching behavior of the resting state[END_REF][START_REF] Deco | Ongoing cortical activity at rest: criticality, multistability, and ghost attractors[END_REF][START_REF] Haimovici | Brain organization into resting state networks emerges at criticality on a model of the human connectome[END_REF]. The data features extracted from the time series provides a link between the empirical data and the model. Here, the functional structure in the brain is carried by the rare high amplitude cofluctuation events as it was previously demonstrated in empirical fMRI data [START_REF] Tagliazucchi | Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis[END_REF][START_REF] Esfahlani | High-amplitude cofluctuations in cortical activity drive functional connectivity[END_REF][START_REF] Sporns | Dynamic expression of brain functional systems disclosed by fine-scale analysis of edge time series[END_REF], and in simultaneous EEG and fMRI measurements [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF]. Similarly, recent modeling study has shown the role of structural modules of the network in shaping the co-fluctuation events [START_REF] Pope | Modular origins of high-amplitude cofluctuations in fine-scale functional connectivity dynamics[END_REF], which is aligned with the brain network as the symmetry breaking gives rise to the low-dimensional dynamics.

The slow time scale fluctuations of the dynamical functional connectivity reflect the movement of the brain activity between the low-and high-activity subspaces of the manifold. The flow in the high-activity subspace supports the cascades, which in turn are reflected in the high activity coactivations. This movement points to the multistable rather then metastable interpretation of the resting state dynamics [START_REF] Cocchi | Criticality in the brain: A synthesis of neurobiology, models and cognition[END_REF], and reflects the observation of switching between a low-amplitude incoherent and high-amplitude coherent states in empirical data [START_REF] Kong | Sensory-motor cortices shape functional connectivity dynamics in the human brain[END_REF]. Furthermore, the slow transitions between the high-and low-activity subspaces is compatible with the reports on the spontaneous infra-slow brain activity [START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF][START_REF] Mitra | Spontaneous Infra-slow Brain Activity Has Unique Spatiotemporal Dynamics and Laminar Structure[END_REF] and the detailed reports on its spatio-temporal structure. For example, the slow traveling waves [START_REF] Gu | Brain Activity Fluctuations Propagate as Waves Traversing the Cortical Hierarchy[END_REF] propagating along the principal gradient of cortical organization [START_REF] Margulies | Situating the default-mode network along a principal gradient of macroscale cortical organization[END_REF][START_REF] Huntenburg | Large-Scale Gradients in Human Cortical Organization[END_REF] would provide a refined description of the trajectory through the manifold subspaces.

The attractive subspace of the low-dimensional manifold and the associated structured flow arise in the presented system from the changes in the fixedpoint structure due to the irregular connectivity. In particular, the network input mediates the modulation of the escape times of the noise-induced transitions. These chain into domino-like sequences [START_REF] Ashwin | Fast and slow domino regimes in transient network dynamics[END_REF][START_REF] Ashwin | Sequential escapes: onset of slow domino regime via a saddle connection[END_REF], which in turn constitute the neuronal cascades. On a network level, our results elaborate the previous Overall, the movement of the system through the subspaces of a low-dimensional manifold is in accordance with empirical and modeling results on recurrence and state clustering of resting state fMRI BOLD recordings. Using clustering algorithms to partition the BOLD time series yields statistically similar and temporally recurrent whole brain spatial coactivation patterns [START_REF] Beim Graben | Metastable Resting State Brain Dynamics[END_REF][START_REF] Cornblath | Temporal sequences of brain activity at rest are constrained by white matter structure and modulated by cognitive demands[END_REF] associated with specific dwell times and transition probabilities. However, compared to the clustering approaches applied to the BOLD time-series, the SFMs allows us to refine the partitioning of the state-space in two aspects: we unfold the subspaces based on the similarity of the coactivations on the level of the BOLD signal, and we provide a detailed description of the flow of the system in these subspaces e.g. in terms of the cascades. The former is in line with the recent advances regarding the low-dimensional representation of meso- [START_REF] Kim | Cortical-hippocampal coupling during manifold exploration in motor cortex[END_REF][START_REF] Chaudhuri | The intrinsic attractor manifold and population dynamics of a canonical cognitive circuit across waking and sleep[END_REF] and macroscopic [START_REF] Favaretto | Subcortical-cortical dynamical states of the human brain and their breakdown in stroke[END_REF][START_REF] Rué-Queralt | Decoding brain states on the intrinsic manifold of human brain dynamics across wakefulness and sleep[END_REF] brain dynamics, but the latter describes the origin of those subspaces as constrained by the connectome. Interestingly, the clustering of phase-locking BOLD states [START_REF] Vohryzek | Ghost Attractors in Spontaneous Brain Activity: Recurrent Excursions Into Functionally-Relevant BOLD Phase-Locking States[END_REF] leads to very similar low-dimensional representation of the resting state dynamics to our approach, with a single dominant global phase locked state and a number of transient partially phase-locked states related to functional networks. Similarly, by embedding the resting state data onto the task manifold extracted with the help of diffusion maps, [START_REF] Gao | Nonlinear manifold learning in functional magnetic resonance imaging uncovers a low-dimensional space of brain dynamics[END_REF] found that resting state time-points concentrate in the task-fixation and transition subspaces, and only a minority of time-points reach to the cognitive subspaces of the task manifold.

The description of the structured flow addresses also the fast time-scale by including the cascades, which we previously showed to relate to the co-activations observed in the BOLD signal [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF]. In EEG literature, the spatio-temporal structure of the resting state dynamics is characterized with the help of microstatessensor-level transient patterns lasting on average for 60-150 ms [START_REF] Michel | EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review[END_REF]. Attempts have been made to relate the microstates to BOLD activation clusters [START_REF] Britz | BOLD correlates of EEG topography reveal rapid resting-state network dynamics[END_REF][START_REF] Bréchet | Capturing the spatiotemporal dynamics of self-generated, task-initiated thoughts with EEG and fMRI[END_REF],

but identifying the sources generating the microstates with clustering or regression analysis has been challenging so far due to unclear relationship between the broadband EEG activity and the BOLD signal fluctuations [START_REF] Drew | Ultra-slow oscillations in fMRI and resting-state connectivity: neuronal and vascular contributions and technical confounds[END_REF]. To advance we propose to reframe the question as a search for a shared manifold of the neuronal activity with specific slow and fast time-scale characteristics which in turn are reflected in the EEG and the BOLD observables.

The manifold we describe is conceptually reminiscent of energy landscapes described in previous works [START_REF] Watanabe | Energy landscapes of resting-state brain networks[END_REF][START_REF] Gu | The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure[END_REF]. However, previous energy landscape models, such as in [START_REF] Gu | The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure[END_REF], implicitly assume energy minimization and thus, by construction, encode the hypothesis that the activation of two brain regions that are connected via a direct structural connection is more energetically favorable than that of two regions that are not directly connected. We make no such assumption here and, instead, the effective energy landscape emerges, in the form of a low dimensional manifold, out of the interplay of the non-linearity in the local neural mass model and the connectome, thus fully embracing the the network impact, beyond the pair-wise interactions. In addition, previous energy landscape analysis [START_REF] Watanabe | Energy landscapes of resting-state brain networks[END_REF] assumed that the network changes only gradually by flipping one region at a time, and did not account for transitions in which several regions flip simultaneously. Treating the brain as a whole, the BNM that we presented here instead allows for such latter transitions of the system in state space, which may very well be due to strongly connected regions that are able to simultaneously influence their nearest neighbors during coactivation events.

It is worth pointing that our framework covers only one part of the mechanisms that shape the brain's manifold and the flow on it, that is the connectome.

We have assumed identical parameters for each region, ignoring the known structural hierarchies [START_REF] Wang | Macroscopic gradients of synaptic excitation and inhibition in the neocortex[END_REF], which have been shown to improve the predictive value of the BNMs [START_REF] Kong | Sensory-motor cortices shape functional connectivity dynamics in the human brain[END_REF][START_REF] Chaudhuri | A Large-Scale Circuit Mechanism for Hierarchical Dynamical Processing in the Primate Cortex[END_REF][START_REF] Wang | Inversion of a large-scale circuit model reveals a cortical hierarchy in the dynamic resting human brain[END_REF]. While we observed differential functional properties of the nodes across the cortical hierarchy [START_REF] Mesulam | From sensation to cognition[END_REF], we didn't recover the exact spatial correspondence to the established functional gradients [START_REF] Margulies | Situating the default-mode network along a principal gradient of macroscale cortical organization[END_REF]. Neuromodulation and the subcortical drives [START_REF] Shine | The thalamus integrates the macrosystems of the brain to facilitate complex, adaptive brain network dynamics[END_REF] are another missing aspect that similarly improve the performance of BNMs [START_REF] Kringelbach | Dynamic coupling of whole-brain neuronal and neurotransmitter systems[END_REF]. However, both of these elements are not yet established in the framework of BNMs, as is the impact of the connectome [START_REF] Breakspear | Dynamic models of large-scale brain activity[END_REF].

Thus our goal here is not to generate in silico observables that are as close as possible to the empirical one, which nevertheless differ a lot depending on the preprocessing, e.g. see [START_REF] Aquino | On the intersection between data quality and dynamical modelling of large-scale fMRI signals[END_REF], but to focus on the generative mechanisms for the key data features across time-scales and neuroimaging modalities that render functional activity identifiable across subjects [START_REF] Battaglia | Dynamic Functional Connectivity between order and randomness and its evolution across the human adult lifespan[END_REF][START_REF] Petkoski | White-matter degradation and dynamical compensation support age-related functional alterations in human brain[END_REF][START_REF] Van De Ville | When makes you unique: temporality of the human brain fingerprint[END_REF].

A natural next step will be to extend the analysis to include the impact of the data-informed regional variance [START_REF] Amunts | Julich-Brain: A 3D probabilistic atlas of the human brain's cytoarchitecture[END_REF] which is now reachable by TVB through Ebrains [START_REF] Schirner | Brain simulation as a cloud service: The Virtual Brain on EBRAINS[END_REF]. Similarly intriguing direction for the extension of the framework presented here is in more refined inclusion of the subcortical structures, especially their impact through the neuromodulation. Notably, recent works [START_REF] Shine | The thalamus integrates the macrosystems of the brain to facilitate complex, adaptive brain network dynamics[END_REF][START_REF] Müller | Diffuse neural coupling mediates complex network dynamics through the formation of quasi-critical brain states[END_REF][START_REF] Munn | The ascending arousal system shapes neural dynamics to mediate awareness of cognitive states[END_REF] exploring the role of thalamus, locus coeruleus, and basal nucleus of Meynert in shaping of the dynamical landscape of the cortical activity are already formulated in the dynamical systems language while incorporating carefully the detailed anatomical and cytoarchitectural knowledge. Integrating these advances in the SFM framework is a natural next step towards the original motivation of SFM, that is to link the mesoscopic neuronal activity to the behaviour, as the intricate interactions between cortex and the subcortical areas are one of the organizing principles of the underlying the biological mechanisms supporting behaviour [START_REF] Cisek | Resynthesizing behavior through phylogenetic refinement[END_REF].

Parcellation-induced variation of empirical and simulated brain connectomes at group and subject levels is another issue that needs to be considered [START_REF] Domhof | Parcellationinduced variation of empirical and simulated brain connectomes at group and subject levels[END_REF].

Nevertheless, we focus on general mechanisms without going on regional level specificities, so the choice of parcelation should not play such a role.

In conclusion, our results show how the low-dimensional dynamics arises from breaking the symmetry in the brain on the level of the connectome. Describing these dynamics as structured flows on manifolds allows us to bridge the gap between the observational measures and the state-space trajectories of the system. As such, this object is well suited for comparison across different models, scales, and neuroimaging modalities, and provides means for integration of the diverse descriptions of the resting state dynamics.

4 Materials and methods

Brain network model

Computational brain network model [START_REF] Sanz-Leon | Mathematical framework for large-scale brain network modeling in The Virtual Brain[END_REF] is used to simulate resting state activity under varying values of network coupling scaling parameter G. The dynamics of each of the network nodes were governed by the neural mass model (NMM) derived analytically as the limit of infinitely all-to-all coupled θ-neurons [START_REF] Montbrió | Macroscopic description for networks of spiking neurons[END_REF],

namely for i-th node for the firing rate r i and membrane potential v i as:

τ c ṙi = ∆ πτ c + 2r i v i , τ c vi = v 2 i + η -(τ c πr i ) 2 + Jτ c r i + I i , (3)
where I i is the input current, η is the average neuronal excitability, J is the synaptic weight, ∆ is the spread of heterogeneous noise distribution, and τ c is the characteristic time.

The N nodes are then coupled with a connectome derived from empirical data as

I i (t) = G j W ij r j (t -D ij ), ( 4
)
where G is the network scaling parameter, W ij is the connection weight, D ij = L ij /S is the delay caused by propagation of the signal on a tract of length L ij with finite speed S. We picked the speed S = 2m/s from the biologically plausible range [START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF], and a connectivity matrix of a subject from the Human Connectome Project [START_REF] Van Essen | The WU-Minn human connectome project: an overview[END_REF] in the Desikan-Killiany parcellation [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] with 84 regions including subcortical structures.

The equations 3 and 4 comprise the drift a(Ψ, t) in the stochastic delay differential equation formulation with linear additive noise reading:

dΨ(t) = a(Ψ(t))dt + b(Ψ(t))dW (t), ( 5
)
where Ψ is the state vector [ψ 1 , ...ψ n ] with ψ n = [r n , V n ], dW (t) is a differential of a Wiener process with Gaussian increment with variance σ 2 , and b(Ψ, t) = 1

is the diffusion coefficient-here constant yielding the noise term additive.

The model was implemented in The Virtual Brain [START_REF] Sanz-Leon | The Virtual Brain: a simulator of primate brain network dynamics[END_REF] and equipped with BOLD forward solution comprising the Balloon-Windkessel model applied to the firing rate r(t) [START_REF] Stephan | Comparing hemodynamic models with DCM[END_REF].

The model parameters η = -5.0, J = 15.0, τ c = 1.0, and ∆ = 1.0 were selected to set the nodes in the bi-stable regime in the absence of coupling [START_REF] Montbrió | Macroscopic description for networks of spiking neurons[END_REF].

We then varied the global coupling parameter G and the noise variance σ, and simulated 10 minutes of resting state BOLD activity with T R = 720ms after discarding 10 seconds of the initial transient from random initial conditions.

Functional connectivity dynamics

In order to track the time-dependent changes in the functional connectivity, we compute the windowed dynamical functional connectivity dF C w [START_REF] Hansen | Functional connectivity dynamics: modeling the switching behavior of the resting state[END_REF][START_REF] Allen | Tracking whole-brain connectivity dynamics in the resting state[END_REF]

(i, j) = corr(F C(w i ) , F C(w j ) ), (6)
where F C(w) is the vectorized upper part of the F C matrix.

For the window-less dF C e [START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF] we start from the edge time-series [START_REF] Esfahlani | High-amplitude cofluctuations in cortical activity drive functional connectivity[END_REF] defined as E nm (t) = z n (t)z m (t) for n, m = 1 . . . N where z n (t) = Bn-µn σn is the z-scored BOLD time-series of a node n. The edge dynamical functional connectivity is then computed as correlation between the edge vectors at each pair of time points t 1 , t 2 :

dF C e (t 1 , t 2 ) = corr(E nm (t 1 ), E nm (t 2 )).

The co-fluctuation events (CF) are defined as time points in the edge timeseries E nm (t) during which the root sum squared RSS = nm E 2 nm (t) crosses a given threshold, here chosen as 95th percentile. Time points where RSS is below the threshold are then labeled as non-events (nCF).

The avalanches were computed on the binary mask a(t) on the r(t) such that a i (t) = 1 ⇐⇒ z(r i (t)) > 3 where z(r i (t) is the z-score of firing rate r of a node i.

Manifold subspaces

As a first step in the analysis of the local dynamics specific to a particular attractive subspace, we have extracted the time-points belonging to these subspaces with k-means clustering applied to the edge time series E nm (t). We varied the number of clusters k and selected k = 5 at which the co-fluctuation events separated to distinct cluster.

To extract the segments of r(t) corresponding to the E nm (t) time points we first estimated the BOLD signal lag l = 2500 ms as optimal peak-to-peak alignment with r(t) smoothened by a Gaussian filter with same effective width (σ = 700). Then for all BOLD time points in a given cluster c we selected the 2000 corresponding time points in r(t) and concatenate these to get the fast time-scale activity r c (t) in the subspace corresponding to cluster c. Each of the r c (t) was then projected to space spanned by the first two PCA components of the whole r(t) time-series to evaluate how much of the overall state-space is covered by individual clusters.

The local trajectory for a given event e was computed by selecting interval r e (t) corresponding to BOLD timepoints above the RSS threshold and three timepoints before and after the event. Local P CA e of was then computed from r e (t), and the smoothened trajectory was computed by convolving r e (t) with a Gaussian filter (σ = 100).

Manifold sampling

To identify the fixed point scaffold of the manifold as traced by the trajectory resulting from integrating the Equation 5, we sample the segments from the simulated trajectory (r i (t), v i (t))| ts+τmax ts , and use them as initial conditions for integration of the deterministic interpretation of Equation 5, i.e. dΨ(t) = a(Ψ(t))dt. From each such an initial condition, we integrated the system to steady state equilibrium corresponding to a fixed point (r * , v * ).

The number of stable fixed points (r * , v * ) of system with G = 0 is 2 N reflecting all the combinations of up-and down-states of the N nodes. To sample the stable fixed points of the system with G > 0 we solve repeatedly the system of equations:

0 = ∆ πτ c + 2r * i v * i , 0 = v * 2 i + η -(τ c πr * i ) 2 + Jτ c r * i + I i (8)
using Newton-Raphson method with the initial conditions chosen randomly as a vector of up-and down-state fixed points of the isolated nodes, i.e. (r

* 0 i , v * 0 i ) ∈ {(r * ↑ , v * ↑ ), (r * ↓ , v * ↓ )}, ∀i where (r * ↑ , v * ↑
) and (r * ↓ , v * ↓ ) are the up-and down-state fixed points for the isolated node. For each initial condition (r * 0 , v * 0 ) we then check if the corresponding solution of Equation 8 is equivalent up to the composition in terms of up-and down-states. If not, it is discarded, otherwise we evaluate the stability of the found fixed point using linear stability analysis.

As a low-dimensional projection of the sampled manifold we have used the two slowest eigenmodes of the structural connectivity. These are computed as eigendecomposition of the graph Laplacian L = W-I, that is LU = UΛ, where eigenvalues λ k can be interpreted as structural freqencies and the eigenmodes u k as structural connectome harmonics [START_REF] Preti | Decoupling of brain function from structure reveals regional behavioral specialization in humans[END_REF].

Linear stability analysis

We perform a linear stability analysis to identify the fixed-points obtained from the NR method. If each fixed-point (r * , v * ) is perturbed by ( r , v ) , then the evolution of the perturbations depend on the Jacobian matrix (J) and are given by:

              ˙ r 1 ˙ r 2 . . . ˙ r N ˙ v 1 ˙ v 2 . . . ˙ v N               =              
              .               r 1 r 2 . . . r N v 1 v 2 . . . v N               (9)
The stability of a fixed-point depends on the eigenvalues of the Jacobian evaluated at the fixed-point. The fixed-point is stable iff all the eigen-values of J are negative. Therefore, we numerically evaluate the largest eigenvalue of Jacobian for each fixed-point and label the point as stable if its real-part is negative.

Fixed point sampling from simulated trajectory

From a given trajectory of the system given as 10 minutes of ψ(t) we have selected a restart point t each 50 ms (12000 starting points altogether). For each of the restart point t we extracted the segment Ψ(t)| t -τmax t where τ max is the length of the longest delay, and used as initial condition to a equivalent system to Equation 5 with b = 0: dΨ(t) = a(Ψ(t))dt.

(10)

Integrating this system to equilibrium yielded then for each restart point t r a fixed point Ψ * = (r * , v * ). The stability of each of the fixed points Ψ * was then evaluated using the linear stability analysis as the largest eigenvalue of the respective Jacobian matrix.

Escape time analysis

The switching behaviour of a single node is driven by the stability of the upand down-state fixed points in the presence of noise. We employed escape time analysis [START_REF] Yamapi | Global stability analysis of birhythmicity in a self-sustained oscillator[END_REF] to measure the stability of these fixed points for range of values of external input I. In detail, for a single node of the system given by Equation 3we found the up-and down-state stable fixed points (r * , v * ) ↑ and (r * , v * ) ↓ , and the unstable saddle node (r * , v * ) × . Next we computed the separatrix between the two basins of attraction by integration of the model backwards in time resulting in an closed curve ω. To find the characteristic escape time for a fixed point (r * , v * ) we have integrated the system from the initial condition (r 0 , v 0 ) = (r * , v * ) for a given value of I 100 times, measuring the time t E at which the trajectory crosses ω for the first time. The values of I were drawn from the range given by [0, I max ] where I max = max{I i (t), ∀i} is the largest value of I i encountered in the integration of the full system in the working point.

Empirical data and spatial analysis

The functional gradient on empirical data was computed from the group connectivity matrix of the Human Connectome Project dataset using the brainspace toolbox [START_REF] Vos De Wael | BrainSpace: a toolbox for the analysis of macroscale gradients in neuroimaging and connectomics datasets[END_REF]. For a simulated resting state session with G w , the time in avalanche was computed for each node as total time for which the r i (t) was above the threshold of 3 standard deviations, and the event z-score as a sum of z-scored BOLD signal in time-points marked as events. The nodes were then groupped according to the cortical hieararchy [START_REF] Mesulam | From sensation to cognition[END_REF] projected to the Desikan-Killiany parcellation.

A parcellation-based BOLD signals of a resting-state session from a subject from the Human Connectome Project [START_REF] Domhof | Parcellationbased resting-state blood-oxygen-level-dependent (BOLD) signals of a healthy cohort[END_REF] were used to validate the separation of the events in the low-dimensional embedding. The data consisted of 1200 time points sampled at 720 ms in the Desikan-Killiany parcellation [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] with 70 cortical regions.

Figure 1 :

 1 Figure 1:Structured flows on manifolds as focus of resting state characterization. With respect to the structure of the connectivity of the dynamical system, we consider spectrum defined by the two symmetrical limit cases: fully connected and fully disconnected network. Driven by noise, the disconnected system exhibits fully statistical, high-dimensional dynamics -it explores the whole state space in a equidirectional manner. On the other hand, the dynamics of the fully connected system is fully constrained corresponding to a SO(n) hypersphere with zero flow. The dynamics on the sparsely connected system leads to an object in between -a low-dimensional attractive manifold with an associated flow (SFM). It is this object we wish to put in the center of interest and characterize. While the SFM object remains the same, connections are made to data of various modalities with the help of suitable data features.

ForFigure 2 :

 2 Figure 2: Brain network model and symmetry breaking. The brain network model is simulated for varying levels of global coupling parameter G and noise variance σ to produce both time-series of the state space variables r(t), V(t), and the BOLD signal. For each combination of G and σ we compute the sliding window dF C w matrix from the simulated BOLD signal, and quantify the "switching index" of the dFC as the variance of the upper triangle (A). Kolmogorov-Smirnov distance of the centered (mean-subtracted) distributions of the values of the upper triangle of the dFC computed from empirical and simulated resting state BOLD time series. The region of parameter space where the distributions are closest overlaps with the region with high fluidity of dFC (B). For selected values of (G, σ) we show the sliding window dF C w (C), edge based dF C e (D) and the projection of r(t) time series in the first two PCA components (E) annotated with corresponding fractional variance accounted for (VAF). In the working point around G = 0.54 and intermediate values of σ the system exhibits recurrence in the large-scale dynamics as captured by nonzero switching index, and reduction of dimensionality as captured in the increase in explained variance by the first PCA components and the asymmetry in the respective projection. For values of G below or above the working point, the systems loses the fluidity property as reflected in the absence of the off-diagonal blocks on the dFC, and exhibits high-dimensional dynamics.

Figure 3 :

 3 Figure 3: Network dynamics. (A) The network input I modulates the probability of a noise-driven transition between the down-and up-state. (B) Example of a cascade-coordinated increase in activity translating to a delayed correlated peak in the BOLD signal. Below we compare the network dynamics in and outside the working point, and the empirical data. In both empirical data and the working point (G = 0.54), the BOLD co-activations follow the neuronal cascades with a latency (C), and show distinct spatial profiles which are recurrent in time (D): edge time series on the (top panel) captures the spatial profiles of the co-activations, the similarity across time is captured by the dF C e matrix (middle panel), and the distributions of correlation between co-activation events (CA) and non-events (nCA) is compared (bottom panel).

Figure 4 :

 4 Figure 4: Manifold subspaces and characteristic dynamics. (A) The edge-based dynamic functional connectivity dF C e of a simulation of the model in the working point (top) shows the off-diagonal structure of similarity of the system's activity across time. The edge time series (middle) shows the time evolution of the functional connectivity of the simulated BOLD signal between each of the node pairs, and exhibits the characteristic co-activation events defined as time points with the root sum squared (RSS, bottom) crossing the thershold of 95th percentile. Dividing the edge-time series into 5 clusters (k-means, shown in the colorbar under the dFC) has separated the event and non-event time points, and also differentiated the events based on their respective similarity. (B) The time intervals in r(t) corresponding to the 5 clusters were selected; in the first panel the centroids of the time points of the individual clusters are marked with a cross in the projection to the first two principal components of the whole time series, following panels show the projection of the r(t) intervals of particular clusters. Cluster #2 captures the high-activity subspace, cluster #4 corresponds to the low-activity state, and the clusters #1, #3, and #5 capture the co-activation events. (C) Local trajectories in the manifold subspaces: the time series of the three example events (a,b,c, marked in the panel A bottom)was projected to the first two components of PCA applied to each time segment individually. The smoothened trajectory marks the advance of the system through the event and out of it, and is colored by the value of RSS (yellow at the peak of the event blue at the beginning and the end). (D) The event trajectories on the manifold. The trajectory of the simulated BOLD signal is projected in the space defined by the first three PCA components with the events colored by the RSS value (yellow at the peak of the event).

Figure 5 :

 5 Figure 5: Mechanistic structure of the manifold. (A) Composition of the sampled stable fixed points in terms of number of nodes in the up-state as a function of G, normalized to G = 0. (B) Projection of the stable fixed points into the first two leading eigenmodes of the network Laplacian E 1 , E 2 , color coded with the value of the largest eigenvalue in the linear stability analysis. (C) Fixed points (colored) derived by noise-free integration to equilibrium from the trace of a simulation (black) in the working point G w , color coded by the value of the largest eigenvalue λ 1 . (D) Frequency peak in the simulated source activity of each of the regions plotted against the node structural connectivity in-strength. (E) Empirical BOLD time series projected into the first two PCA components with the events colored by the RSS value (yellow in the peak). (F) Across the cortical hierarchy, the time spent in avalanches of the r(t) timeseries (top) decreases, as does the cumulative z-scored simulated BOLD from the event time-segments (middle). The spatial distribution of the principal functional gradient extracted from empirical fMRI is also aligned along the cortical hierarchy (bottom).

  analytical results of increased entropy of the attractors in an Ising-spin network model for intermediate values of coupling strength[START_REF] Deco | How anatomy shapes dynamics: a semianalytical study of the brain at rest by a simple spin model[END_REF]. The relationship to the dimensionality of the exhibited dynamics is such that for the low values of coupling strength G, where the Ising model is in the trivial state with all spins equal to 0, the model presented here is also driven by noise to the all-down state due to the significantly larger basin of attraction of the down-state, and the nodes make uncoordinated noise-driven excursions to the up-state reflected in the high-dimensionality of the dynamics. For high values of G the situation is opposite, and for intermediate values of G the Ising model exhibits high entropy of attractors, which is in our case reflected in the available states organized in the low-dimensional manifold with the structured flow governed by the stability of these states.

  and edge dynamical functional connectivity dF C e[START_REF] Rabuffo | Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest[END_REF]. Starting from the regional BOLD time-series B n (t) for each node n, we compute functional connectivity matrices F C(w) for each time window w = 1...W defined as B n (t)| tw+τ

tw with window length τ = 60s and window step size t (w+1) -t w = 2s. Next we compute the dF C w matrix of order W as dF C w