Mallory Carlu

Omar Chehab

Dalla Leonardo

Damien Porta

Charlotte Depannemaecker

Maciej Héricé

Elif Jedynak

Paolo Köksal Ersöz

Selma Muratore

Cristiano Souihel

Capone

L Dalla Porta

D Depannemaecker

C Héricé

M Jedynak

E Köksal Ersöz

P Muratore

S Souihel

C Capone

Y Zerlaut

A Destexhe

& M Di Volo

A mean-field

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

Introduction

Brain dynamics can be investigated at different scales, from the microscopic cellular scale, describing the voltage dynamics of neurons and synapses [START_REF] Markram | Reconstruction and simulation of neocortical microcircuitry[END_REF] to the mesoscopic scale, characterizing the dynamics of whole populations of neurons [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF], up to the scale of the whole brain where several populations connect together [START_REF] Bassett | On the nature and use of models in network neuroscience[END_REF][START_REF] Deco | Rethinking segregation and integration: contributions of whole-brain modelling[END_REF][START_REF] Sanz Leon | The virtual brain: a simulator of primate brain network dynamics[END_REF].

In their pioneering work [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF], Wilson and Cowan describe the dynamics of a population of neurons through a well-known differential equation where the input-output gain function is described by a sigmoid function. This approach inspired a long lasting research in neuroscience where population models, usually called rate models, permit a qualitative insight into the dynamics of population of neurons [START_REF] Di Santo | Landau-ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization[END_REF][START_REF] Hopfield | Neurons with graded response have collective computational properties like those of two-state neurons[END_REF][START_REF] Sompolinsky | Chaos in random neural networks[END_REF][START_REF] Sussillo | Generating coherent patterns of activity from chaotic neural networks[END_REF].

Moreover, a large effort has been made in order to derive population descriptions from the specificity of the network model under consideration. This bottom-up approach permits to obtain a dimensionally reduced mean f ield description of the network population dynamics in different regimes [START_REF] Amit | Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex[END_REF][START_REF] Brunel | Fast global oscillations in networks of integrate-and-fire neurons with low firing rates[END_REF][START_REF] Capone | Sleep-like slow oscillations improve visual classification through synaptic homeostasis and memory association in a thalamocortical model[END_REF][START_REF] Di Volo | Heterogeneous mean field for neural networks with short-term plasticity[END_REF][START_REF] El Boustani | A master equation formalism for macroscopic modeling of asynchronous irregular activity states[END_REF][START_REF] Montbrió | Macroscopic description for networks of spiking neurons[END_REF][START_REF] Ohira | Master-equation approach to stochastic neurodynamics[END_REF][START_REF] Renart | Mean-field theory of irregularly spiking neuronal populations and working memory in recurrent cortical networks[END_REF][START_REF] Schwalger | Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size[END_REF][START_REF] Tort-Colet | Attractor competition enriches cortical dynamics during awakening from anesthesia[END_REF][START_REF] Tsodyks | Rapid state switching in balanced cortical network models[END_REF][START_REF] Van Vreeswijk | Chaos in neuronal networks with balanced excitatory and inhibitory activity[END_REF][START_REF] Vreeswijk | Chaotic balanced state in a model of cortical circuits[END_REF]. On one hand mean f ield models permit a simpler, reduced picture of the dynamics of a population of neurons, thus allowing to unveil mechanisms determining specific observed phenomena (di [START_REF] Di Volo | Transition from asynchronous to oscillatory dynamics in balanced spiking networks with instantaneous synapses[END_REF][START_REF] Jercog | Up-down cortical dynamics reflect state transitions in a bistable network[END_REF][START_REF] Reig | Synaptic transmission and plasticity in an active cortical network[END_REF]. On the other hand, they enable a direct comparison with imaging studies where the spatial resolution implies that the recorded field represents the average over a large population of neurons (i.e. a mean-field) [START_REF] Capone | Slow waves in cortical slices: How spontaneous activity is shaped by laminar structure[END_REF][START_REF] Chemla | Suppressive traveling waves shape representations of illusory motion in primary visual cortex of awake primate[END_REF].

During awake states, cortical dynamics generally show asynchronous spiking activity, where individual neurons are characterized by an irregular (typically Poissonian) firing pattern [START_REF] Burns | The spontaneous activity of neurones in the cats cerebral cortex[END_REF][START_REF] Dehghani | Dynamic balance of excitation and inhibition in human and monkey neocortex[END_REF][START_REF] Softky | The highly irregular firing of cortical cells is inconsistent with temporal integration of random epsps[END_REF]. In this dynamical regime, so-called Asynchronous Irregular (AI), the correlation of the network activity decays relatively quickly in time, making it possible to develop a Markovian formalism in order to obtain mean-field equations.

The application of such a theory to binary neurons led to the derivation of dynamical equations for population rates [START_REF] Ginzburg | Theory of correlations in stochastic neural networks[END_REF][START_REF] Ohira | Master-equation approach to stochastic neurodynamics[END_REF]. More recently, such a theory has been extended to spiking neurons, permitting to obtain differential equations for neurons average activity and for higher-order moments [START_REF] Buice | Systematic fluctuation expansion for neural network activity equations[END_REF][START_REF] Dahmen | Correlated fluctuations in strongly coupled binary networks beyond equilibrium[END_REF][START_REF] El Boustani | A master equation formalism for macroscopic modeling of asynchronous irregular activity states[END_REF]. In the first order, these equations are formally the same as the rate models, like the Wilson-Cowan approach, although the function linking input-output properties of populations of neurons, namely the transfer function, is more complex than a sigmoid. Indeed, in this formalism, it encompasses internal properties of the neuronal models, together with the type of synaptic interactions under consideration, to yield a population scale description. In general, such function cannot be expressed in a closed form for complex neurons, especially if some realistic ingredients like conductance based interactions are taken into account.

2

In this article we present a general approach to determine the transfer function for complex models, from the Adaptive Exponential Integrate-and-Fire (AdEx) to the Hodgkin-Huxley (HH) and the Morris-Lecar (ML) models. As a result, we obtain mean-field equations for the population dynamics in AI regimes as observed in cortical regions for highly detailed models, creating a bridge between electrophysiology at the microscopic scale and the details of the famous transfer function first used by Wilson and Cowan as a sigmoid.

Finally, we test not only the ability of our mean-field models to describe spontaneous activity of the considered neuronal populations, but also their predictive power for network response to external stimuli. We show that, provided the stimuli are fairly slow, the mean field model gives good quantitative predictions.

Materials and Methods

We describe here the neuronal and network models used in this study. We also introduce mean-field equations describing population dynamics and the template to estimate the transfer function that we apply to all the neuronal models under consideration.

Network of spiking neurons

We consider a random directed network of N = 10 4 cells, among which 80% are regular spiking (RS) excitatory (E) and 20% are fast spiking (FS) inhibitory (I) neurons. The connections between pairs of neurons are set randomly with a fixed probability (p = 0.05). Unless otherwise stated, the same network and synaptic constants are used for all the neuronal models (Hodgkin-Huxley, Adaptive Exponential Integrate-and-Fire and Morris-Lecar). The dynamics of each node k follows:

ẋk = F (x k ) + I syn , (1)
where x and F (x) represent the neuronal state and dynamics, the latter depending on the specific model (see Section 2.2). Note the notation xk which indicates that, in general, each neuron is characterized by a vector of variables. The synaptic current impinging on the postsynaptic neurons k, I syn , is modeled as:

I syn = (E e -v k )G e syn + (E i -v k )G i syn , (2)
where E e = 0 mV (E i = -70 mV) is the excitatory (inhibitory) reversal potential and G (e,i) syn is the conductance modeled as a decaying exponential function:

G (e,i) syn (t) = Q (e,i) n Θ(t -t sp (n))e t-tsp (n) τ , (3)
where

Q e (Q i ) is

Single neuron models

We describe here the neuronal models used in the rest of the paper, starting from the Integrateand-Fire up to the Morris-Lecar and Hodgkin-Huxley models.

Adaptive Exponential Integrate-and-Fire model

The dynamics of each of the AdEx neurons i is described by the following 2D (here xi = (v i , w i ))

differential equations [START_REF] Brette | Adaptive exponential integrate-and-fire model as an effective description of neuronal activity[END_REF] :

c m dv i dt = g L (E L -v i ) + g L ∆e v i -v t ∆ -w i + I syn , (4)
dw i dt = - w i τ w + b tsp(i) δ(t -t sp (k)) + a(v i -E L ), (5)
where c m = 150pF is the membrane capacity, v i is the voltage of neuron i and, whenever v i > v t = -50 mV at time t sp (i), v i is reset to the resting voltage v rest = -65 mV and fixed to this value for a refractory time T refr = 5 ms. The leak term has a fixed conductance of g L = 10 nS and the leakage reversal E L = -65 mV, if not stated otherwise. The exponential term has a different strength for regular-spiking (RS) and fast-spiking (FS) cells, i.e. ∆ = 2 mV (∆ = 0.5 mV) for excitatory (inhibitory) cells. The variable w mimicks the dynamics of spike frequency adaptation.

Inhibitory neurons are modeled according to physiological insights as the FS neurons with no adaptation while the excitatory RS neurons have a lower level of excitability due to the presence of adaptation. Here we consider b = 60 pA, a = 4 nS and τ w = 500 ms, if not stated otherwise.

Hodgkin-Huxley

The dynamics of the Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] is given by the following 5-dimensional system of differential equations [START_REF] Pospischil | Minimal hodgkin-huxley type models for different classes of cortical and thalamic neurons[END_REF]:

