Andrew P Davison

Sandra Diaz

Petruţ Antoniu Bogdan

Nathalie Azevedo Carvalho

Nathalie Azevedo Carvalho
email: nathalie.azevedo-carvalho@inria.fr

Sylvain Contassot-Vivier

Laure Buhry

Dominique Martinez

Simulation of Large Scale Neural Models With Event-Driven Connectivity Generation

Keywords: brain simulation, Hodgkin-Huxley neurons, time-stepping method, event-driven connectivity generation, Runge-Kutta method, Parkinson's disease, large scale networks

Accurate simulations of brain structures is a major problem in neuroscience. Many works are dedicated to design better models or to develop more efficient simulation schemes. In this paper, we propose a hybrid simulation scheme that combines time-stepping second-order integration of Hodgkin-Huxley (HH) type neurons with event-driven updating of the synaptic currents. As the HH model is a continuous model, there is no explicit spike events. Thus, in order to preserve the accuracy of the integration method, a spike detection algorithm is developed that accurately determines spike times. This approach allows us to regenerate the outgoing connections at each event, thereby avoiding the storage of the connectivity. Consequently, memory consumption is significantly reduced while preserving execution time and accuracy of the simulations, especially the spike times of detailed point neuron models. The efficiency of the method, implemented in the SiReNe software 1 , is demonstrated by the simulation of a striatum model which consists of more than 10 6 neurons and 10 8 synapses (each neuron has a fan-out of 504 post-synaptic neurons), under normal and Parkinson's conditions.

INTRODUCTION

Major projects such as Blue Brain [START_REF] Markram | The blue brain project[END_REF], the Human Brain Project [START_REF] Einevoll | The scientific case for brain simulations[END_REF], the BRAIN Initiative or Mindscope [START_REF] Hawrylycz | Inferring cortical function in the mouse visual system through largescale systems neuroscience[END_REF], aim at simulating a brain or brain structures. The simulation of brain structures demands not only computing resources but also a lot of memory to store the connectivity that grows as the power 1.4 of the number of neurons (see Figure 1). Then, the simulation of the human brain would require ≈ 100 terabytes just for storing the Boolean connectivity pattern (connection/no connection). Therefore, large-scale simulations of realistic cortical networks have been undertaken by using supercomputers with huge memory space [START_REF] Izhikevich | Large-scale model of mammalian thalamocortical systems[END_REF][START_REF] Migliore | Parallel network simulations with NEURON[END_REF][START_REF] Chatzikonstantis | Multinode implementation of an extended hodgkin-huxley simulator[END_REF][START_REF] Eliasmith | The use and abuse of large-scale brain models[END_REF].

Here, we propose a hybrid approach that combines a time-stepping approach for the numerical part of the simulation with an event-driven updating of the synaptic currents for complex Hodgkin-Huxley (HH) type neurons. This approach is particularly well-suited to studies of real-time neural mechanisms requiring high accuracy, such as pathological behaviors like the Parkinson's disease. Our event-driven approach completely avoids the storage of the connectivity pattern by regenerating the connectivity on the fly, when needed, after spiking events. This event-driven FIGURE 1 | Number of synapses N s vs. number of neurons N n for different animal species (vertebrates and invertebrates). The best fit is for the power law N s = N 1.4 n (red line) and not for N s = N n (black line), as claimed in [START_REF] Lansner | Virtues, pitfalls, and methodology of neuronal network modeling and simulations on supercomputers[END_REF], or the full connectivity N s = N 2 n (blue line). All fits are consistent in the sense that N s = 0 for N n = 0.

generation of the connectivity makes use of pseudo-random generators and consistent seeds. The exact computation of spike times is not possible for HH neurons because the model is continuous and the membrane voltage is approximated by timestepping methods on discrete time [START_REF] Koch | Methods in Neuronal Modeling: From Ions to Networks[END_REF].

Moreover, when the membrane potential crosses threshold twice during one time-step (the first crossing is upward and the second one is downward), the spike may be missed. A failure to detect a spiking event may cause dramatic changes on the behavior of the system, especially in the case of event-driven connectivity.

In this paper, we consider an event-driven connectivity generation within time-stepping schemes of Runge-Kutta 2 midpoint type. We develop a spike detection method for HH neurons, that accurately determines the spike timings so that the accuracy of the second-order Runge-Kutta methods (RK2) is preserved when connectivity is generated at spiking events. By avoiding the connectivity storage, our method is intended to simulate large-scale models made of Hodgkin-Huxley type neurons on a single computing node. Indeed, the limited memory consumption pushes back the necessity to use multiple machines, whose induced communications often reduce the overall performance. Yet, computing performance is not neglected in our approach as a parallel multi-threaded version has been developed in order to take advantage of multicore/many-core machines.

Our approach is implemented in the SiReNe software whose accuracy and efficiency are exhibited in a series of experiments among which is the simulation of the striatum structure at the rat scale, i.e., with more than 10 6 neurons.

In the next section, we demonstrate the originality and interest of our approach in respect to related works. The different methods used in our neural simulator are presented in section 3. Then, validation and performance analysis are provided in section 4. Section 5 draws a general conclusion and proposes a list of future works.

RELATED WORKS

Event-Driven Connectivity Generation

In neuroscience studies, one often simulates a snapshot of a network over a short period of time with fixed parameters, e.g., for comparing normal vs. pathological neural configurations. In such a context, the causal mechanisms (e.g., plasticity mechanisms) of neural evolution from normal to pathological states is not relevant. Thereby, connectivity storage is not mandatory and can be advantageously replaced by dynamic generation. Originally, the idea of event-driven connectivity generation has been proposed in the case of abstract neurons for which spike timing is exactly known, i.e., rule-based artificial cell units, or finite state machines [START_REF] Lytton | Rule-based firing for network simulations[END_REF]. This approach has then been applied with integrate-and-fire (IF) neurons, i.e., quadratic IF [START_REF] Izhikevich | Large-scale model of mammalian thalamocortical systems[END_REF] and leaky IF over GPU hardware [START_REF] Knight | Larger GPU-accelerated brain simulations with procedural connectivity[END_REF].

To the best of our knowledge, the event-driven connectivity generation approach has never been developed for Hodgkin-Huxley neurons. Our event-driven connectivity generation makes use of pseudo-random generators and consistent seeds. The principle of pseudo-randomly generating the connectivity through an event-driven approach has been reused in a recent US patent [START_REF] Lipasti | Efficient and Scalable Systems for Calculating Neural Network Connectivity in an Event-Driven Way[END_REF]. Yet, the LFSR (Linear-Feedback Shift Register) nature of the pseudo-random number generator (PRNG) used in their approach is not statistically strong as it does not pass classical linear tests provided in TestU01 [START_REF] L'ecuyer | Testu01: a c library for empirical testing of random number generators[END_REF] which is a reference in the domain. Moreover, a single PRNG is used to perform all the dynamic data generations, leading to additional time and memory consumption.

Spike Detection

Traditionally, event-driven strategies are applied at the neuron level [START_REF] Makino | A discrete-event neural network simulator for general neuron models[END_REF][START_REF] Brette | Exact simulation of integrate-and-fire models with synaptic conductances[END_REF][START_REF] Brette | Exact simulation of integrate-and-fire models with exponential currents[END_REF][START_REF] Tonnelier | Event-driven simulation of nonlinear integrate-and-fire neurons[END_REF][START_REF] Rochel | An event-driven framework for the simulation of networks of spiking neurons[END_REF][START_REF] Mattia | Efficient event-driven simulation of large networks of spiking neurons and dynamical synapses[END_REF][START_REF] Ros | Event-driven simulation scheme for spiking neural networks using lookup tables to characterize neuronal dynamics[END_REF][START_REF] Morrison | Exact subthreshold integration with continuous spike times in discrete-time neural network simulations[END_REF][START_REF] Rhodes | sPyNNaker: a software package for running PyNN simulations on SpiNNaker[END_REF]. In pure event-driven strategies the spike timings are analytically given and are calculated with an arbitrary precision (up to the machine precision). This scheme allows for an exact simulation where no spike is missed. Yet, only a limited class of simplified neuron models of integrate-and-fire (IF) type is amenable to exact simulations. Thus, more complex neuron models are simulated on discrete times by using time-stepping methods, e.g., the second order RK2 algorithm. Nevertheless, aligning the spike times on the time grid leads to an accuracy of order one, as for example in [START_REF] Lipasti | Efficient and Scalable Systems for Calculating Neural Network Connectivity in an Event-Driven Way[END_REF]. For IF neurons, the determination of the spike times by linear interpolation is needed to preserve the order 2 of the RK2 method [START_REF] Hansel | On numerical simulations of integrate-and-fire neural networks[END_REF][START_REF] Shelley | Efficient and accurate time-stepping schemes for integrate-and-fire neuronal networks[END_REF]. In contrast to IF neurons, HH neurons do not have an explicit threshold so that determining the spike times by threshold crossing, as in [START_REF] Lobb | Parallel event-driven neural network simulations using the hodgkin-huxley neuron model[END_REF], or by linear interpolation, as for IF neurons, will lead to low accuracy. Here, we show that for HH neurons, a quadratic interpolation (e.g., Bézier curves) is required to be consistent with the order 2 of the RK2 method. In contrast to [START_REF] Morrison | Exact subthreshold integration with continuous spike times in discrete-time neural network simulations[END_REF] in which quadratic and cubic interpolations are considered, it is not necessary to use more than a quadratic interpolation. Actually, cubic (and higher-order) interpolation should be avoided as it implies additional computation cost for no gain in accuracy [START_REF] Hansel | On numerical simulations of integrate-and-fire neural networks[END_REF].