c m dv i dt = g L (E L -v i ) + g N a m 3 i h i (E N a -v i )+ + g K n 4 i (E K -v i ) + g M p i (E K -v i ) + I syn , (6)
dn i dt = α n (v i )(1 -n i ) -β n (v i )n i , (7)
dm i dt = α m (v i )(1 -m i ) -β m (v i )m i , (8)
dh i dt = α h (v i )(1 -h i ) -β h (v i )h i , (9)
dp i dt = (p ∞ (v i ) -p i )/τ p (v i ), (10)
with the gating functions,

α n (v) = -0.032(v -V T -15) exp[-(v -V T -15)/5] -1 , β n (v) = 0.5 exp -(v -V T -10) 40 , α m (v) = -0.32(v -V T -13) exp[-(v -V T -13)/4] -1 , β m (v) = 0.28(v -V T -40) exp[(v -V T -40)/5] -1 , α h (v) = 0.128 exp[-(v -V T -17)/18], β h (v) = 4 1 + exp[-(v -V T -40)/5] , p ∞ (v) = 1 1 + exp[-(v + 35)/10] , τ p (v) = τ max 3.3 exp[(v + 35)/20] + exp[-(v + 35)/20] , (11)
where v i is the voltage and (n i , m i , h i , p i ) are the corresponding gating variables of the i-th neuron.

We set the spike emission times t sp (k) for this model to time steps in which the membrane potential v exceeded a voltage threshold of 10 mV. Unless stated otherwise, the membrane capacitance c m = 200 pF/cm 2 , the maximal conductance of the leak current g L = 10 mS/cm 2 , the sodium current g N a = 20 mS/cm 2 , the delayed-rectifier potassium current g K = 6 mS/cm 2 , the slow non-inactivating potassium current of the excitatory (RS) neurons g M = 0.03 mS/cm 2 and of the inhibitory (FS) neurons g M = 0 mS/cm 2 , with corresponding reversal potentials E L = -65 mV, E N a = 50 mV, E K = -90 mV, the spiking threshold V T = -53.5 mV and τ max = 0.4 s are the fixed parameter values in Eqs. ( 6)-(11).

Morris-Lecar

The dynamics of the Morris-Lecar model [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF] is described by the system of differential equations:

c m dv i dt = g L (E L -v i ) + g Ca M ss (v i )(E Ca -v i ) + g K N i (E K -v i ) + I syn + I 0 , ( 12
)
dN i dt = N ss (v i ) -N i τ N (v i ) , (13)
where c m = 2µ F/cm 2 is the membrane capacitance, v i is the membrane potential in mV, N i and M ss are the fraction of open potassium and calcium channels, respectively. The current I 0 = 0.2 nA/cm 2 is a reference DC external current. Spike emission times are established in the same way as for the HH model. The maximal conductances for the leakage current (L), calcium (Ca) and potassium (K) were fixed to g L = 20 mS/cm 2 , g Ca = 80 mS/cm 2 and g K = 160 mS/cm 2 , respectively. The reversal potentials are E L = -50 mV for excitatory RS neurons and E L = -70 mV for inhibitory FS neurons, E Ca = 120 mV and E K = -84 mV. The quantities M ss and N ss are modeled as:

M ss (v) = 1 2 1 + tanh v -V 1 V 2 , ( 14
)
N ss (v) = 1 2 1 + tanh v -V 3 V 4 , (15)
with

τ N (v) = 1 2 φ cosh v -V 3 2V 4 , (16)
where

V 1 = -1.2 mV, V 2 = 18 mV, V 3 = 2 mV, V 4 = 30
mV are tuning parameters that determine the half activating voltage and slope of the activation curves for calcium and potassium conductances. This choice of parameters is such that the ML neuron is set in a type II excitability class, i.e. its response to a DC current is discontinuous and the neuron firing rate increases very slowly with the injected current (data not shown).

Mean-field formalism

Mean-field theory scales the analysis of interacting point-wise neurons to their macroscopic, collective, dynamics based on the moment-statistics of the system, requiring a self-averaging hypothesis for physical quantities. We make here an additional hypothesis that the biological neural network 5 is set to asynchronous irregular dynamical regime. The latter is chosen for its biological plausibility [START_REF] Destexhe | The high-conductance state of neocortical neurons in vivo[END_REF] as observed in awake cortical states of adult mammalian brains.

We use here the master equation formalism reported by [START_REF] El Boustani | A master equation formalism for macroscopic modeling of asynchronous irregular activity states[END_REF] providing a system of ordinary differential equations that describe the evolution of the mean and variance of the firing rate of excitatory and inhibitory neurons. The central argument for this derivation is to consider the network dynamics as markovian on an infinitesimal (a time resolution T, typically 20ms) scale, as in [START_REF] Buice | Systematic fluctuation expansion for neural network activity equations[END_REF]; [START_REF] Ginzburg | Theory of correlations in stochastic neural networks[END_REF]; [START_REF] Ohira | Master-equation approach to stochastic neurodynamics[END_REF]. Moreover, such a theory is based on the assumption that neurons emit maximum one spike over the Markovian step T , meaning that the theory assumes relatively low firing rate of neurons, lower than 1/T ∼ 50 Hz [START_REF] El Boustani | A master equation formalism for macroscopic modeling of asynchronous irregular activity states[END_REF], as typically is the case in the asynchronous irregular regimes here investigated. The differential equations read:

T dν µ dt = (F µ -ν µ ) + 1 2 c λη ∂ 2 F µ ∂ν λ ∂ν η , (17)
T dc λη dt = δ λη F λ (1/T -F η ) N λ + (F λ -ν λ )(F η -ν η ) + ∂F λ ∂ν µ c ηµ + ∂F η ∂ν µ c λµ -2c λη , (18)
where µ = {e, i} is the population index (excitatory or inhibitory), ν µ the population firing rate and c λη the covariance between populations λ and η. The function

F µ={e,i} = F µ={e,i} (ν e , ν i ) is
the transfer function which describes the firing rate of population µ as a function of excitatory and inhibitory inputs (with rates ν e and ν i ). At the first order, i.e. neglecting the dynamics of the covariance terms c λη , this model reduces to the well known Wilson-Cowan model, with the specificity that the functions F need to be obtained according to the specific single neuron model under consideration. We introduce this procedure in the next section.

Transfer function estimate

The transfer function relates the firing rate of a neuron's response to its presynaptic excitatory and inhibitory firing rates. The particular form of the transfer function is related to the dynamics describing neuronal activity. Deriving an analytical expression for the transfer function is a nontrivial endeavor due to the nonlinear character of the dynamics, e.g. through conductance based interactions. Therefore, we use here a semi-analytic approach to fit a family of plausible transfer functions to the data obtained by means of numerical simulations with the desired neuron types.

This method, developed first by [START_REF] Zerlaut | Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime[END_REF] on data from experimental recordings, is based on the assumption that the transfer function depends only on the statistics of the subthreshold membrane voltage dynamics which is assumed to be normally distributed. These statistics are:

the average membrane voltage µ V , its standard deviation σ V and auto-correlation time τ V . Under these assumptions the neuronal output firing rate F ν is given by the following formula:

F ν = 1 2τ V erfc V eff thre -µ V √ 2σ V , ( 19
)
where erfc is the Gauss error function, V eff thre is an effective or phenomenological threshold accounting for nonlinearities in the single-neuron dynamics. Note that when dealing with extremely high spiking frequencies, e.g. in the case of Hodgkin-Huxley model close to depolarization block, a multiplicative factor α can be added in front of r.h.s. of Eq. ( 19) to permit the fitting procedure to deal with such high frequencies. In the asynchronous irregular dynamical regime, investigated 6 in this work, neurons have relatively low firing rates (smaller than 30 Hz). Accordingly, we never use this extension (i.e. the inclusion of the factor α) apart from the inset of Fig. 2b where we fit the transfer function of the Hodgkin-Huxley model over a broad range of frequencies, including those close to depolarization block where the firing rate is around 500-600 Hz. For this case we used α = 2. In the following section we introduce how the quantities µ V , σ V , τ V can be expressed as functions of the presynaptic excitatory and inhibitory firing rates ν E and ν I .

Table 1: Fit Parameters AdEx Neurons (in mV )

Cell Type We start by calculating the averages (µ Ge,Gi ) and standard deviations (σ Ge,Gi ) of the conductances given by Eq. 3 under the assumption that the input spike trains follow the Poissonian statistics (as is indeed the case in asynchronous irregular regimes here considered). In such case we obtain (Zerlaut and Destexhe, 2017a):

P 0 P µV P σ V P τ V P µ 2 V P σ 2 V P τ 2 V P µ V σ V P µ V τ V P σ V τ V RS-
P µV P σ V P τ V P µ 2 V P σ 2 V P τ 2 V P µ V σ V P µ V τ V P σ V τ V RS-Cell -48
P σ V P τ V P µ 2 V P σ 2 V P τ 2 V P µ V σ V P µ V τ V P σ V τ V RS-
µ Ge (ν e , ν i ) = ν e K e τ e Q e , σ Ge (ν e , ν i ) = ν e K e τ e 2 Q e , µ Gi (ν e , ν i ) = ν i K i τ i Q i , σ Gi (ν e , ν i ) = ν i K i τ i 2 Q i , (20)
where K i,e is the average input connectivity received from the excitatory or inhibitory population (in our cases typically K e = 400 and K i = 100) and in our model τ e = τ i = τ (see Eq. 3).

The mean conductances will control the total input of the neuron µ G and therefore its effective membrane time constant τ eff m :

µ G (ν e , ν i ) = µ Ge + µ Gi + g L , τ eff m (ν e , ν i ) = c m µ G . ( 21
)
Here we make the assumption that the subthreshold moments (µ V , σ V , τ V ) are not affected by the dynamics of the currents coming into play at the spiking time (e.g. sodium channels dynamics or the exponential term of the AdEx model). We thus consider, for all neurons, only the leakage term and the synaptic input in order to estimate subthreshold moments. Accordingly, we can write the equation for the mean subthreshold voltage:

µ V (ν e , ν I ) = µ Ge E e + µ Gi E i + g L E L µ G . ( 22
)
The final formulas for σ V and τ V follow from calculations introduced in Zerlaut et al. ( 2018) and they read:

σ V (ν e , ν i ) = s K s ν s (U s • τ s ) 2 2 (τ eff m + τ s ) , ( 23
) τ V (ν e , ν i ) = s K s ν s (U s • τ s ) 2 s K s ν s (U s • τ s ) 2 /(τ eff m + τ s ) , ( 24
)
where we defined U s = Qs µ G (E s -µ V ) and s = (e, i). Notice that neglecting all the spiking currents becomes a poorer assumption as the neuron activity increases. Nevertheless, we consider here AI dynamics where neurons have typically low firing rates (of the order few Hz). Moreover, as we show in the following sections, the fitting procedure will account for discrepancies in the actual evaluation of voltage moments by permitting an accurate prediction of neuron output firing rate.

From sub-threshold voltage moments to the output firing rate

The quantities µ V , σ V and τ V , obtained in the previous section, can now be plugged into Eq. 19

when an additional relation is taken into account. This relation follows from theoretical and experimental considerations [START_REF] Zerlaut | Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime[END_REF] showing that the voltage effective threshold

V eff thre can be expressed as a function of (µ V , σ V , τ V ). In [START_REF] Zerlaut | Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime[END_REF] the phenomenological threshold was taken as a second order polynomial in the following form:

V eff thre (µ V , σ V , τ N V ) = P 0 + x∈{µ V ,σ V ,τ N V } P x • x -x 0 δx 0 + x,y∈{µ V ,σ V ,τ N V } 2 P xy • x -x 0 δx 0 y -y 0 δy 0 , (25)
where we introduced the quantity τ N V = τ V G l /c m . We evaluated {P } through a fit according to simulations on single neurons activity setting first µ 0 V = -60 mV, σ 0 V = 0.004 mV, (τ N V ) 0 = 0.5, δµ 0 V = 0.001 mV, δσ 0 V = 0.006 mV and δ(τ N V ) 0 = 1. By the fitting procedure we find the values of the P parameters for the three neuronal models considered here (additionally for each model we consider two neuronal types: RS and FS) and we report the results in Table 3. In the first part of the Results section we describe the goodness of this procedure for the three considered neuronal models.

Results

We present here the results of a comparison between mean-field predictions and direct simulations.

We first test the technique to estimate the transfer function of single cells in AdEx, Hodgkin-Huxley and Morris-Lecar models and then compare theoretical predictions of the mean-field to numerical simulation of sufficiently large networks of neurons. 8

Transfer function for integrate-and-fire models

The transfer function of a simple AdEx neuron can be straightforwardly estimated by numerical simulations. As we report in Fig. 1 its shape is very similar to a sigmoidal function (as in the seminal paper by Wilson and Cowan) but its specific parameters follow from a complex combination of microscopic information, e.g. neurons resting potential.

See the black and blue dots in Fig. 1 for different values of the leakage reversal potential E L . Two main spiking modes can be distinguished in the neuronal dynamics. One, characterized by low output firing rate, where spikes are strongly driven by the membrane voltage fluctuations, namely fluctuation driven mode (see the bottom inset in Fig. 1). The second mode is characterized by a highly deterministic and regular firing observed at very high output firing rates (larger than 40 -50 Hz, top inset).

EL=-65mV EL=-80mV

Excitatory input (Hz)

Firing rate (Hz) The agreement remains very good for relatively low neuronal activity (up to 50 Hz). This is a direct consequence of the semi-analytic approach that assumes that neurons fire in an irregular manner (as observed in cortical dynamics) strongly driven by fluctuations around the mean membrane 9 voltage. In this work we only consider Asynchronous Irregular population dynamics for which the activity of neurons is low, irregular and strongly fluctuation driven.

Transfer function for complex models.

We report here the application of the techniques described in the Materials and Methods section to evaluating the transfer function of more complex neuronal models. To this end we consider We consider two kinds of neurons in agreement with neurophysiolgical information about cor-10 tical cells: excitatory neurons modeled as RS cells with a lower gain of the transfer function and inhibitory neurons modeled as FS cells with a higher gain. A different gain of the transfer function can be obtained by changing the excitability of the cells through their resting potential, or increasing the adaptation strength (see Materials and Methods for details).

By comparing the theoretical prediction with numerical simulation we observe that, for the three models considered here, the transfer function is correctly estimated both for inhibitory neurons (FS cells) and excitatory neurons (RS cells). This result shows that, even by considering a much more complicated model than AdEx it is possible to have access to a semi-analytic form of its transfer function and, importantly, to modify neurons excitability in the way allowing to obtain a similar transfer function (of excitatory RS and inhibitory FS cells) between different models.