Parallel Computation

Many studies have been done about the use of parallelism in neural simulations [START_REF] Kunkel | Meeting the memory challenges of brain-scale network simulation[END_REF][START_REF] Lansner | Virtues, pitfalls, and methodology of neuronal network modeling and simulations on supercomputers[END_REF][START_REF] Kunkel | Spiking network simulation code for petascale computers[END_REF][START_REF] Jordan | Extremely scalable spiking neuronal network simulation code: from laptops to exascale computers[END_REF]). Yet, most of them use connectivity storage, implying a huge memory consumption and the resort to distributed parallelism. By using event-driven connectivity generation, the memory requirements are significantly reduced and it becomes possible to simulate very large networks on a single computing node. In such a context, distributed parallelism is interesting merely for the increased computational power it offers, as compared to a single many-core node. However, it is worth noting that distributed parallelism implies additional overhead due to data communications between machines, which may significantly reduce the interest of the distribution. Consequently, we propose in this paper a multi-threaded parallel version of event-driven updating and connectivity generation that can run efficiently on a single multi-core node. Update its post-synaptic currents for time t + t

METHODS

In the following sections, we present the different types of simulation schemes: the common time-stepping approach, the event-driven connectivity generation, as well as the spike detection method we developed. Then, we describe the striatum model and use it to test the different methods. Finally, we present the implementation in the SiReNe software.

General Simulation: Time-Stepping

In the time-stepping approach, the state variables of the neurons are updated at each time-step ( t = 0.005 ms in our simulations). As the dynamics of the neurons is highly non-linear and sharply varies, the choice of the time-step must be thoroughly considered. In Figure 2, we observe that when the time-step is too large (e.g., 0.05 ms) the simulation results are incorrect. With a time-step of 0.03 ms, the potential curve begins to take the right shape, and it gets more accurate with smaller time-steps.

The ordinary differential equations of the Hodgkin-Huxley neurons type (HH) are solved by an explicit iterative temporal discretization method, i.e., Runge-Kutta 2 (RK2) in our case. The RK2 method is a second-order method leading to an error O t 2 of order two.

In classical time-stepping simulations, the synaptic current of each neuron is updated at each time-step. The pre-synaptic neurons are used at each step [START_REF] Lytton | Just in time connectivity for large spiking networks[END_REF] to perform the updating of the synaptic current. However, this approach is very time consuming and, most of the time, partly useless, as only spiking neurons generate post-synaptic currents. We propose in the next section a hybrid simulation scheme with event-driven post-synaptic updating that avoids useless computations.

Event-Driven Connectivity Generation and Post-synaptic Updating

In the hybrid simulation scheme that we propose, an eventdriven strategy is applied at the connectivity level. The idea is to be more pro-active by changing the order of the synaptic updating process. Instead of computing the state of the synaptic currents at time t according to previous time tt, the updating scheme executes the steps given in Algorithm 1.

The advantage of this method is to update the post-synaptic neurons of only the firing neurons, whereas in classical timestepping simulations, pre-synaptic neurons of all neurons are systematically updated, even those that have not fired recently. So, this event-driven updating scheme can significantly reduce the overall computations when only a fraction of the neurons fires in the same time-step. Moreover, this method can be combined to a pseudo-random generation of the connectivity, that significantly reduces the memory cost of the simulation.

Event-Driven Generation of the Connectivity

In order to avoid the storage of the entire connectivity of the neural system, and thereby to limit the memory consumption, the set of post-synaptic neurons as well as the intrinsic synaptic parameters (peak conductance,...) can be pseudorandomly generated after each spike [START_REF] Lytton | Just in time connectivity for large spiking networks[END_REF]. Synaptic parameters are defined for each synaptic model, and individuation can be obtained either by a specific computation based on the neuron number or by the addition of noise.

The pseudo-random nature of the generation ensures that the generated sequence is reproducible for a same seed. Thus, the connectivity of a neuron stays the same during the entire simulation as long as the same seed is used for that neuron. Also, the use of distinct seeds and internal PRNG states between the neurons implies that the connectivity is different from a neuron to one another. Consequently, when a neuron fires, its postsynaptic connections are pseudo-randomly generated according to its corresponding random seed and internal PRNG state, and the synaptic current of each post-synaptic neuron is updated.

In this context, the number of post-synaptic neurons of every neuron is fixed at the simulation initialization. This number can be provided by the user, either as an absolute value or through a connection density. Then, there is a need for an algorithm that uniformly draws M elements (post-synaptic neurons) from the integers set [1, N] (all possible candidates) without repetition. One practical solution is to sequentially parse the set of candidates and to perform a pseudo-random selection according to a probability defined in function of the number r c of remaining candidates and the number r s of elements still to be selected. Therefore, when considering candidate neuron i, the probability to select i is given by P(i) = r s r c . It can be checked that, by construction, this value is defined in [0, 1]. Also, the probability distribution of the selected connections obtained by this process is definitely uniform. In fact, for drawing the M elements, other solutions are possible that are theoretically more accurate in terms of probabilities. Indeed, the smallest value of k reals randomly drawn in the interval [x, 1] (0 ≤ x ≤ 1) is theoretically given by 1

-(1 -R) 1/k × (1 -x),
where R is a random real number in [0, 1]. Then, an additional step is required to obtain integer values. Obviously, such a method is much more computationally expensive than the one we use. Indeed, our choice is a good compromise between computation efficiency and distribution quality of the candidates selection.

In its current form, our simulator uses different types of neurons, such as excitatory vs. inhibitory neurons within the same neural group or distinct neurons belonging to different groups. A type of synapse is thus defined (with peak conductance and time constant parameters) between two types of neurons (source type and destination type). Consequently, the memory storage devoted to the synapses scales at most as the square of the number of neuron types and not as the square of the number of neurons. Also, a neuron in one group may be connected to any neurons in its own group and in the other groups. Its post-synaptic connections are generated by selecting for each post-synaptic group of the neuron, the number of post-synaptic neurons defined by the connection density between the neuron group and that particular post-synaptic group. Then, the parsing and pseudo-random selection process described above is applied independently to each post-synaptic group, using the pseudorandom draw presented in Algorithm 2.

In this algorithm, the function Seed() generates a unique seed from the neuron number. Although this is not mandatory, it may be designed to ensure a minimal distance between seeds Algorithm 2: Pseudo-random connectivity generation Input : NeuronNum is the neuron from which post-synaptic neurons are generated Input : NbNeuronsPostGroup is the total number of neurons in the post-synaptic group Input : NbToSelect is the number of post-synaptic neurons to select Output: List of selected neurons

1 State ← InitPRNG(Seed(NeuronNum)) 2 List ← ∅ 3 NbSelected ← 0 4 for i ← 0 to NbNeuronsPostGroup -1 do 5 r c ← NbNeuronsPostGroup -i 6 r s ← NbToSelect -NbSelected 7 if RandomReal(State) × r c < r s then // RandomReal() returns a value in [0,1] 8 List ← List ∪ {i} 9 NbSelected ← NbSelected + 1 10 end 11 end
of neurons with similar connectivity parameters, in order to avoid similarities in the distances between selected post-synaptic neurons. The function InitPRNG() initializes the internal state of the PRNG so that distinct states are produced from distinct seeds. Then, the internal state is used and updated in function RandomReal() that draws a random value in [0, 1]. Also, the list of selected post-synaptic neurons is not stored in memory in practice, but the corresponding synaptic currents are directly updated on the fly. When the post-synaptic group is the same as the neuron group, the neuron NeuronNum can be excluded from the list of post-synaptic neurons in order to avoid autapses.

In such a context, the loop in line 4 must be divided in two loops, according to the neuron indices that are lower or greater than NeuronNum.