Notice that the ML model shows a decrease of firing rate at frequencies higher than 8 Hz (i.e. no voltage oscillations and thus no firing activity), as reported previously for this model by [START_REF] Kim | The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons[END_REF]. This is a consequence of the depolarization block (DB) observed at high input frequencies (i.e. high average external current). Accordingly, we obtain a bell shaped transfer function, well predicted by our semi-analytical formalism. In previous studies this effect was taken into account in the context of Wilson-Cowan equations by using a Gaussian transfer function, instead of a sigmoidal [START_REF] Meijer | Modeling focal epileptic activity in the wilson-cowan model with depolarization block[END_REF], permitting to study the effect of the depolarization block in focal seizures at the population scale. In our model this shape, resembling a Gaussian, follows directly from Morris-Lecar equations, through the semi-analytical fitting. More specifically, in [START_REF] Meijer | Modeling focal epileptic activity in the wilson-cowan model with depolarization block[END_REF] the DB was studied in the Hodgkin-Huxley model. Indeed, also the HH model shows a DB but, at variance with the Morris-Lecar case, it appears in our parameters setup at very high firing rates, around 600-700 Hz (see the inset of Fig. 2b). In our simulations we do not consider this dynamical regime, being far from the dynamics typical of neurons in asynchronous irregular regimes. Moreover, as described in Fig. 1, the semi-analytic fitting procedure works well for low firing rates and discrepancies appear at high rates, also in the case of the simpler AdEx model. Nevertheless, we report here that, when performing the fitting over a wide range of input and output rates (see Methods), it is possible to obtain an overall good fit of the bell-shaped transfer function (see the inset of Fig. 2b).

Beyond the methodological point, our results show that even if the details of the mechanisms that generate a specific transfer function are very different, it is possible to adjust neurons parameters (e.g. excitability) in a way allowing to obtain similar transfer functions (at least in the region before entering a depolarization block). As a consequence, according to the mean-field theory,

where what matters for the population dynamics is only the transfer function, we expect different models to have a comparable emergent dynamics at the population (collective) scale.

Asynchronous irregular dynamics and mean-field predictions

In this section we compare the mean-field predictions of the emergent dynamics of networks of AdEx, HH and ML neurons, in the setup described in the previous section. In particular, we simulate a sparse network of RS and FS cells (see Fig. 2) coupled through conductance based interactions (see Materials and Methods for details). By looking at Fig. 2 we observe that, before reaching the DB, all three models have similar transfer functions, with FS neurons having a higher gain with respect to RS neurons, approximately of factor 3-4. As a result, we expect the population activity in the three models to fall within a similar dynamical regime, as a natural consequence of the mean-field's sole dependence on transfer functions, previously stated.

Indeed, by looking at Fig. 3 we observe that in the different networks the dynamics stabilizes on an asynchronous regime. In all cases, this regime is characterized by irregular microscopic dynamics (neuron's spiking statistics are Poissonian, data not shown) and represents the typical spiking patterns recorded during awake states in cortical regions (the autocorrelation function of population rate decreasing to zero in the time scale of tens of milliseconds). Moreover, as expected, 11 inhibitory FS cells fire at a higher frequency with respect to RS cells. Through the mean-field model it is possible to measure both the average population rate and its covariance (second order meanfield, see Materials and Methods). As reported in Fig. 3 we show that the mean-field model gives a good quantitative prediction of both quantities when they are compared to the histogram obtained by sampling the population rate in the network simulation. The higher discrepancy we observe for the complex neuronal models (e.g. HH and ML case) is related to a higher mismatch of the transfer function linked to the higher complexity of the model.

Network response to external stimuli

In order to complete the comparison between the mean-field model and the network dynamics, we study the response of the system to external stimuli. In particular, we consider an incoming

Poissonian train of spikes characterized by time-varying frequency and targeting both excitatory and inhibitory cells according to the following equation:

12

ν(t) = A Θ(t 0 -t)e -(t-t 0 ) 2 T 2 1 + Θ(t -t 0 )e -(t-t 0 ) 2 T 2 2 , ( 26
)
where Θ is the Heaviside function, and T 1 and T 2 are the rise and decay time constants, respectively. In Fig. 4 we report the comparison between the mean-field prediction and the network dynamics. By looking first at the AdEx and HH models, we observe that both mean-field models under investigation compare favorably with their corresponding network dynamics. We also verified, as it has been shown in (di Volo et al., 2019), that the faster the input dynamics is, the worse the agreement becomes. Indeed, for the Markovian hypothesis to hold, we need the time scale T to be much larger than the autocorrelation time in the spontaneous activity T ∼ τ m ∼ 10 ms.

Considering now the case of the ML model, we observe that by looking at Fig. 2 a relatively 13 strong input would bring single neurons to a depolarization block, which appears at relatively low activity levels. According to this difference with respect to AdEx and HH models, we would expect the population dynamics to show different properties in response to external perturbations. Indeed, as reported in Fig. 4c, the response to an external stimulus is very different from the one observed in the HH and AdEx models. In fact, in this case the excitatory stimuli turns out to inhibit both population activities. This anti-correlation between population input and output is well captured in its time course also by the mean-field model. This result shows that also for a more complex and highly non-linear setup the mean-field model is capable of predicting the ongoing activity and the time course of the response of a network of neurons operating in the asynchronous irregular dynamical regime.

Finally, we compare the results of the first and second order mean field on average population rates. In Fig. 4 we superimpose the continuous green (red) line for excitatory (inhibitory) rate obtained with the second order mean file with the results obtained with the first order (black dots).

We observe that the two quantities almost overlap (the difference is too small to be appreciated at this scale). Nevertheless, it is worth noticing that the second order mean field permits to obtain non-trivial information on the population dynamics and its fluctuations in time, with good quantitative predictions of the covariance of population rates (see the histograms in Fig. 3 and shadows in Fig. 4)).

Discussion

In this manuscript, we have reviewed a formalism to derive mean-field models from networks of spiking neurons and we have applied it to different complex neuronal models. The key to derive such "biologically realistic" mean-field models is to be able to obtain the transfer function of complex neuronal models. The approach we followed used a mean-field formalism based on a