Finally, in order to get a good distribution quality of the PRNG while preserving the performance constraint (generation speed), a fast and robust generator must be chosen. In our simulator, the Lehmer64 generator is used (derived from [START_REF] Lehmer | Mathematical methods in large-scale computing units[END_REF], as it is one of the fastest generator that passes the Big Crush of TestU01 (L'Ecuyer and Simard, 2007), a battery of statistical tests that is a reference for quality evaluation of PRNGs. Moreover, as this PRNG works on 64 bits, it provides 2 64 ≈ 10 19 seeds, which is much more than the number of neurons in the human brain (≈ 10 11 ).

Event-Driven Computation of Synaptic Currents

Time evolution of the synaptic conductance g jk between presynaptic neuron j and post-synaptic neuron k is modeled by the following differential equation

dg jk dt = - g jk τ jk (1)
with τ jk , the time constant of the synapse. Integrating the equation between t (current time) and t sp (spike time) leads to

g jk (t) = ḡjk exp - t -t sp τ jk (2)
with ḡjk = g jk t sp , the peak conductance. The total synaptic current I syn,k for neuron k is the sum of the contributions of the pre-synaptic neurons j and all their pre-synaptic spikes

I syn,k (t) = j g jk (t) V k -E jk = j ḡjk V k -E jk n j i=1 exp - t -t sp i τ jk Fact jk t spn j (3)
where E jk is the reversal potential of the synapse from neuron j to k and V k is the membrane potential of the neuron k. The number of spikes received from neuron j is denoted n j and t sp i represents the firing time of the i th spike. Fact jk is updated at each event (i.e., when neuron j emits a spike) as follows

n j ←-n j + 1 Fact jk ←-Fact jk + exp - t -t sp n j +1
τ jk (4)

Event Detection

Spikes are detected when the membrane potential exceeds a given threshold. In classical methods, the firing time is then aligned to the time-step. Yet, this trivial detection method leads to an error O ( t) that is not consistent with the RK2 method [START_REF] Hansel | On numerical simulations of integrate-and-fire neural networks[END_REF]. From a numerical point of view, the firing time should be at the maximum of the membrane potential. Herein, we propose two interpolation methods to obtain a more accurate estimation of the firing time. The former is based on the intersection between two linear interpolations and the latter is based on a Bézier curve.

As shown in section 4, only the latter leads to an error O t 2 that is consistent with the RK2 method.

Intersection Between Linear Interpolations

Here, the spike time t S is found at the intersection between two lines deduced from the membrane potential derivatives at the time-step frontiers (Figure 3A). In our case, the slopes of the two lines are defined by the derivatives dV 0 /dt and dV 1 /dt where V 0 and V 1 are the membrane potential at t and t + t, respectively. Note that when the time-step is small enough, the top value of the spike takes place inside a time-step whose derivatives at start and end times (t and t + t) are positive and negative, respectively. If we define b 0 = V 0 -dV 0 dt t and b 1 = V 1 -dV 1 dt (t + t), the equation of the two lines can be written as

y = dV i dt x + b i with i ∈ {0, 1} , ∀x ∈ [t, t + t] (5) FIGURE 3 | (A)
Spike time tS obtained by lines intersection issued from derivatives at t and t + t (full green line). (B) Spike time tS detected with the maximum of a quadratic Bézier curve (full green line) defined between t = t 0 and t + t = t 2 . In both figures, the reference time spike t S is obtained with 10 -7 simulation time-step (dashed red line).

where y ∈ R. The spike time estimate tS found at the intersection of the two lines is given by

tS = b 1 -b 0 dV 0 dt - dV 1 dt (6)
Although this method provides good results in some cases, it has not been selected for use in our simulator as it is not accurate enough in general. As a consequence, we have chosen a method of higher order that makes a more meaningful use of the three points given above (intersection and the two time-step boundaries), by interpreting the two derivatives at t and t + t as the tangents of a quadratic Bézier curve.

Bézier Curve

The Bézier curve of order two is a polynomial curve specified by three points P 0 , P 1 , and P 2 that define the tangents P 0 P 1 and P 1 P 2 , respectively at starting point P 0 and ending point P 2 . Its parametric form B (x) with x ∈ [0, 1] is given by

B (x) = (1 -x) 2 P 0 + 2 (1 -x) x P 1 + x 2 P 2 (7)
As illustrated in Figure 3B, the spike time tS is estimated from the interpolation of the membrane potential with a quadratic Bézier curve, the points P i = (t i , V i ), with i ∈ {0, 1, 2}, are defined at times t i and potentials V i . The end points of the curve, P 0 and P 2 , are given by the potentials V 0 and V 2 at times t 0 = t and t 2 = t + t. The point P 1 is defined as the intersection of the two lines following the derivatives in P 0 and P 1 , similarly to section 3.3.1. The time t 1 is determined as tS in Equation ( 6), and V 1 is deduced from Equation ( 5) and ( 6), as

V 1 = b 2 dV 0 dt -b 0 dV 2 dt dV 0 dt - dV 2 dt (8)
Once P 0 , P 1 , and P 2 are obtained, tS is analytically computed as the time at which the Bézier curve reaches a null derivative (curve peak), i.e.,

tS = 1 -x 2 t 0 + 2 x 1 -x t 1 + x2 t 2 (9) where x ∈ [0, 1] is defined such that dB dx (x) = 0, leading to x = V 0 -V 1 V 0 -2V 1 + V 2
. The particular case where V 0 -2V 1 +

V 2 = 0 does not occur in practice as it would mean that the spike extremum is constant over the time step. Indeed, the Bézier interpolation is performed only when the derivatives at the time step bounds are not null.

Simulated Striatum Model

The striatum is a brain structure presumably involved in generating pathological β-oscillations observed in Parkinson's disease [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF][START_REF] Corbit | Pallidostriatal projections promote β oscillations in a dopamine-depleted biophysical network model[END_REF]. The striatum is composed in its vast majority (≈ 95%) of Medium Spiny Neurons (MSN) [START_REF] Kemp | The structure of the caudate nucleus of the cat: light and electron microscopy[END_REF][START_REF] Corbit | Pallidostriatal projections promote β oscillations in a dopamine-depleted biophysical network model[END_REF].

Our MSN model equations and the ionic channel are derived from the [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF] model.

The voltage V at each time-step is described by:

C m dV dt = -ḡK n 4 (V -E K ) I K -ḡNa m 3 h (V -E Na ) I Na -ḡM p (V -E K ) I M -ḡ L (V -E L ) -I syn + I app (10)
where C m is the membrane capacitance, ḡX is the maximal conductance of ion X. The values h, m, n, and p are the respective gating variables of the different ions channels (activation/inactivation). Also, E X is the reversal potential of the ion channel X. The fast potassium current I K has four activation gates and no inactivation gate. The sodium current I Na has three activation gates and one inactivation gate. M-current is a non inactivating potassium current which has one activation gate and no inactivation gate. The leak current is denoted by I L .

The synaptic current I syn is a GABAa inhibitory current given in Equation ( 3) with E jk = -80 mV and ḡjk is between 0.1 N mS and 0.6 N mS, where N is the number of synaptic connections. The applied current I app is a step current, i.e., at the beginning it is -10 µ A and at t = 500 ms it goes up to 0.4 µ A. We also add a uniform noise between -b and b to the applied current.

Since the idea of our article is to promote the simulation of large-scale networks, we focus on the simulation of the rat's MSN network at scale one. This network is composed of about 1.3 million neurons [START_REF] Oorschot | Total number of neurons in the neostriatal, pallidal, subthalamic, and substantia nigral nuclei of the rat basal ganglia: a stereological study using the cavalieri and optical disector methods[END_REF]. We consider that a MSN neuron is surrounded by ≈ 2800 other MSN neurons and that the connection density in this neighborhood is of 18%, as in [START_REF] Taverna | Recurrent collateral connections of striatal medium spiny neurons are disrupted in models of parkinson's disease[END_REF] and Lindahl and Hellgren Kotaleski (2017)), i.e., the number of post-synaptic neurons per MSN would be ≈ 504 neurons.

We thus simulate a MSN model of 1.3 million neurons with 504 connections per neuron. To demonstrate the efficiency of our simulation approach, we compare our results to [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF], during a 4 s biological time simulation. In all our simulations, the time-step is fixed to t = 0.005 ms and the Bézier's curve is used to interpolate the spike times. The Local Field Potential (LFP) is the sum of the synaptic currents of each neuron at each time-step [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF]. The LFP power spectrum was deduced from a standard Fourier Transform Function.

Further tests are done to compare classical time-stepping and event-driven updating approach, in terms of memory consumption and execution times. In this context, we simulate 100, 500, 1 000, 5 000 and 10 000 neurons for 100 ms biological time.

Finally, a last experiment is done to evaluate the performance of the parallel multi-threaded version with OpenMP. We compare execution times for a 1 s biological simulation time with 10 000 MSN neurons as a function of the number of threads (1, 2, 4, 8, 16, and 32).

Software Implementation in SiReNe

SiReNe ("Simulateur de Réseaux de Neurones, " "Neural Network Simulator" in English) software, is a C program that has been developed in our laboratory for several years. Originally, it is a pure time-stepping simulator, but it has been extended with the event-driven updating method presented in this paper. So, it is now able to simulate large neural networks either with a classical time-stepping method, or with the event-driven one. In both approaches, the non-linear differential equations of the HH model are updated with the help of the RK2 integration method. After updating all state variables, it is checked whether the neurons have fired. The potential of the neuron has to cross a threshold during this time-step and the derivatives of the potential at times t and t + t must be respectively positive and negative. In order to detect the spikes, the time delay of each neuron has to be larger than the time-step, otherwise some spikes may be missed. In the pure time-stepping version, the presynaptic currents are systematically retrieved for every neuron, whereas in the event-driven version, neurons are updated first, then only the post-synaptic connections of the spiking neurons are generated and their currents are updated. Figure 4 shows the general algorithmic scheme of the SiReNe software.

In addition, two variants were implemented, with or without storage of the neurons connectivity. The most efficient combination of time-stepping approach and connectivity management, according to execution time as well as memory consumption, is the event-driven (spike-driven) approach with connectivity generation.

Finally, parallel multi-threading has been added via the OpenMP API. The parallel strategy mainly consists in distributing the computations related to the neurons over the threads. In particular, parallelism is exploited in the computation of the derivatives, the spike detections and the synaptic current updates. For the parallel generation of neuron connections, it is mandatory to use a thread-safe random generator. Thus, we have added the required software layer over the random generator described in section 3.2.1.

Experimental Context

The simulations presented below have been done on a Dell R720 server under Linux Debian 4.9 amd64 with 2 Intel(R) Xeon(R) CPU E5-2640 v2 @ 2.00 GHz with 8 cores each, and 128GB RAM. Times are measured with the OpenMP function omp_get_wtime() and the program is compiled with gcc 6.3.0 and the -O3 optimization level. All simulations presented in this section are done in sequential (one thread) save for the ones related to the performance evaluation of the parallel version (multiple threads).

RESULTS

Order of the Method

The order of the simulation method is estimated on the last spike of a 20 ms simulation fired by a single neuron over repeated trials (N t = 10) with random initialization. We simulate different time-steps between t = 5e -5 ms and t = 0.01 ms. Three different methods are evaluated. In addition to the lines intersection (RK2Lines) and the maximum of Bézier curve (RK2Bezier) described in section 3.3, we add the simplest method in which the spike is aligned to the time-step (RK2Threshold) when the threshold is exceeded. A complete list of the simulation parameters related to this experiment (named Method order) is given in Supplementary Material.

Although the RK2 numerical scheme is used, the simulation order is not systematically two, as it depends on the spike detection method [START_REF] Hansel | On numerical simulations of integrate-and-fire neural networks[END_REF]. As there is no analytical solution to calculate the error of our model, we use as reference a simulation with a very small time-step t = 1e -7 ms. Indeed, for such small time-step, the error can be considered negligible. Thereby, in this context the error over the spike time (ǫ t ) is defined as

ǫ t = 1 N t N t i=0 | ti -t i | (11)
where N t is the number of trials, ti is the time of the last spike of the neuron within the simulated period for a given t, and t i is the corresponding time for the reference simulation with t = 1e -7 ms.

The error ǫ t over the spike time is depicted in Figure 5 for the three interpolation methods as a function of the time-step t (log scales). For each method, linear regression provides a line whose slope represents its order. Indeed, RK2Threshold and RK2Lines are first order methods (error is O ( t)). The RK2Bezier detection method is more accurate and preserves the order two of the original RK2 integration method (error is O t 2 ).

Performance Comparison of the Different Methods

As described before, two distinct simulation approaches are available in SiReNe: the classical time-stepping approach and the event-driven updating one. In addition, it is possible either to store the connectivity or to generate it on-demand. When combining those two aspects, we obtain four algorithmic variants. In order to determine which of the four approaches is the most efficient, we compare them in terms of memory consumption and execution time. Also, in order to compare our simulator to a reference from the community, we add simulations obtained with BRIAN2, which is one of the most used simulator for spiking neurons [START_REF] Goodman | Brian: a simulator for spiking neural networks in python[END_REF].

In Figure 6 (resp. Figure 7), memory consumption (resp. execution time), are given as a function of the size of a MSN network with 100% of connectivity (A) or 30% of connectivity (B) during 100 ms biological time. The five compared approaches are the classical Time-Stepping approach with connectivity Storage (TS-S), the same approach with connectivity Generation (TS-G), the Event-Driven updating with connectivity Storage (ED-S), the same approach with connectivity Generation (ED-G), and the approach used in the BRIAN2 simulator. A complete list of the simulation parameters related to this experiment (named Memory and time) is given in Supplementary Material.

From these results, it appears that storing the connectivity is definitely inappropriate, as it sharply increases the memory consumption (Figure 6). As expected, the method with storage has a memory consumption that scales in a polynomial way with the number of neurons while the scaling is linear for the method without storage. Concerning the performance, we observe a relationship between the execution time and the number of neurons. The time varies as a power of the number of neurons that is t ∝ n k , where t is the execution time, n the number of neurons and k the exponent of the power law. In Figure 7, the slopes of the ED-S, ED-G methods and BRIAN2 are of order one (k = 1), whereas for the other two methods the slopes are of order two (k = 2). It means that the execution time increases linearly for ED-S, ED-G, and BRIAN2 and quadratically for the two other methods. Nevertheless, BRIAN2 has a significantly larger simulation time compared to the EDmethods. The ED-S and ED-G methods have similar simulation times. However, for full connectivity (Figure 7A), we observe a slight gain in execution time with the non-storing version that comes from the possibility to avoid the on-demand generation in that particular case. When the connection probability is less than 1 (Figure 7B), the connectivity generation induces a slight overhead on the execution time. However, we observe that this overhead tends to decrease when the number of neurons increases.

For a last comparison between SiReNe and Brian, we implemented in SiReNe the COBAHH network model described in article [START_REF] Brette | Exact simulation of integrate-and-fire models with exponential currents[END_REF]. The COBAHH model is a network of 4,000 excitatory-inhibitory neurons (80% excitatory and 20% inhibitory Hodgkin-Huxley-type neurons with full connectivity). In [START_REF] Brette | Exact simulation of integrate-and-fire models with exponential currents[END_REF], this benchmark model was simulated in Brian with Euler integration (0.01 ms step-size) and the spikes were detected by threshold crossing at -20 mV with 3 ms refractory period and alignment of the crossing events on the step-size. Although this detection method associated with Euler integration leads to a precision of order one, we implemented them in SiReNe for a fair comparison with Brian. We simulated this benchmark model for 1 s of biological time on the same machine (Dell R720 server, section 3.6) for both Brian and SiReNe. The simulation time is more or less the same, ∼ 1 min. The memory consumption was 18 MB for SiReNe vs. 375 MB for BRIAN.

Parallel Computing Performance

In order to evaluate the performance of the parallel multithreaded version of SiReNe, the memory consumption and execution times are measured as a function of the number of threads for the best approach (event-driven). The test case is the simulation of 10 000 MSN neurons with 504 incoming connections each, during a biological time of 1 s. A complete list of the simulation parameters related to this experiment (named Parallel performance) is given in Supplementary Material.

As can be seen in Figure 8A, although the memory consumption increases slightly with the number of threads, it stays very limited. Indeed, the additional consumption with 32 threads compared to 1 thread is smaller than 3 Mb (less than 12 % of the initial consumption). This comes from the fact that only the conductance factors between neurons (see Equation 4) are duplicated to support their concurrent updates by the threads. Concerning the execution time, we observe in Figure 8B a significant decrease when the number of threads increases. This result shows that the parallel version provides a significant gain of time. However, the time decrease is less important than expected, leading to moderate speed-ups for large numbers of threads (≥ 16), and a faster decrease of the parallel efficiency (speedup over the number of threads) than expected (see Figure 9). This is due to the irregular inter-dependencies between neurons, that induce irregular memory access patterns, penalizing the parallel accesses. This issue deserves a complete detailed study that is planed as a future work.

Simulation of the Striatum Model

Before going on to the simulation of the striatum, we simulate three different sizes of neural networks (500, 1,000, and 2,000 neurons) and we compare the average number of spikes per neuron and the simulation time as a function of connectivity. A complete list of the simulation parameters related to this experiment (named Network comparison) is given in Supplementary Material. In Figure 10, we can see that for a sparse connectivity, the neural network has an asynchronous regime and for a dense connectivity, the neural network has a synchronous regime. For the respective neural network sizes, abrupt transitions of the number of spikes can be observed at different percentages of connectivity. The simulation times are quite stable when the connectivity varies.

For the comparison of the results with [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF], two types of plots are produced: the raster plot of the spiking neurons and the power spectral density of the LFP model. The former displays the set of spiking neurons at each time step. The latter draws the LFP power as a function of the spiking frequency.

The test case is a MSN network with 1.3 million neurons, simulated during 4 s biological time, either under normal or Parkinson conditions. The simulation time is 10 h per second of biological time with a time-step of 0.01 ms. The memory allocation is less than 1 GB. The normal conditions correspond to an applied input I app = 1.19 µA and a maximal Mcurrent conductance ḡM = 1.34 ms. To obtain Parkinson conditions, ḡM is decreased to 1.1 ms. A complete list of the simulation parameters related to this experiment (named Reference (striatum)) is given in Supplementary Material.

The raster plot in Figure 11A shows a MSN simulation under normal conditions. The average spiking rate for the neurons is 0.94 ± 0.63 Hz, which is coherent with the average MSN spiking rate in vivo (1.1 ± 0.18 Hz, [START_REF] Kish | Multiple single-unit recordings in the striatum of freely moving animals: effects of apomorphine and d-amphetamine in normal and unilateral 6-hydroxydopamine-lesioned rats[END_REF]0.96 ± 0.03 Hz, McCarthy et al., 2011). Figure 11B shows the simulation of a parkinsonian striatum. As expected, there is a pathological synchronization of the neurons, which is representative of Parkinson's disease. The average spiking rate raises up to 3.71 ± 1.03 Hz (2.11 ± 0.43 Hz in [START_REF] Kish | Multiple single-unit recordings in the striatum of freely moving animals: effects of apomorphine and d-amphetamine in normal and unilateral 6-hydroxydopamine-lesioned rats[END_REF]4.9 ± 0.15 Hz in McCarthy et al., 2011). The mean firing rates are significantly different between healthy and Parkinson conditions (p < 1e -3, t-test). Concerning the LFP spectrum, it can be seen that we obtain similar qualitative behaviors than those reported in [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF] under normal and parkinsonian conditions. In Figure 12A (normal conditions), a small peak can be observed (max of 30 dB at 17 Hz), revealing a weak synchronized state. In 2011). It has to be mentioned that no comparison could be done with Brian2 for such a large network, as it exceeds the memory capacity of our server (128 GB RAM) to store the connectivity.

CONCLUSION

An event-driven updating approach for the simulation of neural networks has been presented. It mainly consists in updating post-synaptic currents only for spiking neurons, thus reducing the computational cost. Moreover, the addition of a pseudorandom generation of the neurons connectivity allows for a significant reduction of memory cost, passing from O(n 2 ) to O(n) complexity, for n neurons in the system. Finally, two interpolation methods (linear intersection and Bézier curve) have been compared for the computation of spike times.

A set of experiments has been conducted to evaluate the validity of the approach as well as its accuracy and performance. The validity is confirmed by a comparison to the reference work of [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF] modeling a basal ganglia MSN network in healthy and Parkinsonian conditions. Concerning the accuracy, we conclude that to preserve the order 2 of the Runge-Kutta method, it is mandatory to use an order 2 interpolation method of the spike times. In our case, we have opted for quadratic Bézier curves as they are controlled by the tangents at end points and they have limited computational cost. Finally, the performance study is 2-fold. It shows that the eventdriven updating approach is not only faster than the classical time-stepping one, but that it has a better performance scaling according to the number of neurons. Also, the connectivity generation seems to have a limited impact on the performance, especially when the number of neurons increases.

One limitation in the event-driven connectivity generation concerns the inclusion of plasticity rules, such as STDP (Spike-Time Dependent Plasticity). Indeed, such rules are difficult to implement in an event-driven strategy because the synaptic conductances have to be stored in order to be updated at each spike time, which would cause a large increase in memory consumption. However, in many neuroscience studies, the plasticity is not a crucial key-point as a snapshot of a neural network is simulated over a short period of time in order to analyze its functional behavior rather than its evolution. This is true in the simulations of our Striatum network, as pathological oscillations in the Parkinson state are not obtained from synaptic plasticity but by changing the intrinsic M-channel conductance of the MSN neurons [START_REF] Mccarthy | Striatal origin of the pathologic beta oscillations in parkinson's disease[END_REF]. Other limitations of the current version of SiReNe, like the event-driven generation of distance-dependent connectivity or synaptic delays, are less difficult to tackle and should be addressed with a limited impact over performance.

Another contribution concerns the overall numerical methods used to simulate Hodgkin-Huxley based neuron models. This type of neurons is sensitive to the time-step value because with large time-steps (> 10 -2 ms), the significance of the derivatives used to deduce next step values are much weaker than with small time-steps (10 -7 ms) [START_REF] Moore | On numerical integration of the hodgkin and huxley equations for a membrane action potential[END_REF]. In fact, when there is no sharp variation of the membrane potential inside a time-step, the derivative at the beginning is close to the average value of derivatives inside the step, leading to a good approximation at the end of the step. However, when the potential sharply varies, typically when the neuron is spiking, the derivative at the beginning of the step is much lower than the average value of the derivatives inside the step, leading to an underestimation of the potential at the end of the step. In order to stretch the limits induced by this problem, a specific process is under consideration for inclusion in the SiReNe software. Although a complete study is necessary to fully evaluate the gain and interest of such corrective process, preliminary experiments show promising results. A future work will be dedicated to this subject.

In this paper, the main features of the event-driven updating and connectivity generation have been studied. The combination of this approach with time-stepping numerical integration of the equations forms a very pertinent solution to efficiently and accurately simulate large neural networks with limited computing resources (single server). However, there is still room for improvements and extensions. As a few examples, the reduction of irregular memory accesses in the multi-threaded version, extensions of the parallel version to support multiple machines and/or GPUs, the inclusion of distance-dependent connectivity, synaptic delays as well as the corrective process of the numerical method when using large time-steps would deserve further studies. All those aspects will be considered as future works.

FIGURE 2 |

 2 FIGURE 2 | Time-stepping simulation of a single MSN neuron with different time-steps t = 0.05, 0.03, 0.01 and 0.005 ms (plot of the first action potential).

FIGURE 4 |

 4 FIGURE 4 | Schematic representation of the SiReNe software (implemented in C).

FIGURE 5 |

 5 FIGURE 5 | Error ǫ t over spike time in function of time-step t, for the three spike detection methods.

FIGURE 6 |

 6 FIGURE 6 | Memory consumption as a function of the size of a MSN network with 100% of connectivity (A) or with 30% connectivity (B) during 100 ms biological simulation, for the TS-S, TS-G, ED-S, and ED-G simulation approaches and the BRIAN2 simulator.

FIGURE 7 |

 7 FIGURE 7 | Execution time as a function of the size of a MSN network with 100% of connectivity (A) or with 30% connectivity (B) during 100 ms biological simulation, for the TS-S, TS-G, ED-S, and ED-G simulation approaches and the BRIAN2 simulator.

FIGURE 8 |

 8 FIGURE 8 | Memory consumption (A) and execution time (B) as a function of the number of threads for the simulation of a MSN network with 10 k neurons (504 incoming connections per neuron) during 1 s biological time.

FIGURE 9 |

 9 FIGURE 9 | Speedup and efficiency of a 1 s biological time simulation of 10 k MSN neurons with 504 incoming connections, as a function of the number of threads.

Figure 12B (

 12B Parkinson's conditions), the LFP reaches higher β band oscillations (≈80 dB at 23 Hz), which is representative of Parkinson's pathological oscillations with strong neural synchronization (McCarthy et al., 2011).

FIGURE 10 |

 10 FIGURE 10 | Three different neural networks are simulated, i.e., 500, 1,000, and 2,000 neurons. The average number of spikes per neuron and simulation time are compared as function of the connectivity. For a sparse connectivity, we observe an asynchronous regime [raster plot indicated with (AS)]. For dense connectivity, we observe synchronous regime [raster plot indicated with (S)]. An abrupt transition occurs at a certain percentage of connectivity.

FIGURE 11 |

 11 FIGURE 11 | Raster plot of the large MSN test case for neurons in interval [100,000 : 100,100]. Spiking neurons at each time-steps under normal conditions (A) and Parkinson's conditions (B) for the simulation of a MSN network at scale 1 of the rat, during 4 s biological time.

FIGURE 12 |

 12 FIGURE 12 | LFP power in function of the frequency under normal conditions (A) and Parkinson's conditions (B) for the simulation of a MSN network at scale 1 of the rat, during 4 s biological time.

Source code is available at: https://sirene.gitlabpages.inria.fr/sirene.

Frontiers in Neuroinformatics | www.frontiersin.org

October 2020 | Volume 14 | Article 522000

DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS

NA, SC-V, LB, and DM contributed to the conception of the simulator. NA, SC-V, and DM contributed to the development of the Sirene software. NA conducted the experiments to validate the models. All authors contributed to manuscript revision, read, and approved the submitted version.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf. 2020.522000/full#supplementary-material

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Yifan Zhu
email: yifan.zhu@univ-lorraine.fr

Shi-Wang Fan

Liyun Cao

Krupali Donda

Badreddine Assouar
email: badreddine.assouar@univ-lorraine.fr

Acoustic meta-equalizer

Keywords: Acoustics, Metamaterials, Equalizer, Tunable filter, Ultra-broadband feature

Acoustic equalization is the process of adjusting the frequency response of a broadband sound signal, which is widely used in communication system and audio acoustics. Here, we introduce the unique concept of acoustic meta-equalizer, viz., a passive acoustic metamaterial-based filter that has tunable frequency response within ultra-broadband (200Hz -20000Hz), and capable to generate a frequency equalization on the input sound signal. By exploiting two kinds of acoustic resonant elements, viz., Helmholtz and Fabry-Perot resonators, we realize and demonstrate ultra-broadband and tunable filters, within an integrated meta-structure. In analogy with the conventional equalizer concept in signal processing, we numerically, analytically, and experimentally demonstrate functional filters, signal reproductions, and sound effect controls by the conceived acoustic meta-equalizer.

I. INTRODUCTION

Equalization is an important concept in signal processing, that is, adjusting the frequency response of a broadband signal. It has been and still widely used in telecommunications [START_REF] Qureshi | Adaptive equalization[END_REF][START_REF] Falconer | Frequency domain equalization for single-carrier broadband wireless systems[END_REF][3], and various acoustic applications including audio device [3][START_REF] Everest | The master handbook of acoustics[END_REF][START_REF] Howard | Acoustics and Psychoacoustics[END_REF], room acoustics [START_REF] Radlovic | Equalization in an acoustic reverberant environment: Robustness results[END_REF][START_REF] Talantzis | Robustness of multichannel equalization in an acoustic reverberant environment[END_REF][START_REF] Spors | Active listening room compensation for massive multichannel sound reproduction systems using wave-domain adaptive filtering[END_REF][START_REF] Kuo | Active noise control systems[END_REF][START_REF] Bouchard | Multichannel affine and fast affine projection algorithms for active noise control and acoustic equalization systems[END_REF] and underwater communications [START_REF] Stojanovic | Adaptive multichannel combining and equalization for underwater acoustic communications[END_REF][START_REF] Chitre | Underwater acoustic communications and networking: Recent advances and future challenges[END_REF]. The resulting filter device is known as equalizer (EQ). The acoustic equalizer device can adjust the loudness curve of a broadband sound signal to improve the quality of audio sound, because the tone color of sound is influenced by different frequency ingredients [START_REF] Everest | The master handbook of acoustics[END_REF][START_REF] Howard | Acoustics and Psychoacoustics[END_REF]. In room acoustics [START_REF] Radlovic | Equalization in an acoustic reverberant environment: Robustness results[END_REF][START_REF] Talantzis | Robustness of multichannel equalization in an acoustic reverberant environment[END_REF][START_REF] Spors | Active listening room compensation for massive multichannel sound reproduction systems using wave-domain adaptive filtering[END_REF], different complex acoustic environments will influence the uniformity of the frequency response, so that acoustic equalization can be used to reproduce the sound signals [START_REF] Radlovic | Equalization in an acoustic reverberant environment: Robustness results[END_REF][START_REF] Talantzis | Robustness of multichannel equalization in an acoustic reverberant environment[END_REF][START_REF] Spors | Active listening room compensation for massive multichannel sound reproduction systems using wave-domain adaptive filtering[END_REF], or to realize active noise control [START_REF] Kuo | Active noise control systems[END_REF][START_REF] Bouchard | Multichannel affine and fast affine projection algorithms for active noise control and acoustic equalization systems[END_REF]. In hydroacoustics, similarly, the complex environment makes the acoustic equalization process necessary in underwater acoustic communications [START_REF] Stojanovic | Adaptive multichannel combining and equalization for underwater acoustic communications[END_REF][START_REF] Chitre | Underwater acoustic communications and networking: Recent advances and future challenges[END_REF].

On the other hand, during the last decade, acoustic metamaterials have been actively studied [START_REF] Ma | Acoustic metamaterials: From local resonances to broad horizons[END_REF][START_REF] Cummer | Controlling sound with acoustic metamaterials[END_REF], providing different design strategy for various acoustic devices. For example, acoustic metamaterials have been proposed to actively control specific interesting acoustic phenomena such as cocktail effect [START_REF] Xie | Single-sensor multispeaker listening with acoustic metamaterials[END_REF], acoustic reverberation [START_REF] Ma | Shaping reverberating sound fields with an actively tunable metasurface[END_REF], and acoustic hologram [START_REF] Zhu | Fine manipulation of sound via lossy metamaterials with independent and arbitrary reflection amplitude and phase[END_REF]. Acoustic metamaterials are also exploited for underwater communications [START_REF] Shi | High-speed acoustic communication by multiplexing orbital angular momentum[END_REF][START_REF] Jiang | Twisted acoustics: metasurface-enabled multiplexing and demultiplexing[END_REF]. For the designs of sound absorber, acoustic metamaterials can bring a real added value to the absorption performance with large bandwidth and reduced sample thickness [START_REF] Yang | Optimal sound-absorbing structures[END_REF]. To further decrease the thickness of acoustic resonant units, the concept of acoustic metasurfaces was proposed [START_REF] Assouar | Acoustic metasurfaces[END_REF] to design unconventional ultra-thin acoustic devices. For instance, a metasurface-based Schroeder diffuser design [START_REF] Zhu | Ultrathin acoustic metasurface-based Schroeder diffuser[END_REF] is only 1/10 thickness of the conventional ones [START_REF] Cox | Acoustic absorbers and diffusers: theory, design and application[END_REF]. For sound absorber, deep-subwavelength metasurface-based absorber was evidenced [START_REF] Li | Acoustic metasurface-based perfect absorber with deep subwavelength thickness[END_REF], which is much thinner than the absorber from classical materials, showing the great flexibilities and potentials of acoustic metamaterials and metasurfaces.

In this research, we extend the general concept of acoustic metamaterials to equalizer, by proposing the unique concept of acoustic meta-equalizer (AMEQ), viz., a tunable and ultra-broadband passive acoustic metamaterial filter. Previous designs of acoustic filters [START_REF] Ma | Low-frequency narrow-band acoustic filter with large orifice[END_REF][START_REF] Liu | An acoustic dual filter in the audio frequencies with two local resonant systems[END_REF][START_REF] Assouar | Broadband plate-type acoustic metamaterial for low-frequency sound attenuation[END_REF][START_REF] Song | Broadband fractal acoustic metamaterials for low-frequency sound attenuation[END_REF] have shown limited bandwidth (narrow band) or limited amplitude values (only 0 and 1).

Therefore, an ultra-broadband and tunable acoustic filter is highly desired for practical acoustic applications. In our proposed concept, the designed AMEQ can realize different specific amplitude controls with different values (-2dB to -10dB, with a step of -2dB) within ultra-broadband (200Hz-20000Hz). This design goal is realized by using two kinds of acoustic resonant elements, viz., Helmholtz resonators (HRs) [START_REF] Zhu | Multifunctional acoustic metasurface based on an array of Helmholtz resonators[END_REF][START_REF] Zhu | Systematic design of multiplexed-acoustic-metasurface hologram with simultaneous amplitude and phase modulations[END_REF] and Fabry-Perot resonators (FPRs), that is, an array of narrow tubes with different depths [START_REF] Yang | Optimal sound-absorbing structures[END_REF][START_REF] Jiang | Ultra-broadband absorption by acoustic metamaterials[END_REF]. We realize both low-frequency and high-frequency equalizations with these two acoustic elements, respectively, in an integrated acoustic meta-structure. Although these acoustic resonant elements have been demonstrated for sound absorption [START_REF] Yang | Optimal sound-absorbing structures[END_REF][START_REF] Jiang | Ultra-broadband absorption by acoustic metamaterials[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF] before, realizing ultra-broadband filter is challenging due to the very large span of the targeted frequencies (7 octaves). This makes the size of whole structure sub-wavelength for low frequencies, but larger than wavelength for high frequencies. Therefore, the low frequency filters are owing to the coupled resonance of the whole structure, while the high frequency filters are determined by both the resonance of FPRs and geometry acoustics effect. This complex case requires deliberate design (including the geometric feature of the channel, the size and position of each resonator). It is distinct from conventional acoustic absorber designs, which will be demonstrated numerically, analytically, and experimentally in the following sections of this article.

We illustrate in Fig. 1 the classical EQ concept in different domains of applications. Figure 1(a) shows an audio EQ device, which is widely used in sound effects manipulations. Figure 1(b) exhibits the EQ software which can adapt the sound frequency ingredients of an audio on computer. The amplitude response of different frequencies can be freely modulated within specific range. Figure 1(c) shows that underwater communications relay on the frequency equalization process since the complex hydroacoustic environment needs frequency ingredients correction. Our proposed acoustic meta-equalizer is shown in Fig. 1(d). It consists of an acoustic metamaterial with several modulators. The latter are designed for specific amplitude control at different frequencies. By geometrically optimizing the metamaterial frame size and the parameters of HRs and FPRs, we can achieve high precise manipulations of the transmission coefficient for each octave. An integrated physical model is proposed to calculate the transmission for low and high frequencies, respectively. The corresponding analytical results agree well with the simulated and measured ones. Here, in analogy with conventional EQ device, we will show the realization of different functional filters, signal reproduction, and the modulation of specific sound effects by the AMEQ. Our finding integrating multiple functionalities will lead to real-world applications of acoustic metamaterials in room acoustics and architectural acoustics, etc.

II. RESULTS

A. Meta-structure design and analytical model

The conceived AMEQ is a two-port structure as schematically shown in Fig. 2(a). The acoustic wave is incident from one port and transmit to another port as the arrows indicate. The whole structure can be regarded as a "black box" with frequency equalization function. The two-dimensional (2D) structure is designed as a furcate channel with several modulators on both sides of the channel. The main difficulty and challenge for AMEQ design is that the modulators should cover distinctly different resonant frequencies within ultra-broadband. Inspired by previous works, such as acoustic absorber and diffuser designs [START_REF] Yang | Optimal sound-absorbing structures[END_REF][START_REF] Zhu | Ultrathin acoustic metasurface-based Schroeder diffuser[END_REF], we introduce both HRs and FPRs as the low and high frequency modulators, respectively. Furthermore, in order to combine these two elements in an integrated structure, the furcate channel is deliberately designed with optimized geometry.

The modulators at the bottom and the top are used to control low frequencies (three center frequencies with f1=250Hz, f2=500Hz, and f3=1000Hz). The modulators in the middle are exploited to modulate high frequencies (four center frequencies with f4=2000Hz, f5=4000Hz, f6=8000Hz, and f7=16000Hz). different sizes [START_REF] Zhu | Multifunctional acoustic metasurface based on an array of Helmholtz resonators[END_REF][START_REF] Zhu | Systematic design of multiplexed-acoustic-metasurface hologram with simultaneous amplitude and phase modulations[END_REF]. The red structure is corresponding to high frequencies, and the corresponding filter structure is based on FPRs with different lengths [START_REF] Yang | Optimal sound-absorbing structures[END_REF][START_REF] Jiang | Ultra-broadband absorption by acoustic metamaterials[END_REF]. The whole structure is symmetrical in the four quadrants (marked by A, B, C, and D in Fig. 2(c)), so that each resonator has four same duplicates in the panel. Figure 2(e) shows the effective circuit model of the AMEQ, which is electrical element analogy [START_REF] Zhu | Multifunctional acoustic metasurface based on an array of Helmholtz resonators[END_REF][START_REF] Zhu | Systematic design of multiplexed-acoustic-metasurface hologram with simultaneous amplitude and phase modulations[END_REF] corresponding to the analytical derivations we will present in the following section. The arrows indicate the input and output ports.

We achieve the precise modulation of the transmitted amplitude by optimizing the resonant frequencies distributions fnm of the different resonators as shown in Fig. 2(f), in which y-axis is scaled with log function. Red marks in Fig. 2(f) are the resonant frequencies fnm of m th HR unit cell of n th (n=1,2,3) supercell contributing to n th center frequency (fn), while blue marks are the resonant frequencies fnm of m th FPR unit cell of n th (n=4,5,6,7) supercell contributing to fn. For example, for 6 th center frequency (n=6, f=8000Hz), 6 resonators are designed with the resonant frequencies of f61=6400Hz, f62=7000Hz, f63=7640Hz, f64=8420Hz, f65=9040Hz and f66=9720Hz to ensure the broadband feature, as the blue marks in Fig. 2(f) display. The whole structure formed by these resonant frequency distributions can achieve precise modulations of amplitude with desired values.

To theoretically approach the properties of our concept, we investigate the specific resonators design with analytical expressions. For low frequency filters, a HR array is designed.

As shown in Fig. 2(d), the three supercells of the HRs correspond to the three center frequencies (f1-f3), having 3, 4, and 10 HRs, respectively. In principle, the center frequency for each filter band is the average of resonant frequencies of different HRs, viz., fn=`fnm, (n=1,2,3). The resonant frequency of HR is expressed as: [START_REF] Qureshi | Adaptive equalization[END_REF] where Mnm=ρ0lnm/wnm and Cnm=Lnmdnm/ρ0c0 2 are acoustic masses and acoustic capacitances of the HR, respectively. ρ0=1.21kg/m 3 and c0=343m/s are the mass density and the sound speed of air, respectively. Lnm and dnm are the height and the width of the cavity as marked in Fig. 2(d). lnm and wnm are the height and the width of the neck, respectively.

Then, the effective impedance of HRs in effective circuit model is calculated by [START_REF] Zhu | Systematic design of multiplexed-acoustic-metasurface hologram with simultaneous amplitude and phase modulations[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF]:

Znm=i(ωnmMnm-1/ωnmCnm)+Rnm, ( 2
)
where ωnm=2πfnm are the angular frequencies, Rnm is the acoustic resistance induced by effects of

On the other hand, for the high frequency filters (red), an FPR array is designed. Similar to low frequency structure, the center frequency for each high frequency filter band is the average of all resonant frequencies of different FPRs fn=`fnm, (n=4,5,6,7). The length and the width is hnm and wnm, respectively. The filter frequency of FPR is expressed as [START_REF] Yang | Optimal sound-absorbing structures[END_REF]:

(

The high frequency acoustic filter is directly decided by the sound abospotion in the FPR unit cells. On the other hand, the high frequency acoustic filter will be also influenced by multiple sound reflections in the integrated meta-structure. We characterize the high frequency acoustic filter by acoustic resistance Rn=`Rnm (n=4, 5, 6, 7) in the effective circuit. The effective impedance of FPRs is calculated by [START_REF] Yang | Optimal sound-absorbing structures[END_REF]:

0 0 0 1 0 2 2 8 4 2( ) 2 8 1 , nm nm nm nm nm nm nm nm nm nm nm L c J k w R w w w kw µw r r µw r p p p é ù » + + - ê ú ë û -1 1 ( 1,2,3 1 ( ) , ) m j nm Z n Z n = ae ö = = ç ÷ è ø å 0 4 nm nm c f h = 9 (6)
Here, γ is the ratio of the surface area occupied by FPR's cross sections to the effective surface area exposed to sounds. M is total number of the FPR for n th center frenquency. β describes the weak system dissipation for the acoustic metamaterial.

) ) We, then, can modulate the transmission coefficient (acoustic presure amplitude) by controlling the turn on and turn off states of each resonantor using the modulators. Since each resonator has four same duplicates in A, B, C, D quadrants in the panel. The number of the activated (turned on) resonantors can be chosen as δ(n)=0, 1, 2, 3, or 4. Thus, based on the effective acoustic impedance Z(n), the transmission T(n) of the AMEQ for n th center frequencies are related to the number of δ(n), calculated as:

M nm m nm c Z i M n n i r w w w w bw p g é = ù = ê ú -- ë û å FIG. 2.
(7a) (7b) (7c) 0 ( ) , ( ( ) 0) n T n T d = = 0 0 0 ( ) ( ) , ( )+ ( ( ) 1) Z n T n T Z n c n r d é ù = ê ú ë = û 0 0 0 0 0 ( )- ( ) , ( )+ ( ( ) 2) Z n c T n T Z n n c d r r é ù = ú û = ê ë (7d) (7e)
Here, T0 is the transmission of a pure channel without any resonators as shown in Fig. 3(a).

Because the channel has an unusual shape, the transmission curve of T0 is not planar, but related to the working frequency, as shown in Fig. 3(b), with the averaged closed to 0.8. Here, we assume T0=0.8 to simplify the following analytical calculations. From the curve of T0, one can see the transmission above 1890Hz is oscillating because it is the approximate cutoff frequency of the channel. Therefore, the low frequency filters from 180Hz to 1890Hz are owing to the resonance of the whole structure. The high frequency filters from 1890Hz to 22000Hz are determined by both the resonance of FPR element as well as geometrical acoustics effects. This is reflected in the parameter γ in Eq. [START_REF] Radlovic | Equalization in an acoustic reverberant environment: Robustness results[END_REF], that is, an effective filling rate related to the geometric feature of the channel and the distribution of FPRs. In this case, the filling rate γ is smaller than the previous case for a compact structure [START_REF] Yang | Optimal sound-absorbing structures[END_REF].

B. Amplitude Control

Owing to the precise modulations by the optimized resonator arrays, the transmitted amplitude is controlled with Levels 0-4 corresponding to amplitudes A=0.8, 0.63, 0.5, 0.4, 0.3 (viz., -2dB, -4dB, -6dB, -8dB, -10dB, respectively). Figure 3(c) shows the schematic diagram for Levels 0-4. In this diagram, the 28 circles represent 28 modulators for different fn (each frequency has 4 tunable modulators). The red circles represent turn off states, for which we fix a red slice covering the opening of the resonator to prevent the effect of the resonator. The gray

0 0 0 2 0 0 ( )( ( )-) ( ) , ( ( ) ( ( ) + ) 3) Z n Z n c T n T Z n c n d r r é ù = ê ú ë û = 2 0 0 0 0 0 ( )- ( ) , ( (
( ) 4) )+ Z n c T n T Z n n c r r d é ù = ê ë = ú û
circles represent turn on states, indicating no slice is fixed. The photographs of the samples for Levels 0-4 are also shown in Fig. 3(c). The sample is tunable by introducing the red slices as modulators [START_REF] Fan | Tunable broadband reflective acoustic metasurface[END_REF]. It has a size of 56cm×56cm×4cm (Horizontal size 56cm is about 1/4 wavelength of the lowest working frequency). The samples and the slices are made by 3D printer with polylactic acid (PLA) having the mass density ρ=1250kg/m 3 , Young's modulus E=3.2×10 9 Pa, and Poisson's ratio ν=0. [START_REF] Lee | Ultrasparse acoustic absorbers enabling fluid flow and visible-light controls[END_REF].

The simulations were performed using commercial finite element analysis software COMSOL Multiphysics 5.4a with the "Acoustic-Thermoviscous Acoustic Interaction, Frequency Domain" and "Thermoviscous Acoustic-Solid Interaction, Frequency Domain" modules. The considered mass density and the sound speed of background medium air are ρ0=1.21kg/m 3 and c0=343m/s, respectively.

The simulated, analytical, and measured filter curves for 5 samples are shown in Fig. 3(c).

The transmission for one center frequency, such as 8000Hz, is obtained by calculating the averaged acoustic pressure amplitude within one octave from 5657Hz to 11314Hz (with a step of 1/24 octave, viz., 25 frequencies). For the experiment, the measurements are carried out in 2D parallel waveguide experimental system [START_REF] Fan | Tunable broadband reflective acoustic metasurface[END_REF]. The distance between two parallel waveguides is 4cm. The size of the parallel waveguides is 0.6m×0.8m. Sound-absorbing cotton is placed at the boundaries of the parallel waveguides to prevent undesired reflections. A 10cm-diameter loudspeaker (matches with the input side size 9.6cm) is fixed at the input boundary of the sample.

The acoustic wave is generated continuously and the static acoustic field is measured.

The output cross-section of the sample is scanned to obtain the output average acoustic wave amplitude, with 1/4-inch-diameter Brüel&Kjaer type-4961 microphone and Brüel&Kjaer PULSE Type 3160. Brüel&Kjaer software is used for measurement, in which "Constant Percentage Bandwidth Analyzer" (CPB) is chosen that the incident frequency is sweeping from 180Hz to 22000Hz. One can conveniently measure the averaged acoustic pressure amplitude in real-time for every octave from f1=250Hz to f7=16000Hz. The results in Fig. 3(c) show that, by controlling the number of the modulators (0, 7, 14, 21, and 28, respectively), we can achieve predesigned amplitude response (0.3-0.8) at the output side.

The analytical results agree with the simulated and experimental ones. It is noted that the simulated amplitudes for 1000Hz and 2000Hz are slightly higher than the predesigned ones. This is an inevitable error, due to the fact that Level 0 is corresponding to the pure channel case, whose transmission T0 is not exactly planar as shown in Fig. 3

C. Functional filter

Different functional filters are important elements in signal processing [START_REF] Qureshi | Adaptive equalization[END_REF], such as the classical ones, i.e., high-pass, low-pass, band-pass, and stop-band. In the following, we show the realizations of the above four functional filters. The schematic diagram and the samples are shown in Fig. 4(a) and 4(d). The distributions of the modulators are optimized to achieve the desired filter curve. The high-pass filter is designed as high-transmission for 4 octaves with center frequencies from 2000Hz to 16000Hz, while low-transmission is designed for 3 octaves with center frequencies from 250Hz to 1000Hz. In Fig. 4(e), the numerical and experimental results show that the amplitude response for 200Hz-1400Hz is about A=0.3, and for 1400Hz-20000Hz is about A=0.6. Similar designs are applied to other filters in Fig. 4

D. Signal reproduction

Signal reproduction is an important concept proposed in room acoustics [START_REF] Spors | Active listening room compensation for massive multichannel sound reproduction systems using wave-domain adaptive filtering[END_REF], that uses an equalization process to reproduce the input signal with random or unbalanced frequency ingredients, and reconstructs it as a balanced one with a planar frequency response curve.

E. Sound effect control

Sound effect control [START_REF] Howard | Acoustics and Psychoacoustics[END_REF] is another fascinating application of EQ device which requires a fine modulation of each frequency ingredient. We predesign four EQ curves (denoted as EQ1-EQ4) whose sound effects are famous in music productions, viz., EQ1 (Hip-hop), EQ2 (Popular), EQ3 (R&B), and EQ4 (Dance), respectively. The AMEQ structures as well as the simulated, analytical, and measured data are illustrated in Figs. 6(a-h). By the optimizing the modulators distribution, we obtain the desired EQ curves at the output side of the AMEQ. The results show that the low frequency control is more precise than the high-frequency one (such as the results at 16kHz). This is due to the fact that high frequency filter is more sensitive to the plane wave incident waveform (direction). The ability of the real sound effect control by AMEQ is significant for the applications in architectural acoustics, when the high-quality hearing effect is desired.

III. DISCUSSION AND CONCLUSION

We have proposed the unique concept of acoustic meta-equalizer capable of providing acoustic wave equalization with a passive acoustic metamaterial. The equalization requires tunable and ultra-broadband frequency filter. By using two kinds of acoustic resonant unit cells (Helmholtz resonators and Fabry-Perot resonators), we have realized ultra-broadband filters within a fully integrated structure. The optimization of the resonators design ensures highly precise amplitude modulations. The obtained bandwidth can reach 7 octaves, viz., 200Hz-20000Hz, while the amplitude can be modulated with 5 levels (-2dB to -10dB with the step of -2dB) realized by tuning 28 modulators. It's worth mentioning that the tunability of the current system is manual via designed slices. For further improving the tunability of the system, a possible method is to introduce origami structure [START_REF] Zhu | Reconfigurable origami-inspired metamaterials for controllable sound manipulation[END_REF] to switch the "turn on" and "turn off" states of the system by using folding slices. The bandwidth and the degree of freedom of amplitude outperform conventional acoustic filters. In analogy with the equalization applications in communication systems, room acoustics, and audio device, we have numerically and experimentally demonstrated and provided the proof-of-concept of functional filters, signal reproduction, and sound effect control, respectively. The proposed concept of acoustic meta-equalizer aims to connect acoustic metamaterials to the practical acoustic applications.

Although the precision of the proposed meta-equalizer is not as high as the matured electronic counterpart, the proposed meta-equalizer may have unique applications in architectural acoustics.

For example, a passive structure with ultra-broadband feature is a good candidate for noise control with controllable frequency spectrum. The concept of equalizer extended from electronic device to acoustic meta-device will lead to more possible applications, such as noise control with ventilation property [START_REF] Lee | Ultrasparse acoustic absorbers enabling fluid flow and visible-light controls[END_REF][START_REF] Sun | Broadband acoustic ventilation barriers[END_REF], beyond the applications of conventional electronic equalizer.

FIG. 1 .

 1 FIG. 1. The concept of equalizer. (EQ) with different potential applications. (a) EQ device in

  Figure 2(b) illustrates the three-dimensional (3D) structure of the AMEQ. The filters for different center frequencies (fn) are marked in the figure.

  Figure 2(c) shows the air channel which is the complementary part of the solid structure in Fig.2(b). Sound is propagating in the air channel.

Figure 2 (

 2 d) exhibits the whole 2D meta-structure of the AMEQ. Dashed box is the enlarged view of the solid box, showing the structural schematic diagram of each unit cell, with different geometrical parameters marked in the figure. The blue structure is corresponding to low frequencies. The corresponding filter structure is based on HRs with

  radiation and thermal-viscosity. It can be formulated as: (3) where µ=1.983×10 -5 Pas is the dynamic viscosity of air, and J1 is the Bessel function of the first kind. Combining the effective circuit model, we can explain the low frequency acoustic filter as the result of the sound reflection and the thermal-viscous loss (corresponding to the energy dissipation at Rn=`Rnm) at the resonances of the HRs (corresponding to the resonances of the circuit branch).After obtaining the acoustic impedance Znm for each individual HRs, we can get the total effective impedance Z(n) of the AMEQ for 3 different center frequencies.

  FIG. 2. The schematic diagram of the acoustic meta-equalizer (AMEQ). (a) The AMEQ device

  FIG. 3. Amplitude Control. (a) The pure channel without any resonators. (b) The transmission of

  FIG. 4. Functional filter. (a-d) The schematic diagram of modulators distribution and the

Figure

FIG. 5 .

 5 FIG. 5. Signal reproduction. (a) The input signals with unbalanced frequency amplitude

FIG. 6 .

 6 FIG. 6. Sound effect control. (a-d) The schematic diagram of modulators distribution and the

ACKNOWLEDGMENTS

This work is supported by the Air Force Office of Scientific Research under award number FA9550-18-1-7021, and by la Région Grand Est and Institut Carnot ICEEL.