Master Equation and which is applicable to every neuron, provided the transfer function is known [START_REF] El Boustani | A master equation formalism for macroscopic modeling of asynchronous irregular activity states[END_REF]. More recently, we have shown that the usual mathematical form of the transfer function, known analytically for the Integrate-and-Fire model, can capture more complex neuronal models [START_REF] Zerlaut | Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons[END_REF][START_REF] Zerlaut | Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime[END_REF]. This gave rise to a "semi-analytic" approach, where the transfer function is parameterized and fit numerically to the neuron model, while the mean-field remains analytic as only the parameters are obtained from the fitting. This approach was applied to the AdEx model (di [START_REF] Volo | Dynamical ventral tegmental area circuit mechanisms of alcohol-dependent dopamine release[END_REF][START_REF] Zerlaut | Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons[END_REF], and we extend it here to more complex models, namely the Morris-Lecar and the Hodgkin-Huxley models.

It is important to note that we limited here to "simple" firing patterns, i.e. neurons fire tonically in response to an external stimulus. In this setup the transfer function is well defined as the neuron's firing rate defines completely the spiking pattern. In cases where neurons exhibit different kind of activity, e.g. bursting, a different approach needs to be employed (see [START_REF] Ostojic | From spiking neuron models to linear-nonlinear models[END_REF]). Nevertheless, in the context of tonic neuronal activity the method is shown to be able to capture the response function of highly realistic models. We have studied here the predictions of the considered mean-field models on networks dynamics of excitatory RS and inhibitory FS cell populations during asynchronous irregular regimes, as observed in awake cortical activity. The results positively compare in the case of Morris-Lecar and Hodgkin-Huxley models for both the average and the variance of network population activity.

The good predictions at the population levels in the framework of the asynchronous irregular regimes is strongly dependent on the goodness of the fitting procedure for single neurons transfer functions. Even if such procedure works very well for neurons working in a low rate regime, whenever the firing rate becomes very high (higher than 100 Hz) the quantitative agreement gets worse. A more refined technique for the evaluation of the transfer function in different states (low and high rates activity) is an important topic for future research (recent work has addressed this 14 issue in AdEx model (Capone et al., 2019a)). A step forward in this direction can be important when dealing with neurons entering depolarization block at high firing rates, a mechanism playing an important role in focal seizures [START_REF] Meijer | Modeling focal epileptic activity in the wilson-cowan model with depolarization block[END_REF] or in dopaminergic neurons under normal condition or under drugs assumption (Di Volo et al., 2019;[START_REF] Dovzhenok | Exploring neuronal bistability at the depolarization block[END_REF]. Both in Morris Lecar and Hodgkin-Huxley model, the semi-analytic fitting is found to give quite good predictions on the presence of depolarization block, especially in the Morris Lecar case as this setup does not consider very high spiking frequencies. Even if work remains to be done in order to extend this framework to obtain a more reliable quantitative prediction on the depolarization block at very high frequencies, these preliminary results indicate the possibility to use these mean field techniques to connect the physiology at the cellular scale with pathological states at the population level, as the case of focal seizures.

We also reported that, in the framework of the Asynchronous regimes here considered, corrections to first order mean field due to second order terms (see Eq. ( 17)-( 18)) were relatively small but gave a good quantitative indication on the covariance of population rates (see histograms in Fig. 2). Nevertheless, in the case of dynamical regimes with higher neuronal correlation with respect to the ones here considered, we expect the second order mean field (taking explicitly into account the dynamics of covariances), to play an important role in the prediction of population average collective dynamics.

Let us point out that the goodness of the mean-field prediction depends also on the emergent dynamics of the network, i.e. in a highly synchronous dynamical regime the Markovian hypothesis fails and the mean-field model cannot give accurate predictions. Nevertheless, even if light synchronization is considered, e.g. during slow waves sleep, the mean-field models has been shown

to correctly predict such collective oscillations (di [START_REF] Volo | Dynamical ventral tegmental area circuit mechanisms of alcohol-dependent dopamine release[END_REF]. In this case it is however necessary to consider a mean-field model that includes the slow dynamics of spike frequency adaptation or that of the I M current in the case of Hodgkin-Huxley model. The possibility to include a conductance based adaptation to this formalism, e.g. by considering the slow dynamics of I M current, is a stimulating perspective for future works and will permit to obtain mean-field models for realistic neuronal models beyond the asynchronous irregular regime.

Moreover, beyond the input-output transfer function used here, a more complex transfer function has been used in order to take into account other features of neuron response, e.g. response in frequency [START_REF] Ostojic | From spiking neuron models to linear-nonlinear models[END_REF]. The addition of variables to account for a richer spiking pattern is an interesting direction, in case one is interested in modeling brain regions characterized by non-tonic firing of neurons (e.g. bursting cells in the thalamus). The general framework presented here could be extended in this direction, as it has been done to account for spike frequency adaptation yielding slow oscillations at the population scale.

Another possible extension is to apply the same formalism to complex neuronal models that include dendrites. A first attempt has been made in this direction (Zerlaut and Destexhe, 2017b) by considering simple "ball and stick" neuron models, where some analytic approximation is possible.

In principle, it should be possible to apply this approach to models based on morphologicallyreconstructed neurons, and calculate the transfer function of such models (work in progress).

This will lead to mean-field models based on morphologically realistic neuronal models. However, the presence of dendritic voltage-dependent currents complicates this approach, and should be integrated in the formalism. This constitutes an exciting future development of our approach.

Finally, the positive results obtained here for complex models, by showing the generality of our approach, motivate the future step of the application of this technique directly to experimental data. To this end, neurons must be recorded intracellularly in the absence of network activity (as typically in vitro), and many combinations of excitatory and inhibitory inputs must be injected as conductances (using the dynamic-clamp technique). The first attempt of this sort was realized on the Layer 5 neurons from mouse primary visual cortex [START_REF] Zerlaut | Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime[END_REF], where the transfer function could be reconstructed for a few dozen neurons. The same dynamic-clamp experiments should be done in the future to characterize the transfer function of inhibitory interneurons. Based on such experiments, it will be possible to obtain a mean-field model based on the properties of real neurons. Such a model will evidently be more realistic than the models we have presented here, which must be considered as a first step towards a quantitative population modeling of cerebral cortex and other brain regions.

  the excitatory (inhibitory) quantal conductance. The variable τ = 5 ms is the decay timescale of excitatory and inhibitory synapses and Θ is the Heaviside step function. The summation runs over the over all the pre-synaptic spiking times t sp (n). For both Hodking-Huxley and Adaptive Exponential Integrate-and-Fire models we set Q e = 1.5 nS and Q i = 5 nS, while for Morris-Lecar model Q e = 4 nS and Q i = 10 nS. On top of inputs coming from other neurons in the network, each excitatory and inhibitory neuron receive an external drive in the form of a Poissonian excitatory spike train at a constant firing rate ν drive = 4 Hz, if not stated otherwise.3

Figure 1 :

 1 Figure 1: Transfer function for an Exponential Integrate-and-Fire model Dots indicate the results of the numerical simulation of the Exponential Integrate-and-Fire model (FS cell, see Materials and Methods). The continuous line illustrates the results based on the semi-analytic fitting. The inhibitory Poissonian spike train used here has a fixed rate, r I = 8 Hz, while we show neuron average output Firing rate as the function of the Poissonian excitatory input spike train of rate r E . In the insets we show two exemplary voltage time traces corresponding to high (top inset) and low (bottom inset) firing rate. Colors stand for different values of the leakage reversal potential as indicated in the bottom-right corner of the figure.

  the well-known Hodgkin-Huxley (HH) model and the Morris-Lecar (ML) models (see Materials and Methods). These models permit to describe the details of Sodium and Potassium channels dynamics neglected in the simpler Integrate-and-Fire model and reproduce time evolution of the action potential. The semi-analytic approach to fit the numerical transfer function can be applied exactly in the same way as for AdEx models (as discussed in the Materials and Methods section).

Figure 2 :

 2 Figure 2: Transfer function for RS and FS cells: AdEx, HH and ML. We report the output firing rate for excitatory RS ( green) and inhibitory FS (red) cells obtained from numerical simulation (dots) and from the semi-analytic approach for the transfer function (continuous line). The inhibitory Poissonian spike train has a fixed rate, r I = 8 Hz. The bottom panels shows the time traces of the membrane voltage of an RS cell for an excitatory input equal to 4 Hz. Left column is obtained for the AdEx model, middle column for the HH model and right column for the ML model (see Materials and Methods). In the inset of panel (b) we report the transfer function for the RS cell estimated over very large values of input rates. In this case a separate fit by considering a broad input frequency range has been performed ( see Materials and Methods).

Figure 3 :

 3 Figure 3: Mean-field predictions and spontaneous activity: AdEx, HH and ML Top panels show raster plots for excitatory (green dots) and inhibitory (red dots) neurons, i.e. the spiking times for each neuron. Bottom panels shows histograms (obtained on a time length T w = 10 s) of population firing rated for excitatory (green) and inhibitory (red) populations. The Gaussian distribution has been plotted from mean-field predictions giving access to average firing rate and its variance. The left column (a) is obtained for the AdEx model, the middle column (b) for the HH model and the right column (c) for the ML model (see Materials and Methods).

Figure 4 :

 4 Figure 4: Population response to external stimuli: AdEx, HH and ML Top panels show the raster plot for excitatory (green dots) and inhibitory (red dots) neurons in response to an external excitatory stimulus (black dashed line in lower panels). Lower panels show the corresponding population rate (noisy line) together with the mean and standard deviation over time predicted by the the second order mean field model (red for inhibition and green for excitation). Superimposed the result obtained for the mean-field at the first order (black dots), which are almost coincident with results at the second order. Left column is obtained for the AdEx model, middle column for the HH model and right column for the ML model. Parameters are the same as in Fig. 3 and the external input (see Eq. (26)) has parameters A = 2 Hz, T 1 = 100 ms, T 2 = 150 ms for AdEx and HH and A = 2 Hz, T 1 = 100 ms, T 2 = 150 ms for ML, with t 0 = 2 s.

Table 3 :

 3 Fit ParametersMorris-Lecar Neurons (in mV )

		.1 3.2	10.9 -0.32 0.98	1.1	-1.2e-3	-1.4	3.9	-0.11
	FS-Cell	-51.2 1.8	-6.1 -0.86 1.6 -0.70	-11	-0.18	1.2	-1.2
	Cell Type	P 0	P µV

Acknowledgments

This research resulted from student projects during the Spring School of the European Institute of Theoretical Neuroscience (www.eitn.org), of which the participating students signed as co-first authors here. Research supported by the CNRS and the European Union (Hunan Brain Project H2020-720270 and H2020-785907). M.J. acknowledges support from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 616268 F-TRACT and the European Unions Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreement No. 785907 (Human Brain Project SGA2).


Michel Tallon
email: mtallon@obs.univ-lyon1.fr

Éric Thiébaut

Maud Langlois

Bernard Gelly

Richard Douet

Clémentine Béchet

Loïc Denis

Gil Moretto

The wavefront sensing making-of for THEMIS solar telescope

Keywords: adaptive optics, wavefront sensing, centroiding method, registration method, solar telescope

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Adaptive optics (AO) is spreading among solar telescopes. [START_REF] Rimmele | Solar Adaptive Optics[END_REF] In this context, an adaptive optics system is being installed on THEMIS * 90 cm solar telescope at Teide Observatory, Tenerife, Canaria Islands. This new developpement offers the opportunity to implement new methods of reconstruction and control, driven by the requirement that the AO runs unsupervised at its best performance in any atmospheric conditions. The AO system is combined with a major refurbishment of the telescope optics, from the telescope secondary down to the spectrograph entrance, in order to keep the unique spectropolarimetric capabilities of this telescope. [START_REF] Gelly | New life for the THEMIS solar telescope[END_REF] These changes in the optical path should lead to an AO corrected telescope fit for high quality polarimetry.

As shown by figure 2, the AO system is based on a deformable mirror with 97 actuators (11×11), from ALPAO, and a Shack-Hartmann wavefront sensor with 76 subapertures (10 × 10) with a Fried geometry 3 (Fig. 3). In order to keep the cost as low as possible, the Real Time Computer is a single PC with 4 cores at 4 GHz, and run the control loop at 1 kHz. This paper will focus on the new wavefront sensing method at the root of the control loop. The method aims at measuring optimally the displacements of the sub-images delivered by the wavefront sensor on the solar granulation seen with 2% contrast and a digitization on 8 bits only.

Since the day time atmospheric conditions may be quite instable, with sub-images randomly blurred in space and in time, some sub-images even fading away, the proposed method also estimate the errors corresponding to the measurements obtained at each frame, including their covariances. These errors will allow wavefront reconstruction to optimally take into account the variations of the errors in the pupil and will allow the loop to optimize itself to fast varying conditions.

Since the measurements now come with their associated precision matrix (invert of the covariance matrix) in real time, the computation of the commands for the adaptive optics correction must also be able to take them into account. The usual method relies on a command matrix obtained with a pseudo inverse of the interaction matrix. This inversion is usually computed with truncated SVD and may include modal filtering. [START_REF] Gendron | Astronomical adaptive optics. i. modal control optimization[END_REF] Such an estimator may include fixed weights, for instance to take laser guide star elongation into account. [START_REF] Tallon | Wavefront reconstruction with elongated spots[END_REF] But here, this scheme would implies to compute a new command matrix at each frame since precision matrix changes at this rate. The solution is to compute the commands by an iterative method, without matrices, [START_REF] Thiébaut | Fast minimum variance wavefront reconstruction for extremely large telescopes[END_REF] since it only uses the direct model (sparse interaction matrix, weights, and priors) that can be changed on the fly. The algorithm to compute the commands from the wavefront sensor measurements is out of the scope of this paper where we focus on wavefront sensing on open-loop real data acquired with THEMIS wavefront sensor.

The wavefront sensing method presented in the following is based on the fit of a reference image used as a model, on the pixels of the sub-images. The covariance matrix of the measurements can be obtained if the error on the pixel values is known. So the method starts with a specific calibration of the detector in order to estimate the error of each pixel value.

The next sections will overview successively the pre-processing of the pixels, with a particular step to compute an estimate of the error on each pixel value, the estimation of the slopes and their covariances from the preprocessed pixels, and the estimation and the update of the reference image used as a model for the estimation of the slopes.

PRE-PROCESSING OF RAW PIXELS

Pixel values are at the root of the whole processing. In order to estimate the covariance of the errors on the slope measurements, we need to estimate the variances of the pixel values. This implies an additional step to the usual dark / flat correction process, which is presented in this section. For each pixel i, we seek for the incident flux in arbitrary unit and the standard deviation of the errors on the estimated flux. In order to simplify the derivation of the equations, we assume here that all the necessary calibration frames (see after) are available with the same exposure time as the data (wavefront sensor frames). This assumption is quite practicable here since the flux from the Sun is fairly stable. With this assumption, a pixel raw value in ADU, r i , can be written as:

r i = t i φ i g i + b i + n i , (1)
where t i is the transmission from the sky down to the pixel i, including the quatum efficiency of the detector, φ i is the flux of the source in an arbitrary unit proportional to photons, g i is the gain (e -/ADU), and b i and n i are respectively the bias and the noise, both in ADU.

We use tree types of calibration frames obtained with the same exposure time:

• dark frames: φ i = 0, i.e. the pixel does not receive any flux. A large set of dark frames are recorded to compute the mean and variance on each pixel, denoted E(r dark i

) and Var(r dark i ) respectively.

• flat frames:

E(φ i ) = E(φ flat ), i.e.
we assume that all the pixels receive the same flux on average, even if this flux varies in time. This flat frame is obtained on the Sun while the telescope is speedily wandering other the solar granulation. We only need the mean of a large set of frames, denoted E(r flat i ). • static frames: Var(t i φ i ) = E(t i φ i ), i.e. we can assume Poisson statistics of the counts, which holds even if the pixels receive different amounts of light. In this case, we use an internal stabilized source that shed some, non-uniform but stable in time, light on the detector (a Mikrotron EoSens 4CXP camera). Here again, a large set of frames is used to compute the empirical mean and variance on each pixel, denoted E(r stat i ) and Var(r stat i ) respectively.

Compared to the usual detector calibration, we thus need the variance of the dark frames, and a set of socalled static frames. In principle, a source could deliver both flat and static frames, but the "flattening method" generally introduces fluctuations of the source, as it is the case with solar granulation.

The flat frames are obtained through the field-stop of the wavefront sensor; so their mean is also used to determine the regions of interest (ROI) on the detector for each sub-image. The flat correction also account for the non-uniform transmission of the optics (e.g. vignetting by the microlenses).

The calibrated value of the pixel i is obtained from a linear combination of raw value r i , as usual:

7 d i = α i (r i -β i ).
(

) 2
By assuming E(d i ) = E(φ i )/E(φ flat ), we get:

α i = 1 E(r flat i ) -E(r dark i ) , (3)
β i = E(r dark i ). (4)
Thus if the flat frame is obtained on the Sun itself, the calibrated values of the pixels are around unity.

Assuming that the value d i is an estimation of the expected value (we have a single sample for one pixel!), the variance of the pixel value is also obtained from a linear combination:

Var(d i ) ≈ max(d i , 0) + v i u i , (5)
with

g i = E(r stat i ) -E(r dark i ) Var(r stat i ) -Var(r dark i ) , (6)
u i = g i E(r flat i ) -E(r dark i ) , (7)
v i = g i Var(r dark i ) E(r flat i ) -E(r dark i ) . ( 8
)
Figure 3. Superimposition of a wavefront frame pre-processed as explained in Sec. 2, and the result of the measurement process as explained in Sec. 3. The 76 sub-images show a 10 field-of-view on 30 × 30 pixels, obtained with a ∼ 200µs exposure time at 1kHz. The contrast of the granulation seen in each subaperture is ∼ 2%. The crosses show the estimated locations of the sub-images for this frame, and the ellipses the estimated standard deviation of these locations. This frame shows that the accuracy along the x-axis is better than the one along the y-axis, and that subapertures at the edges of the pupil are less accurate. The reference image shown at the center is estimated from the previous frames on a larger field-of-view (dashed square) and with more details than the field-of-view of the subapertures (cf. Sec. 4).

The estimation of the errors on the pixels thus needs the additional computation of Var(r dark i

), E(r stat i ) and Var(r stat i ). In the following, we gather all the pixels d i of a sub-image k in a single vector d k . In the same way, their corresponding weights 1/Var(d i ) are gathered in the precision matrix W k of sub-image k. Since the covariance terms are neglected, the matrix W k is diagonal.

MEASUREMENT WITH A KNOWN REFERENCE IMAGE

Assuming Gaussian noise, maximum likelihood amounts to minimizing the following objective function:

ψ(θ) = n k=1 η k d k -α k R k (s k ).r 2 W k -m k log η k + µ D.r 2 2 (9)
with d k the calibrated pixel values in sub-image k, W k the corresponding precision matrix, m k the number of pixels in the sub-image k, R k (s k ) an operator that shifts the reference image r by an offset s k , α k a rescaling of the reference image for modeling the scintillation and η k a rescaling of W k to take modeling errors into account. Thus we aim here at fitting each sub-image d k with a known reference image, r, that can be shifted by an offset s k and rescaled by a factor α k . Because of registration errors, we expect that sub-images do not perfectly match the reference image, for that reason we allow the precision matrix to be rescaled by a quality factor η k for each sub-image. As the quality factors are part of the parameters to fit, the associated logarithm term must be kept in the objective function.

In Eq. ( 9), the term µ D.r 2 2 is meant to regularize the estimation of the reference image r that will be addressed in the next section. In this section, we assume the reference image r fixed and known, this term has no effect.

We are only looking for the offsets s k and their covariances, but need to fit in real time all the parameters θ = {η k , α k , s k } k=1,...,n ∪ r.

The operations applied to the reference image are non linear with respect to the shifts s k . A linearization of α k R k (s k ).r is necessary to derive an algorithm fast enough to be applied in real time. This linearization step relies on a first order expansion and a change of variable around an expected value of s k , as detailed in Thiébaut et al. [START_REF] Thiébaut | Innovative realtime processing for solar adaptive optics[END_REF] The likelihood term for one sub-image now writes:

ϕ k (u k ) = d k -H k .u k 2 W k = γ k -2 b T k .u k + u T k .A k .u k , (10)
with

u k = α k        1 x k y k        (11)
where (x k , y k ) is the sought relative displacement from the predicted position, H k is a m k × 3 matrix obtained from a first order expansion 8 of the shifted reference image, and:

γ k = d k 2 W k , (12)
b k = H T k .W k .d k , (13)
A k = H T k .W k .H k . (14)
The solution of Eq. ( 10) is:

u k = A -1 k .b k , (15)
and the corresponding value of ϕ k (u k ) at the minimum is:

ϕ k ( u k ) = γ k -b T k . u k (16)
The next step is to estimate the rescaling of the weights W k from Eq. (9). After this linearization and solving for u k for each sub-image k, Eq. (9) now writes:

ψ(θ) = n k=1 {η k ϕ k ( u k ) -m k log η k } , (17)
the value of the likelihood term, ϕ k ( u k ), being known from Eq. ( 16). Minimizing ψ(θ) for each sub-image, yields:

η k = m k ϕ k ( u k ) . (18)
For each sub-image k, the covariance matrix, C k , of the parameters (α k , x k , y k ) is approximately given by: 8

C k ≈ 1 η k J T k .A -1 k .J k , (19)
with J k the Jacobian matrix of the partial derivatives of the non-linear relationship between u k and (α k , x k , y k ).

Figure 3 shows an example of a processed frame. On top of a frame from the wavefront sensor, pre-processed as detailed in Sec. 2, the crosses show the estimated locations (x k , y k ) of the sub-images, and the ellipses denote the estimated standard deviation around these locations, from the coefficients of C k . The ellipses have been scaled so that the size of a sub-image correspond to the size of a pixel, so that the standard deviation of the position is ±0.5 pixel when the diameter of the an ellipse equals the size of a sub-image. We can see variations of the errors across the pupil, with larger errors from the subapertures at the edge of the pupil. The ellipse are also elongated because of the particular structure of the reference image shown in the middle of the image: the errors are not the same along x and y axis, with some slight correlation between the coordinates.

UPDATE OF THE REFERENCE IMAGE

In the previous section, by assuming the reference image r is known, the optimization of the general cost fonction ψ(θ) given by Eq. ( 9) allowed all the parameters but r to be determined in the set of the parameters θ = {η k , α k , s k } k=1,...,n ∪ r. In an alternate approach, we now look for the reference image r itself, assuming all the other parameters are known.

Denoting G k = α k R k (s k ) (i.
e. shifting the reference image by an offset s k , and rescaling by factor α k ), we can write the significant terms of Eq. (9) for our purpose as: (20)

ψ (θ) = n k=1 d k -G k .r 2 η k W k + µ D.
Thus the expression is reduced to regularized least-squares. The underlying approximation here for a fast computation when applying G k is to shift r by s k rounded at the nearest pixel, thus avoiding the need of any interpolation.

The regularization term is mandatory here because the reference image is used as the model to be fitted and need to be known everywhere it is needed, with a good accuracy. So the estimation of the reference image must be extrapolated on all the field-of-view, enlarged by the dynamic of the wavefront sensor. We assume that the enlarged field-of-view is twice the field-of-view of the sub-images. The chosen priors, D, is a difference operator that computes finite differences of the reference image.

Minimizing Eq. (20) for r yields:

n k=1 η k G T k .W k .G k + µ D T .D r = n k=1 η k G T k .W k .d k . ( 21
)
This equation is solved for r by using conjugate gradient method. We can notice that the right hand side expression is the weighted sum of the recentered sub-images.

In this process, the chicken or the egg problem is solved at the first frame by bootstrapping the method with setting r as the simple weighted average of all the sub-images. The first set of parameters are determined with this first reference image, up to the computation of a new reference image by solving Eq. ( 21). This new reference image is used for the next frame, and the first one is deleted.

For all the following frames, a new reference image is computed as soon as the measurements are sent to the controller that computes the commands. This new reference is merged with the previous one by using a leaking integrator.

By merging the recentered sub-images spatially in an optimal way, the reference image gets to be known in a field larger than the field-of-view of the sub-images. This allows all the field-of-view of the sub-images to be used to get the measurements of their shifts. An example of such an enlargement is shown on Fig. 3: the reference image is shown at the center of the image. This reference image is much less noisy than one sub-image chosen as a reference. Furthermore, the reference image is fixed and just evolves with the granulation pattern, thus allowing tip/tilt to be continuously measured while maintaining the same line of sight for a long time.

CONCLUSION

The presented wavefront sensing method for solar adaptive optics delivers for each frame, both the slope measurements and their covariance matrix. By using this information, the control loop can be robust to the fast spatial and temporal variations of the measurement errors. The quality of the measurements is improved by using an estimate of the reference image optimally built from the merging of all the sub-images, thus known with more details on a larger field-of-view and with much less noise. Furthermore, the slope measurements obtained by fitting such a reference image on the sub-images has been shown to be optimal. [START_REF] Thiébaut | Innovative realtime processing for solar adaptive optics[END_REF] The bootstrapping of the algorithm has been checked on open-loop data with various samples of the granulation.

At this time, the method is being implemented in the Real Time Computer in order to close the loop at 1 kHz. This goal is yet to be demonstrated on sky.

Figure 1 .

 1 Figure 1. Shack-Hartmann wavefront sensor installed at THEMIS telescope.

Figure 2 .

 2 Figure 2. Geometry of the THEMIS adaptive optics system, with a 10 × 10 Shack-Hartmann wavefront sensor (76 active subapertures) and a ALPAO DM97 deformable mirror with 97 actuators. The red circles denote the shape of the telescope pupil with its central obscuration.

  r

ACKNOWLEDGMENTS

This project has been co-funded by the European Commission's FP7 Capacities Programme under Grant Agreement number 312495, and the Centre National de la Recherche Scientifique.


Ingvald Erfjord
email: ingvald.erfjord@uia.no

Bozena Maj-Tatsis

Esther Levenson
email: levenso@tauex.tau.ac.il

Priska Sprenger
email: priska.sprenger@ph-karlsruhe.de

Marianna Tzekaki
email: tzekaki@auth.gr

Early Years Mathematics: Introduction to TWG 13

Introduction

The Early Years Mathematics working group shares scholarly research revolving around mathematics for children aged 2-8. This includes different areas of mathematics and attention to how mathematics is approached implicitly and explicitly by children and teachers alike. The broad age range affords attention to several transitions within preschool, and from preschool to the early grades of primary school. At each CERME-conference (since CERME 6), the working group on early years mathematics has spent time considering differences in Europe (and outside Europe) regarding the organization of preschool and the cultural approaches of different countries to education at this young age. Some countries have one preschool curriculum for children aged 1-5 (6), while others distinguish between nursery and preschool. The number of children attending preschool in different countries also differs greatly, with consequences for transitions. As an example, Norway has one preschool curriculum and 91.8 % (ages 1-5) and 97.1 % (ages 3-5) of children attended kindergarten in 20181 .

Prior to the conference in Utrecht, paper and poster authors were requested to prepare and send to the leaders a brief text outlining the main points of their contribution. This helped the leaders plan the discussions and presentations, and served as a basis for this introduction. Each of the 22 papers were allocated 30 minutes for presentation and discussion. For each paper, an author of a different paper was requested in advance to prepare questions and comments for the paper presenters. The seven posters were also briefly presented in a group session, followed by a group discussion in a one-hour session.

The number of participants in TWG13 in Utrecht was greater than the number of participants in previous years. In addition to participants from Europe, there were also stimulating contributions from outside of Europe (Japan, Canada, and Malawi). In total, 38 participants from thirteen different countries attended the working group (see Another branch of studies considers approaches that facilitate children's learning of mathematics during the early years. Grimeland and Sikko study children's development of mathematical literacy in a second-grade classroom working with play-coins in an inquiry-based setting. Inquiry-based settings are also the frame for the study by Skoumpourdi. However, her study is related to another area of mathematicsplane geometrical shapes. Skoumpourdi found that cooperative exploration, discussion, and justification processes were fostered in the inquiry-based setting. Pettersen and Volden look at the transition from mathematics in kindergarten to school mathematics. They argue that inquiry-based mathematics education was found to be helpful for children's ability to remember and recognize mathematical activities across the transition.

A study of Matsuo and Nakawa consider whether preschool children can improve their understanding of length and area measurement. In particular, their concern is on direct and indirect measurement, as well as measurement by an arbitrary unit. Their findings pinpoint that area measurement is found hardest to grasp, but partly in play-based settings, while direct comparison is easiest. Vanegas, Prat and Rubio characterize the learning trajectory of children aged six to eight years old when acquiring the notion of length measurement. The group also had a contribution by Bräuning that look at children's development of strategies for a combinatorial task and argues for a process-oriented view.

Palmér and Björklund focus on 2-3 years-old toddlers' ability to explore structural elements while playing "hide the toy-dragon", and pinpoint indications of emergent structural awareness among the small children. The inclusion of toddlers, as young as 2 years, within the scope of focus for TWG13 came as a result of the CERME10 conference and as an overall desire to increase insight into children's learning of mathematics [START_REF] Levenson | Early years mathematics[END_REF]. Nordemann and Rottmann look at how children's repeating pattern competencies can be fostered by physical activity, and report that development can be found especially for children with learning difficulties in mathematics. Nergård studies five-years-olds children's use of mathematical concepts, and found evidence of such use when they compare, explain and argue about natural phenomena.

Erfjord, Hundeland, and Carlsen characterise the mathematical discourse evolving as a teacher facilitates a mathematical activity for five-year-old children in kindergarten. Their findings indicate that mathematical discourse was characterised by engaging only one or two of the six children in longer periods and children's contributions were often of a non-mathematical nature. Gifford and Thouless address how teachers develop young children's pattern awareness.

Maffia and Mancarella consider children's interpretation of equalities represented through images, and found evidence that relational thinking can be supported when children work with number sentences containing equivalence. Maj-Tatsis and Tatsis look at characteristics of tasks that promote mathematical reasoning among young children. Tasks that required a shift from geometrical structures to numerical ones were found most difficult for children.

Methodological approaches

In TWG13, there has been an ongoing discussion regarding methodological approaches to investigating early mathematics learning. These approaches often reflect the culture of educational settings for children. In some countries in Europe, all kinds of testing and even organised activities are considered problematic when it comes to preschool. The role of play, from (total) children initiated free play to strong teacher guided play is another important discussion that has been ongoing in the group [START_REF] Levenson | Early years mathematics[END_REF]. In CERME 11, this discussion again arose and emphasised e.g., by Breive, who considered the balance between freedom and structure in mathematical activities for five-year-old children in kindergarten.

Most of the research reported in TWG13 follows a qualitative data collection approach, with interviews or observations, or lesson observations and video data. However, a variety of ways for analyzing the data were adopted such as: cross-case, process oriented, interpretative, and hermeneutical approaches. There are also reports from a few experiments or design research studies, with pre and post examination of children, mainly with interviews and qualitative analysis.

There are also a few quantitative approaches, utilizing tests or questionnaires, and one mixed study. These differences reflect the scope of that area of research, and that different kinds of research can contribute differently.

  table below). Lüken highlights the role of finger use during numerical considerations, and considers how the role and use changes over time. Sprenger and Benz consider children's ability to perceive and use structures in sets. They argue that eye-tracking, as a data collection instrument, together with observations, provide strong insight into the two processes of perception and determination of the cardinality of sets of objects. Bakos and Sinclair consider the role of a digital app on the development of children's multiplicative awareness. In particular, they consider the relation between bodily movement and mathematical meaning-making. Bakos and Sinclair found two primary gestures used by children which they discuss. Müller and Tiedemann

	Country	Number of papers	Number of posters	Number of Participants
	Canada	1		1
	Germany	5	3	10
	Greece	2		2
	Israel	1		1
	Italy	1		1
	Japan	1		1
	Malawi		2	2
	The Netherlands			1
	Norway	6		10
	Poland	1		1
	Sweden	3	1	5
	Spain	1		1
	UK		1	2
	TOTAL	22	7	38
	A variety of topics
	Several contributions to TWG13 concern classical, but indeed crucial areas in early years
	mathematics, such as numbers and arithmetic skills. Björklund and Runesson consider children's
	way of experiencing

numbers, and found evidence of six different ways children can experience numbers. Bjørnebye investigates children's ability to synchronize multiple representations of numerosity in the embodiment of a counting-on-strategy. He argues that bodily experiences could support early appropriation of a counting-on strategy. report from an ongoing study which focuses on children's development of conceptual subitizing, while Bruhn looks at first graders' creative processes when working on arithmetical open-ended problems. Tzekaki and Papadopoulou consider different levels of generalization that children can develop in arithmetic learning by systematically adopting teaching approaches with relevant tasks and discussion.

  Wernicke et al. look at social training of spatial perception and spatial cognition in an ongoing study based on a training program of children in a day-care centre.In two studies situated in Malawi, Gobede investigates mediation strategies used by early years mathematics teachers and Mandala explores how primary teacher education prepares pre-service teachers to teach early years mathematics. Svensson focuses on pre-school teacher students' ability to capture children's learning of mathematics in activities in preschool while Johansson, Tossavainen, Faarinen, and Tossavainen analyse student-teachers' perceptions regarding teaching mathematics for young children. Johansson et al. found evidence of different pedagogical beliefs and expectations regarding what mathematics teaching in preschool should be. Breive considers the balance between freedom and structure in kindergarten teachers' orchestration of mathematical learning activities. Tsamir, Tirosh, Levenson and Barkai investigated preschool teachers' selfefficacy for teaching patterning, with attention to why teachers report high self-efficacy before the professional program involved. Finally, Vogler focuses on the indirect processes of learning in preschool situations, which may be considered as an opportunity for children to participate and learn in what may be denoted as a "double layer structure' of every day meanings and abstract mathematical meaning".

Published on the website by The Norwegian Directorate for Education and Training https://www.udir.no/tall-ogforskning/statistikk/statistikk-barnehage/tall-og-analyse-av-barnehager-2018/barnehager/ (the website is currently written only in the Norwegian language)