Ulisse Ferrari
email: ulisse.ferrari@gmail.com

Stéphane Deny

Olivier Marre

Thierry Mora

A simple model for low variability in neural spike trains

Keywords: Sensory systems, Refractoriness, Fano factor, Under-dispersion, Spike counter, Retina

Neural noise sets a limit to information transmission in sensory systems. In several areas, the spiking response (to a repeated stimulus) has shown a higher degree of regularity than predicted by a Poisson process. However, a simple model to explain this low variability is still lacking. Here we introduce a new model, with a correction to Poisson statistics, which can accurately predict the regularity of neural spike trains in response to a repeated stimulus. The model has only two parameters, but can reproduce the observed variability in retinal recordings in various conditions. We show analytically why this approximation can work. In a model of the spike emitting process where a refractory period is assumed, we derive that our simple correction can well approximate the spike train statistics over a broad range of firing rates. Our model can be easily plugged to stimulus processing models, like Linear-nonlinear model or its generalizations, to replace the Poisson spike train hypothesis that is commonly assumed.

It estimates the amount of information transmitted much more accurately than Poisson models in retinal recordings. Thanks to its simplicity this model has the potential to explain low variability in other areas.

Introduction

Neural variability imposes constraints on the way neurons transmit and process information [START_REF] Movshon | Reliability of neuronal responses[END_REF] and has been extensively studied in the mammalian visual system [START_REF] Barlow | Three factors limiting the reliable detection of light by retinal ganglion cells of the cat[END_REF]; [START_REF] Heggelund | Response variability and orientation discrimination of single cells in striate cortex of cat[END_REF]; Tolhurst et al. (1981); [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF] and beyond de Ruyter van [START_REF] De Ruyter Van Steveninck | Reproducibility and variability in neural spike trains[END_REF]; [START_REF] Baddeley | Responses of neurons in primary and inferior temporal visual cortices to natural scenes[END_REF]; [START_REF] Deweese | Binary coding in auditory cortex[END_REF]; [START_REF] Churchland | Stimulus onset quenches neural variability: a widespread cortical phenomenon[END_REF]. The Fano Factor, i.e. the ratio of the spike count variance to the mean, is often used to characterize such variability. Many studies have reported various values for Fano Factor depending on brain area [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF] or experimental condition [START_REF] Churchland | Stimulus onset quenches neural variability: a widespread cortical phenomenon[END_REF]. A Fano Factor lower than one suggests that neurons spike regularly, a condition sometimes termed as "under-dispersion" [START_REF] Gur | Response variability of neurons in primary visual cortex (v1) of alert monkeys[END_REF]; [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF]; [START_REF] Barberini | A comparison of spiking statistics in motion sensing neurons of flies and monkeys[END_REF]; [START_REF] Deweese | Binary coding in auditory cortex[END_REF]; [START_REF] Maimon | Beyond poisson: increased spike-time regularity across primate parietal cortex[END_REF]. A Fano Factor above one is termed over-dispersion [START_REF] Baddeley | Responses of neurons in primary and inferior temporal visual cortices to natural scenes[END_REF]; [START_REF] Churchland | Stimulus onset quenches neural variability: a widespread cortical phenomenon[END_REF].

A common model for neural spike count is to assume Poisson statistics, where the variance in the number of spikes emitted is equal to its mean and the Fano Factor is thus one. Like the classical "LNP" (Linear Nonlinear Poisson) [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF], many models describing how stimuli are processed by individual neurons rely on this assumption. New models are therefore needed to account for the deviations from Poisson statistics observed in the data. Several models have been proposed to account for overdispersion in the spike count distribution [START_REF] Goris | Partitioning neuronal variability[END_REF]; [START_REF] Scott | Fully bayesian inference for neural models with negative-binomial spiking[END_REF]; [START_REF] Charles | Dethroning the fano factor: a flexible, model-based approach to partitioning neural variability[END_REF]. For under-dispersed spike count distribution, a few models have been proposed, but they come with specific constraints. One possibility is to use a spike history filter that allows past spikes to inhibit the probability of emitting a new spike at present [START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF]. However, this approach requires defining the probability of spiking over very small time bins (ex: 1ms), and consequently needs several parameters or strong regularization to describe the spike history filter. Other models have been proposed, but some have many parameters (for example, ∼ 7 for each cell in [START_REF] Gao | High-dimensional neural spike train analysis with generalized count linear dynamical systems[END_REF]), and may therefore require very large datasets to learn the model, or make specific assumptions that may not always be verified in the data [START_REF] Sellers | The com-poisson model for count data: a survey of methods and applications[END_REF]; [START_REF] Stevenson | Flexible models for spike count data with both over-and underdispersion[END_REF]. Overall, it is unclear if a general yet simple model can account for the spike count distribution found in the data [START_REF] Charles | Dethroning the fano factor: a flexible, model-based approach to partitioning neural variability[END_REF].

Here we present a simple model that can account for under-dispersion in the spike count variability, with only two parameters. Our starting assumption is that the deviation from Poisson statistics comes from the refractory period, i.e. the fact that there is a minimal time interval between two spikes. We start from the analytic form of the spike count probability taken by a spike generation process composed of an absolute refractory period followed by a Poisson process, and show that it can be approximated by a model with only one parameter to fit to the data. We further simplify the model to make it amenable to log-likelihood maximization. We then relax our assumption by allowing for a relative refractory period [START_REF] Berry | Refractoriness and neural precision[END_REF]; Berry and Meister (1998) and we derive a more flexible model with two parameters that can accurately predict neural variability. The simple form of this model makes it possible to plug it on classical stimulus processing models (e.g. Linear-Nonlinear (LN) model [START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF] or more complex cascade models [START_REF] Deny | Multiplexed computations in retinal ganglion cells of a single type[END_REF]), that can be fitted with log-likelihood maximization. We test our model on retinal data and find that it outperforms the Poisson model, but also other statistical models proposed in the literature to account for under-dispersion. Either our model performs better at describing the data, or other models need more parameters for an equally accurate description. When combined with classical stimulus processing models, our model is able to predict the variance of the spike count over time, as well as the amount of information conveyed by single neurons, much better than classical models relying on Poisson processes. We thus propose a simple model for neural variability, with only two parameters, that can be used in combination with any stimulus processing model to account for sub-Poisson neural variability.

Sub-Poisson behavior of retinal ganglion cells

We used previously published data [START_REF] Deny | Multiplexed computations in retinal ganglion cells of a single type[END_REF] where we recorded ganglion cells from the rat retina using multi-electrode arrays [START_REF] Marre | Recording of a large and complete population in the retina[END_REF]Yger et al. (2016).

Two types of ganglion cells, ON and OFF, with receptive fields tiling the visual space, were isolated in these experiments (n=25 and n=19 cells, respectively). Cells were then stimulated with the video of two parallel horizontal bars performing a Brownian motion along the vertical direction [START_REF] Deny | Multiplexed computations in retinal ganglion cells of a single type[END_REF]. Some stimuli sequences were repeated and triggered reliable responses from ganglion cells (Fig. 1A&B top panels).

We then binned the RGC responses with 16.7ms time windows ∆T locked to the frame update of the video stimulation (60Hz). For each cell i in each time bin t and for each stimulus repetition rep we associated an integer spike count n rep i (t) equal to the number of emitted spikes in that time-window. In order to analyze the cell reliability as a function of the firing rate, for each cell and time bin we computed the mean (λ i (t))

and variance (V i (t)) of this spike count across repetitions of the same stimulus. In Fig. 1A&B (bottom panels) we show the firing rate of two example ON and OFF cells, together with its fluctuation across the repetitions (colored areas). For comparison we also show the fluctuations that can be expected from a Poisson process; that is, with a variance equal to its mean. During transients of high activity, cells are more reliable than what can be expected.

Figs. 1C&D

show that both ON and OFF cells have a sub-Poisson variability. Variances of spike counts are much lower than their means and their ratio (the Fano factor, equal to 1 for Poisson distributions) is not only smaller than one, but decreases with the mean. Note that here and for the rest of the paper, we pool together all the cells of the same type when estimating the parameters of the spike generators. We found that spike statistics were remarkably homogeneous across cells of the same type, and this gives us more data to fit different models. However, the same process could have been applied on a single cell, with enough repetitions of the same stimulus.

These results are consistent with previous findings [START_REF] Berry | Refractoriness and neural precision[END_REF]; Berry and Meister (1998); [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF]: ganglion cells emit spikes more reliably than a Poisson process. Consequently, a model for predicting RGC response that accounts for noise with a Poisson generator, for example the LNP model [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF], largely under-estimates this spiking precision. We then search for a simple model to account for this relation between signal and noise.

3 Spike count statistics for a refractory neuron Absolute or relative refractory periods are known to impose some regularity on the sequence of emitted spikes [START_REF] Berry | Refractoriness and neural precision[END_REF]; [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF] and thus to decrease the neural variability. We aim at looking how refractoriness impacts the spike count distribution of the spike train. We consider a model of refractory neuron where the instantaneous rate is inhibited by an absolute refractory period of duration τ . Such a neuron has the following Inter-Spike interval (ISI) distribution: -τ ) .

ρ Θ (t) = Θ(t -τ )re -r (t
(

) 1
where Θ(u) the Heaviside unit step (Θ(u) = 1 for u > 0 and zero otherwise) and r is the firing rate in absence of refractoriness. From the ISI distribution (1) we can estimate the probability distribution of the spike count n, i.e. the probability distribution of the number of spikes emitted in a time bin ∆t. For the particular case of τ = 0, then λ ≡ n = r∆t and n follows a Poisson distribution:

P Pois n λ = λ n n! e -λ . (2)
For τ > 0, n follows a sub-Poisson distribution P Θ ( n | r, ∆t, τ ), which admits an analytic (albeit complex) expression [START_REF] Müller | Some formulae for a dead-time-distorted poisson process: To andré allisy on the completion of his first half century[END_REF] (see Eq. ( 13) in the Methods and the derivation in the Supplementary Information).

To test if this expression accounts well for the data, we need to adjust the value of the refractory period τ for both cell types. Because the complex expression of

P Θ ( n |

ON-cells

Figure 2: A simple neuron model with an absolute refractory period accounts for the observed spike count statistics. A) Relation between mean and variance predicted by the model (pink) and measured for the OFF population (green, same data as Fig. 1).

B) Same as A, but for the ON population.

4 Simple models for the spike count statistics

We have shown how a rather simple model of refractory neurons accurately reproduces the mean-variance relation of the recorded cells. However, the model has a complicated analytic form, which is not easily amenable to a likelihood-maximization approach, and it cannot be considered as an alternative to models based on Poisson generators such as the LNP model. To overcome this limitation, we propose a further simplification of the refractory neuron that is much more efficient and tractable.

For τ = 0, the refractory distribution P Θ ( n | r, ∆t, τ ) reduces to the Poisson distribution (Eq. 2). The small values of the refractory periods relative to the bin size suggests that an expansion at small τ could capture most of the model's behavior. We expand P Θ ( n | r, ∆t, τ ) around the Poisson distribution to the second order in the small parameter f = τ /∆T (see supp. 8 for details):

P (2nd) Θ ( n | λ, f ) = exp θ λ n -(f -f 2 )n 2 - f 2 2 n 3 -log n! -log Z λ (3) θ λ chosen such that n P (2nd) Θ ( n | λ,f ) = λ . (4)
where Z λ is a normalization constant. The "Second-Order" model (3) has only one free parameter, f or equivalently τ . θ λ is a function of λ and has to be estimated numerically in order to reproduce the mean spike count λ once the refractory period τ has been fixed. For instance, for τ = 0, θ λ = log λ and the Second-Order reduces to the Poisson distribution. Importantly, the coefficients of the n 2 and n 3 terms do not depend on the firing rate r, insofar as the refractory period does not either. The exponential form of the model makes it easy to calculate its derivatives and to use with maximum-likelihood methods.

We infer τ through a log-likelihood maximization (see Methods) using Eq. 3, for OFF and ON populations separately. From these values we can estimate the refractory period for the Second-Order model: τ = 10.8ms for OFF cells and τ = 3.0ms for ON cells. In Fig. 3 we compare the model predictions for spike count variance V with the empirical values. Unlike in the previous section, model parameters were not optimized to best fit this curve, since it was fitted using log-likelihood maximization. Yet the Second-Order model shows a high performance for the OFF population, even at large firing rates. Despite the approximation, the model is therefore still able to describe accurately the mean-variance relation, and can now be fitted using maximum likelihood approaches. For ON cells, the model outperforms the Poisson model but the meanvariance relation is not perfectly predicted. We will then explore why this could be the case and improve our model to get a more flexible one, that can be suited for a broader range of experimental cases.

A simple Effective model to describe spike trains statistics

The previous results hold for a very specific assumption about how refractoriness constrains firing, with just a single fitting parameter -the refractory period. We thus won- We considered a general model in which the instantaneous spike rate is now inhibited by a time-dependent factor, α(r, u)r, where r is the spike rate in absence of refractoriness, and u is the time following the previous spike, with α(r, 0) = 0 and α(r, u ≥ τ ) = 1. For instance, a pure refractory period is defined by α(r, u) = 0 for u < τ and α(r, u) = 1 for u > τ . A soft refractory period, where the spike rate recovers linearly in time, would be given by α(r, u) = u/τ for u < τ and α(u, t) = 1 for u > τ . Calculating the spike count distribution under this assumption is intractable analytically, but a Second-Order expansion such as the one performed in the previous section can still be performed, and yields the following expression at Second-Order (see Supp 8):

P (2nd) α ( n | λ, r, f, ∆t ) = exp θ λ n -γ α n 2 -δ α n 3 -log n! -log Z λ (5) γ α = f -f 2 + (2 + r∆t) f 2 2 -g δ α = f 2 -g = f 2 2 + f 2 2 -g θ λ chosen such that n P (2nd) α ( n | λ,r,f,∆t ) = λ . ( 6
)
where

f ≡ ∞ 0 du (1-α(r, u))/∆t and g ≡ ∞ 0 du u (1-α(r, u))/∆t 2 .
The special case of a pure refractory period, α(r, u) = 0 for 0 ≤ u ≤ τ , gives back Eq. 3 and the previous definition of f . In this case, f 2 /2 -g vanishes. This quantity can thus be considered as an estimate of the deviation from the absolute refractoriness. Again, θ λ should be adjusted to match the average number spikes in each cell and time bin, λ i (t).

This analytic development shows that the coefficients of n 2 and n 3 , γ α and δ α , can have very different forms depending on the exact form of the refractoriness. We thus decided to relax the assumption of a strict dependence between γ α and δ α . We tested if a model with γ α and δ α that does not depend on r shows a good agreement with the data.

In the following we thus treat the coefficients of n 2 and n 3 in Eq. ( 5) as two parameters that are independent of each other, and also independent of the firing rate.

The resulting Effective model is defined for any mean spike count λ as:

P Eff ( n | λ, γ, δ ) = exp θ λ n -γn 2 -δn 3 -log n! -log Z λ (7) θ λ chosen such that n P Eff ( n | λ,γ,δ ) = λ . ( 8
)
where as before θ λ is not a free parameter, but a uniquely defined function of λ with parameters γ and δ. While λ is set to its time-dependent value λ t , γ and δ are assumed 2012) (Eq. 26) has been proposed to account for both under-and over-dispersed mean-variance relation [START_REF] Stevenson | Flexible models for spike count data with both over-and underdispersion[END_REF]. COMP is a one-parameter extension of Poisson, which differs from our Effective model, but it is still a particular case of Gen.Count (see Methods for details). For the sake of completeness we also compare the Refractory model, P Θ of Eq. ( 13) and its second order expansion, the Second-Order model, P 2nd Θ of Eq.(3). We compare the performance of different models as improvement over the Poisson log-likelihood (fig. 4). For the OFF population (fig. 4A), all models outperform Poisson and have similar performance. This is probably because OFF cells show rather small firing rates and rarely emit more than two spikes in the same time-bin. The addition of one parameter (with respect to Poisson) is thus enough for accurately modeling the spike count statistics. However, for the ON population (fig. 4B), while all the considered models outperform Poisson, Effective and Gen.Count show the largest improvement.

Remarkably, the one-parameter Refractory and Second-Order models show very high performance as well, despite the first being learned by fitting the mean-variance relation of the spike count rather than by maximizing the likelihood. Also in this case, the larger

M I(n, λ) = H[ P (n) ] -H[ P (n|λ) ] λ ( 9
)
where H is the entropy function and P (n) is the distribution of the spike count without conditioning on the mean over repetitions. . . . λ is the average over the observed mean spike count and P (n|λ) is the distribution of n at fixed λ, either empirically estimated from the repeated data, or estimated with the model. To avoid under-sampling issues typical of time bins with low activity and due to the finite number of stimulus repetitions, we restrict the averages to all cells and time-bins with λ i (t) > 0.1. In Fig. 6, for both OFF and ON population, Poisson largely under-estimates mutual information, whereas the Effective model predicts well the value of mutual information.

Our model can thus be used to correctly estimate the mutual information thanks to its accurate prediction of the mean-variance relation.

Note that ( 9) is a mutual information for a single time bin, and cannot be used to calculate information rates over spike trains because of its non-additivity across time due to stimulus temporal correlations. See Fig. 8 for further details and explanations.

Improving stimulus processing models

The Effective model, Eq. ( 7), describes efficiently the relation between mean and variance and predicts well the amount of information transmitted by a neuron. Here we show how it can easily be plugged to any stimulus processing model like the Linear-Nonlinear (LN) model, instead of a Poisson process. To estimate the parameters of LN models, a classical approach is to assume a Poisson process for spike generation, and then to maximize the likelihood of the spiking data with respect to the parameters of the model. The major advantage of this method is that the gradient of the likelihood has a very simple form that allows for iterative log-likelihood maximization. Here we show that using the Effective model also leads to a tractable form of the log-likelihood gradient, which can similarly be used for iterative optimization, but with the added advantage that the mean-variance relationship is accurately reproduced.

In general, a stimulus encoding model is defined by a series of computationsparametrized by parameters ψ -that takes the past stimulus S t as input and provides a prediction λψ (S t ) for the spike count mean as a function of time t. Only at the last stage a stochastic spike counter P is introduced to predict the number of spikes n(t) emitted in the time-bin t:

S t → λψ (S t ) → n(t) ∼ P n(t) λψ (S t ) (10
∇ ψ (ψ) = 1 T t d log P Pois n(t) λ dλ λ= λψ (St) ∇ ψ λψ (S t ) = 1 T t n(t) -λψ (S t ) λψ (S t ) ∇ ψ λψ (S t ) , (11)
Once ∇ ψ λψ (S t ) is evaluated, (11) allows for iterative log-likelihood maximization.

The particular structure of the Effective model, Eq. ( 7), allows for an easy estimation of the log-likelihood gradient (see Methods for the calculation details):

∇ ψ (ψ) = 1 T t n(t) -λψ (S t ) V Eff λψ (St) ∇ ψ λψ (S t ) (12)
With respect to the Poisson case, the inference of the parameters ψ requires only the variance V Eff λ as a function of λ, which depends only on the noise model and can be easily estimated from it before running the inference (see Methods for more details). Poisson by the Effective model leads to a significant performance improvement when trying to predict the variance (fig. 7D). For OFF cells, the Second-Order model, Eq. 3, performs as well as the Effective model. The Effective model also predicts well variance for ON cells (fig. 7E). Our Effective model can therefore be plugged in encoding models to describe accurately the variance of the spike count over time.

Discussion

We areas [START_REF] Gur | Response variability of neurons in primary visual cortex (v1) of alert monkeys[END_REF]; [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF]; [START_REF] Barberini | A comparison of spiking statistics in motion sensing neurons of flies and monkeys[END_REF]; [START_REF] Deweese | Binary coding in auditory cortex[END_REF]; [START_REF] Maimon | Beyond poisson: increased spike-time regularity across primate parietal cortex[END_REF], even if their neural variability is not solely determined by the refractory period in the spike generation of the neuron. Previous works [START_REF] Movshon | Reliability of neuronal responses[END_REF]; [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF] suggested that the refractory period present in the phenomenological model may reflect refractoriness at any stage in the circuit and may not directly correspond to the refractoriness of the recorded cell.

We have compared our model with others already proposed in the literature: the Conwey-Maxwell-Poisson (COMP) [START_REF] Stevenson | Flexible models for spike count data with both over-and underdispersion[END_REF] and the Generalized Count (Gen.Count) [START_REF] Gao | High-dimensional neural spike train analysis with generalized count linear dynamical systems[END_REF]. We found the COMP model was quite inefficient at fitting our data.

The Gen.Count model can be considered as encompassing a much larger class of possible models, and includes ours as a special case. However, this comes at the cost of having many more parameters to fit, which can lead to overfitting in some cases, as we have shown in the results. Moreover, in our data, these additional parameters did not allow improving the performance of the model. The model we propose has only two parameters, which makes it easy to fit and usable in many cases, even when the amount of data is limited. It will be interesting to compare our effective models to the Gen.Count model in other sensory structures.

The relation between mean and variance of the spike count strongly depends on the bin size chosen to bin the cell response [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF]. We illustrate this fact on our data by plotting the Fano factor as a function of window size in Fig. 8A. One can distinguish three regimes as a function of the window size relative to the other time scales of the system (Fig. 8B) -the inverse of the spike rate r -1 , the time scale of spike rate variations τ λ , and the refractory period τ refr , the last two of which can be directly read off from the autocorrelation spiking activity (Fig. 8B). If the window size is of the same order compared to the inverse of the spike rate, then the mean spike count λ is small, and the spiking process can be described by a Bernouilli binary variable, with variance V = λ(1-λ) and Fano factor 1-λ ≈ 1. Further, if the window size is smaller or even of the order of the refractory period, then spiking events in subsequent windows cannot be treated independently, and the modeling strategy of focusing on single time windows must be revised. In this case a common solution is to introduce a spike history filter [START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF] that models the spike probability with a dependence on the past activity. At the other extreme, if the window size is large compared to the typical time scale over which the spike rate varies, then each window will contain a mixture of spike rates, leading to a Fano factor larger than one [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF]. We have chosen a bin size between these two extremes, such that the spike count is not a binary variable (∆t ∼ r -1 ), but the firing rate stays roughly constant within a single time bin (∆t τ λ ), and the dependence between spikes in consecutive time bins is relatively weak (∆t τ refr ). The advantage of this choice of bin size is the small number of parameters (only two) needed to describe the spike regularity. We also note that the bin size chosen also corresponds to the timescale of the retinal code Berry and Meister (1998). To test the robustness of our modelling approach to the choice of window size, we have fitted an Effective model (Eq. 7) to the distribution of spike counts taken in a window of 33 ms instead of 16 ms, and still found excellent agreement (Fig. 8D). The interplay between the different time scales discussed here will vary across brain areas and conditions, and the right window size should be picked as a function of these considerations for our approach to be applicable.

Here we have focused on the spike count statistics of single cells. The activity of other cells in the network is also expected to affect spike count variability through noise correlations. Future work should encompass both sources of non-Poissonian variability through self and cross correlations in a single framework [START_REF] Ferrari | Separating intrinsic interactions from extrinsic correlations in a network of sensory neurons[END_REF].

In this study we have pooled cells of the same type (ON and OFF), because they are known to have equivalent response properties [START_REF] Devries | Mosaic arrangement of ganglion cell receptive fields in rabbit retina[END_REF]. However, we checked that fitting the Effective model on single cells rather than on pooled data also gave excellent performance, with even an improvement for about half of the cells (14/25 for OFF cells, 7/19 for ON cells), with a ≈ 11% and ≈ 16% increase in log-likelihood gain relative to the Poisson model for Off and On cells respectively. This demonstrates the potential applicability of our method to brain areas where cells cannot be pooled by types.

Methods

Equilibrium Poisson process with absolute refractory period. A Poisson process with firing rate Θ(u -τ )r, where Θ is the Heaviside step function of the time from the previous spike, has inter-spike interval distribution given by Eq. ( 1). For such a process it is possible to compute the probability distribution of the number a spikes emitted 20 in a time window ∆t, whose starting point is chosen as random [START_REF] Müller | Dead-time problems[END_REF][START_REF] Müller | Some formulae for a dead-time-distorted poisson process: To andré allisy on the completion of his first half century[END_REF], P Θ n r, τ, ∆t = P Θ n ν, f , where f = τ /∆t and ν = r∆t (see SI for a derivation):

P Θ n ν, f = 1 1 + νf Φ(n) +Θ(n Max -2 -n) n j=0 (n + 1 -j) ν j (1 -(n + 1)f ) j e -ν(1-(n+1)f ) j! -2Θ(n Max -n -1) n-1 j=0 (n -j) ν j (1 -nf ) j e -ν(1-nf ) j! + Θ(n Max -n) n-2 j=0 (n -1 -j) ν j (1 -(n -1)f ) j e -ν(1-(n-1)f ) j! (13) Φ(n) ≡          0 , n ≤ n Max -2 n Max (1 + νf ) -ν , n = n Max -1 ν -(n Max -1)(1 + νf ) , n = n Max (14)
where n Max is the smallest integer larger than ∆t/τ and we used the convention Θ(0) = 1. The distribution ( 13) has an expected value given by:

E Θ (n|ν, f ) = ν 1 + νf (15)
and exact variance [START_REF] Müller | Some formulae for a dead-time-distorted poisson process: To andré allisy on the completion of his first half century[END_REF]:

V Θ (n|ν, f ) = 2 n Max -1 n=0 ν(1 -nf ) -n + n-1 j=0 (n -j) ν j (1-nf ) j e -ν(1-nf ) j! -ν -ν 2 1+νf 1 + νf (16)
Inference of the absolute refractory period model. To infer the value of f for the absolute refractory model, see Eq. ( 13), we perform a mean square error (MSE) minimization of the mean-variance relation:

M SE ≡ i,t V i (t) -Var Θ (n|ν i (t), f ) 2 (17) ν i (t) = λ i (t) 1 -λ i (t)f (18)
where λ i (t) and V i (t) are the empirical mean and variance of n for the cell i in the timebin t. We used the exact expression ( 16) for Var Θ (n|ν i (t), f ) although an approximated results can be obtained by the simpler asymptotic expression [START_REF] Müller | Some formulae for a dead-time-distorted poisson process: To andré allisy on the completion of his first half century[END_REF]:

V Θ (n|ν, f ) ≈ ν (1 + νf ) 3 . ( 19
)
ized Count (Gen.Count) distribution, that in our notation and framework reads:

P Gen.Count ( n | λ ) = exp {θ λ [G] n + G[n] -log n! -log Z λ [G]} (25)
where G[n] is a generic real function defined on the non-negative integers n ∈ [0, ∞].

To better characterize its λ dependence, we have rewritten Gen. Equipping stimulus processing model with Second-Order or Effective noise term. In this appendix we detail the calculation for computing the log-likelihood gradient of Eq. ( 12). For a general stimulus processing model, see Eq. ( 10), equipped with noise term P the log-likelihood reads:

P COMP ( n | λ ) = exp {θ λ [η] n -η log n! -log Z λ [η]} (
(ψ) = 1 T t log P n(t) λψ (S t ) (27)
where the summation runs over the duration of the training set. Consequently the log-

likelihood gradient ∇ ψ (ψ) reads ∇ ψ (ψ) = 1 T t d log P n(t) λ dλ λ= λψ (St) ∇ ψ λψ (S t ) . (28)
To estimate ∇ ψ (ψ) we thus need to compute the derivative of the log-probability with respect to the mean spike count. If P is the Poisson distribution P Pois , then:

d dλ log P Pois n(t) λ = d dλ n(t) log λ -λ -log n(t)! = n(t) -λ λ . ( 29
)
If P is instead the Effective model:

d dλ log P Eff n(t) λ = d dθ λ log P Eff n(t) λ d dλ θ λ (30) = (n(t) -λ) d dθ λ λ -1 = n(t) -λ V Eff λ . (31)
The very same expression, with replaced V Eff λ by V 2nd λ , holds for the Second-Order model. Note that Eq. ( 29) has the same form of Eq. ( 31) because for Poisson V Pois λ (t) = λ. Eq. ( 12) is in fact a general result for spike counter P belonging to the exponential family. Note that V Eff λ and V 2nd λ are properties of the noise distribution (respectively P Eff and P 2nd Θ ). These functions can be numerically estimated before running the inference of the stimulus processing model, for example by computing their values for several values of λ and then interpolating with a cubic spline. S1 Appendix. Stimulus. Our stimulus was composed of one or two black bars moving randomly on a gray background. Each bar was animated by a Brownian motion, with additional feedback force to stay above the array, and repulsive forces so that they do not overlap. The trajectories of the bars x 1 and x 2 are described by the following equations [START_REF] Mora | Dynamical Criticality in the Collective Activity of a Population of Retinal Neurons[END_REF]:

dv 1 dt = - v 1 τ + sign(x 1 -x 2 ) R |x 1 -x 2 | 6 -ω 2 0 (x 1 -µ 1 ) + σ W 1 (t) (32) dv 2 dt = - v 2 τ + sign(x 2 -x 1 ) R |x 2 -x 1 | 6 -ω 2 0 (x 2 -µ 2 ) + σ W 2 (t) (33)
where W 1 (t) and W 2 (t) are two Gaussian white noises of unit amplitude, µ 2 -µ 1 = 600µm is the shift between the means, ω 0 = 1.04 Hz, τ = 16.7 ms, R = 655µm and σ = 21.2µm • s -3/2 . The width of one bar is 100µm and the length 2 mm. For receptive field mapping, a random binary checkerboard was displayed for 1 hour at 50 Hz (check size: 60 µm).

S2 Appendix. Equilibrium Poisson process with absolute refractory period. Here we provide a sketch of the derivation to obtain the complete expression (13). We are interested in computing the probability distribution P Θ n r, τ, ∆t of the number of spikes emitted in a time window ∆t, whose starting point is chosen as random, when the inter-spike interval distribution is:

ρ Θ (t) = Θ(t -τ )re -r (t-τ ) . (34)
P Θ can be expressed as the difference between its cumulative distribution:

C(n|∆t) ≡ ∞ k=n P Θ k ∆t, r, τ (35)
P Θ n ∆t, r, τ = C(n|∆t) -C(n + 1|∆t) . (36)
Because C(n|∆t) is the probability of having at least n spikes in the time-bin ∆t it can be computed as:

C(n|∆t) = ∆t 0 dt [ρ E Θ ρ (n-1) Θ ](t) (37)
where is the convolution symbol and (•) n means n-times self-convolution. ρ E Θ is the distribution of the first spike when the beginning of the time-bin is chosen at random (equilibrium process) and its distribution can be computed as [START_REF] Müller | Dead-time problems[END_REF])

ρ E Θ (t) = ∞ t dt ρ Θ (t ) ∞ 0 dt t ρ Θ (t )
= r e -r max(t-τ,0) 1 + rτ = r Θ(t)Θ(τ -t) + Θ(t -τ )e -r(t-τ ) 1 + rτ = 1 1 + rτ rI [0,τ ] (t) + ρ Θ (t) ,

where I [a,b] (t) = 1 if t ∈ [a, b] and zero otherwise. Thanks to the explicit decomposition of ρ E Θ we have:

ρ E Θ ρ (n-1) Θ = 1 1 + rτ rI [0,τ ] ρ (n-1) Θ + ρ n Θ . (39)
In order to estimate the ρ n Θ we introduce the Laplace transform, which for a generic function h(t) reads:

L[h(t)](s) = ∞ 0 dt e -st h(t) , (40)
and use it to get read of the multiple convolutions:

ρ n Θ = L -1 [ L[ ρ n Θ ] ] = L -1 [ L[ ρ Θ ] n ] = L -1 r r + s e -sτ
After some algebra [START_REF] Müller | Dead-time problems[END_REF][START_REF] Müller | Some formulae for a dead-time-distorted poisson process: To andré allisy on the completion of his first half century[END_REF] this calculation provides Eq. ( 13).

S3 Appendix. Small f expansion of the refractory neuron model.

Here we derive the result (5), for which (3) is a particular case. We consider a general model in which the instantaneous spike rate is modulated by a time-dependent factor, α(u)r, where r is the spike rate in absence of refractoriness, with α(u ≥ τ ) = 1, and u is the time following the previous spike. In this case the inter-spike interval distribution is:

ρ α (t) = να(t) exp -ν t 0 α(t) Θ(t) , (44)
and our goal is to expand, for small τ , P α (n|∆t, r, α), the probability of having n spikes within the time-bin ∆t. At first we introduce two useful quantities:

f ≡ 1 ∆t ∞ 0 dt 1 -α(t) (45) g ≡ 1 ∆t ∞ 0 dt t 1 -α(t) (46)
Such that, if α(t) = 1 for all t > 0, then f = g = 0 and we expect to recover the Poisson case, and if α(t) = 0 for all t < τ , then f = τ /∆t and g = τ 2 /∆t 2 /2 and we expect to recover the absolute refractory case, i.e. Eq. ( 3).

Much like the absolute refractory case, we consider the cumulative distribution of P α (n|r, α, ∆t, ):

C α (n|∆t) = ∆t 0 [ρ E α ρ (n-1) α ](t) , (47)
where, as before, ρ E α is the distribution of the first spike for an equilibrium process and ρ (n-1) α is ρ α self-convoluted n -1 times. To perform the expansion in small τ , we decompose ρ E α and ρ α around an exponential distribution:

ρ E α (t) = 1 + ν 2 g ρ(t) + δρ E α (t) (48)
ρ α (t) = e r(1-A)τ ρ(t) + δρ α (t) .

(49)

In the following we will first perform the computation with ρ(t) instead of ρ E α , leaving for the end the corrections due to the factor 1 + ν 2 g and term δρ E α . This is equivalent to consider a shifted process [START_REF] Müller | Dead-time problems[END_REF] where the time-bin starts after the end of the last refractory period, instead of the equilibrium process we are considering here. As for ρ n Θ (t), see Eq. ( 41), we can use the Laplace transform to compute ρ n (t) = γ[n, r](t).

Tolhurst, D., Movshon, J., and Thompson, I. (1981). The dependence of response amplitude and variance of cat visual cortical neurones on stimulus contrast. Experimental brain research, 41(3-4):414-419. Yger, P., Spampinato, G. L. B., Esposito, E., Lefebvre, B., Deny, S., Gardella, C., Stimberg, M., Jetter, F., Zeck, G., Picaud, S., Duebel, J., and Marre, O. (2016). Fast and accurate spike sorting in vitro and in vivo for up to thousands of electrodes. bioRxiv.

Figure 1 :

 1 Figure 1: Sub-Poisson behavior of RGCs. A) Top: Example raster plot for one OFF cell. Bottom: Firing rate behavior for the same cell. Colored (resp. black) area delimits the empirical (resp. Poisson) noise, as mean +/-one s.t.d. B) As panel A, but for an example ON cell. C) Top: histogram of the observed means of the spike count across stimulus repetitions (time-bin of 16.7ms) pooled across all cells of the same type. Timebins with zero mean have been excluded. Bottom: Mean of spike count across stimulus repetitions plotted against its variance. Each point corresponds to one cell in one time bin. Multiple points are superimposed. Green full line: average of points with similar mean. Black line: prediction from a Poisson distributed spike count. D) As in C but for ON cells. Note the increase of activity with respect to OFF cells.

Figure 3 :

 3 Figure 3: Second-Order and Effective models account well for observed spike count statistics. A) The empirical behavior for the OFF population estimated from the repeated two bars stimulus (green, same data as Fig. 1) is compared with the prediction of Second-Order (purple) and Effective (red) models. Lines superimpose as the two models take a very similar form for the OFF population. For comparison, Poisson prediction is shown in black (equality line) Inset: the function θ λ of the two models (see text), is compared with θ λ = log λ for the Poisson model (black line) B) Same as A, but for the ON population. In all cases, models are learned on a separate training set.

  to be constant parameters of the model. The probability distribution (7) belongs to the class of weighted Poisson distributions and its mathematical properties have been already studied elsewhere del[START_REF] Del Castillo | Overdispersed and underdispersed poisson generalizations[END_REF].We infer γ and δ through a log-likelihood maximization on the OFF and ON cell separately. For the OFF population we obtain similar values as the Second-Order model (fig.3A). By contrast, for the ON population the Effective model takes advantage of the additional free parameter and uses it to improve its performance. We obtained the values γ * ON = -0.52 and δ * ON = 0.15, while the equivalent parameter values in the Second-Order model are f * ON -f * 2 ON = 0.15 and f * 2 ON /2 = 0.02. This Effective model is therefore a simple model able to describe accurately the mean-variance relation for different types of refractoriness, with only two parameters. 6 Benchmark of proposed models The Effective model outperforms a Poisson model at predicting the empirical spike count variance. In this section we compare its performance with two other spike count models proposed in the literature. The Generalized Count (Gen.Count) model Gao et al. (2015) can be seen as a generalization of our Effective model (with arbitrary dependence on n in the correction to Poisson, not limited to a linear combination of n 2 and n 3 , see Methods, Eq. 25) and thus offers a larger flexibility to model the spike count statistics. This however comes at the price of introducing more parameters to fit, and could potentially lead to overfitting. The Conwey-Maxwell-Poisson (COMP) model Sellers et al. (

Figure 4 :Figure 5 :

 45 Figure4: Second-Order (eq. 3) and Effective (eq. 7) model performance compared with known spike count models: COMP (eq. 26) and Generalized Count (eq. 25). Log-likelihood improvement over the Poisson for several models, all learned on a separate training set, for OFF cells (A) and ON cells (B).

Figure 6 :

 6 Figure 6: Second-Order and Effective model predicts empirical information transmission. Empirical estimation of the mutual information is compared with prediction for Poisson (black), Second-Order (pink) and Effective (red) models. Using Poisson as noise distribution for the spike count leads to a strong under-estimation, biasing the prediction of information transfer of stimulus processing models. Error-bars are standard deviation of the mean over λ, see Eq. (9).

  ) where for example, P is a Poisson distribution. The classical example for this is the Linear-Nonlinear-Poisson (LNP) model Chichilnisky (2001), but many generalizations have been proposed, especially for the retina McFarland et al. (2013); McIntosh et al. (2016); Deny et al. (2017). One of the major advantages of using a Poisson spike counter P Pois is that it allows for a straightforward optimization of the model parameters ψ Chichilnisky (2001); Mc-Farland et al. (2013). Thanks to the explicit expression of the log-likelihood (ψ) (see Methods), the log-likelihood gradient ∇ ψ (ψ) ≡ d (ψ)/dψ takes a very simple form:

Figure 7 :

 7 Figure 7: Effective model improves performance of stimulus processing models. A) PSTH for an example OFF cell (Green) is compared with predictions by the LN 2 models equipped with Poisson (black) or Effective (red) spike counter. Note that the discrepancies in the spike rate predictions are those of the LN 2 model, and not of the Poisson or Effective model, which inherit them by construction. B) As in panel A, but for the spike count variance. Poisson model largely overestimates the empirical variance during transients of high activity. C) Fano Factor as function of time: Effective model (red) accounts for the empirical behavior (green). Poisson prediction (black constant line) does not. D) Improvement with respect to Poisson in the mean square error (MSE) between empirical and model variance, plotted as a function of the MSE between empirical and model mean spike count. Each point represents a cell of the OFF population. Circled point refers to the example cell in panels A-C. E) Same as D, but for the Effective model applied on the ON population.

  have shown that a simple model taking into account refractoriness in the spiking process explains most of the deviation from Poisson that we observe in the spike count statistics. The model has only two parameters, can easily be plugged into any encoding model to be fitted to sensory recordings. It allows for an accurate estimation of the relation between the mean spike count and its variance, but also of the amount of information carried by individual neurons. The form of this model is inspired by the regularity imposed on spike trains by the refractory period. However, it can potentially work for data with other sources of regularity. In the retina, this model works for two different types of cells. Thanks to its simplicity and generality, this model could potentially be used to account for mean-variance relation in neurons recorded in other sensory

  Several solutionshave been proposed to model this over-dispersion[START_REF] Scott | Fully bayesian inference for neural models with negative-binomial spiking[END_REF];[START_REF] Goris | Partitioning neuronal variability[END_REF];[START_REF] Charles | Dethroning the fano factor: a flexible, model-based approach to partitioning neural variability[END_REF].

Figure 8 :

 8 Figure 8: Range of applicability of the model. (A) Fano factor as a function of the window size (time bin). For short windows, the Fano factor is smaller than 1 because of refractoriness, while for long windows spike rate fluctuations yield Fano factors larger than 1. (B) The autocorrelation function of the spiking activity (green) carries the footprint of both the refractory period (τ refr ) and of the time scale of spike rate fluctuations (τ λ ). The autocorrelation computed from shuffled spike trains across repetitions of the stimulus (black) only reflects the latter, as expected. (C) The Effective model we have introduced is expected to be valid when the window size ∆t is larger than the refractory period (τ refr ), but smaller than the typical time scale of rate variations τ λ . (D) Meanvariance relation for a window of 33 ms. Green curve is data, red line is the Effective model prediction.

  Count introducing the proper θ λ [G] function. For G[n] = -γn 2 -δn 3 the distribution (25) reduces to the Effective model, see Eq. (7). The Gen.Count model is thus a generalization of our model, and as such is potentially more flexible in modeling the spike count statistics. This however comes at the price of introducing more parameters to fit. In practice one needs to define G[n] for n = 0, 1, . . . , n Max , where n Max is the maximal value of n observed during the experiment. It can be shown that the model has n Max -1 free parameters Gao et al. (2015). If n Max = 1, 2, 3 then Gen.Count is equivalent to Poisson, Effective with δ = 0 and Effective models respectively. Otherwise it offers a potentially interesting way to generalize our cubic model at the price of inferring the parameters G[n] for large values of n. Thanks to its exponential form, we easily infer the Gen.Count parameters with log-likelihood maximization and we have set n Max = 4 and 5 for respectively OF F and ON populations. Conwey-Maxwell-Poisson model. In Ref. Stevenson (2016) the Conwey-Maxwell-Poisson (COMP) model Sellers et al. (2012) has been proposed to account for both under-and over-dispersed mean-variance relation. In our notation, the COMP model reads:

  26) where Z λ [η] is a normalization constant. For η = 1 the COMP reduces to the Poisson model, whereas for G[n] = -(η -1) log n! the Generalized Count reproduces the COMP. Also for the COMP, the exponential form allows for log-likelihood maximization.

  t -nτ ) (n-1) e -r(t-nτ ) (n -1)! Θ(t -nτ ) ≡ γ[n, r](t -nτ )(41)where we have introduced the Gamma distribution γ[n, r](t) ≡ Θ(t) r n t (n-1) e -r /(n -1)!. From (37) it follows: τ ] (t -t )γ[n -1, r](t -(n -1) can be computed by integrating several times the gamma distribution using the following relation:1 (n -1)! a 0 dt t (n-1) e -t Θ(t) = 1 -e -a

Acknowledgments

We like to thank R. Brette, M. Chalk, C. Gardella, B. Telenczuk and M. di Volo for useful discussions. This work was supported by ANR TRAJECTORY, the French State program Investissements d'Avenir managed by the Agence Nationale de la Recherche [LIFESENSES: ANR-10-LABX-65], a EC grant from the Human Brain Project (FP7-604102)), NIH grant U01NS090501 and AVIESAN-UNADEV grant to OM.

Code availability. Code is available online at github.com/UFerrari/Effective_ Model.

All predictions of the model were done using cross-validation. The dataset was divided into two sets, each made from one half of a 10s film repeated 54 times: a training set, on which the fit was done, and a testing set, on which the model predictions were tested.

Inference of the Second-Order refractory model. To infer the value of f for the Second-Order refractory model, see Eq. ( 3), we perform log-likelihood maximization with a steepest descent algorithm. Thanks to exponential form of P 2nd Θ (n|λ, f ), the derivative of the log-likelihood ( ) with respect to f takes a simple form:

where

At each iteration, we update the value of the parameter using the log-likelihood gradient and we adjust the function θ λ accordingly.

Inference of the Effective model. To infer the value of γ and δ for the Effective model, see Eq. (7, we perform log-likelihood maximization with Newton's method, where at each iteration we update the values of the parameters using the log-likelihood gradient:

and the Hessian:

where now • model i (t) means average with the model distribution P Eff (n|λ = λ i (t), γ, δ).

After updating the parameters, we adjust the function θ λ and iterate until convergence.

Generalized Count model. In Ref [START_REF] Gao | High-dimensional neural spike train analysis with generalized count linear dynamical systems[END_REF], the authors define a General-This allows us to perform the following expansion:

≈ e (n-1)νf γ[n, r](t)

To perform the integration of (51) we use the following approximation:

which after some algebra we obtain for P S α (n|ν), the distribution for the shifted process:

with

We need now to account for the full ρ

is of the order τ 2 we can estimate P α (n|∆t, r) = P α (n|ν) as: which is equivalent to (5).


L Meyrand
email: meyrand@lma.cnrs-mrs.fr

E Sarrouy

B Cochelin

G Ricciardi

A full PGD/HBM-based algorithm for a low-dimensional NNM continuation with modal enrichment

This work proposes a full algorithm to compute Nonlinear Normal Modes (NNMs) using a Proper Generalized Decomposition (PGD) approach combined to a continuation method. Because periodic solutions are sought, a harmonic balance method (HBM) is implemented to process the time-dependent functions. The interest of this study is first to couple this approach to a continuation method with a progressive enrichment of the PGD modes family and second to improve some of the PGD steps to get a robust algorithm. In particular, the initializations of solvers are chosen consistently with the available physical data through the use of LNMs shapes. The reduced basis is kept as small as possible since it is computed with on the fly enrichment, possibly from a unique mode. No particular assumptions on the non linearity have been made and hence a wide class of mechanical problems can be processed. The resulting method is applied on several beam models with localized nonlinearities which allows for an accurate analysis of its performance.

Introduction

The computing power currently achieved allows to take into account nonlinear effects in physical models and allows to deal with systems with a very large number of degrees of freedom. As systems became more and more realistic and sophisticated, model reduction techniques had to be developed in order to save digital resources while obtaining a physically satisfactory result.

Nonlinear Normal Modes (NMMs) are objects that can be exploited as part of modal reduction methods, similar as with Linear Normal Modes (LNMs). Their properties are studied to estimate and interpret more precisely the dynamics of mechanical systems [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Laxalde | Dynamics of a linear oscillator connected to a small strongly non-linear hysteretic absorber[END_REF]. In this paper, only the Rosenberg's definition of NNMs [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF] will be considered, i.e. a family of periodic solutions of the underlying free and conservative mechanical system. The reader can refer to the work of Shaw and Pierre [START_REF] Shaw | Non-linear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF] or Kerschen and Renson [6,[START_REF] Renson | Nonlinear Normal Modes of Nonconservative Systems[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF] for more information about NNMs.

From a numerical point of view, the algorithms used to compute the NNMs of high-dimensional structures do not have yet reached maturity [START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF] and the nonlinear effects causes convergence issues. NNMs require large sets of variables to be properly described which make them hard to compute in short times and hard to reuse even if some work still take advantage of their computations [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF]. This remark serves as a starting point for the work presented in this article.

The reducted model developed in this paper falls within the framework of Proper Generalized Decomposition (PGD) methods. This class of reduction techniques is based on the separation of the unknowns and a fixed point algorithm which allows to deal with smaller systems in order to reduce the computation times. PGD was already used by Grolet and Thouverez [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] to compute free and forced responses of nonlinear systems. The reader can refers to Chinesta and Nouy [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving timedependent partial differential equations[END_REF] to learn more about the possibilities offered by this approach.

In this work, a full algorithm is proposed in order to compute NNMs using a combination of a PGD approach and classical continuation schemes. A quick reminder of the NNM framework investigated here is given in Sec. 2. Then the PGD method is detailed in Sec. 3, including the use of a harmonic balance method (HBM) to deal with temporal aspects of the problem. The contributions of this work are highlighted in this section: first the PGD approach is coupled with a continuation method and a progressive enrichment of the PGD modes family is made and second some of the PGD steps are improved to get a robust algorithm. The PGD/HBMbased continuation scheme is described in Sec. 4. The last section, Sec. 5, is dedicated to illustrative cases to show some possibilities introduced by this algorithm.

Mechanical framework and Nonlinear Normal Modes (NNMs)

The equations of dynamics describing the behavior of a set of general coordinates can be obtained after a spatial discretization process using for example the Finite Elements method. It usually takes the following form:

M ü + C u + Ku + f nl (u, u) = f (t) (1)
Eq. ( 1) is a set of N second order nonlinear differential equations, the vector of the unknowns u -or vector of the degrees of freedom (dof) -is a vector of R N . M , C and K are respectively the linear mass, damping and stiffness matrices, with {M , C, K} ∈ M N (R) 3 . f (t) is the vector of excitation forces and lies in R N too. Finally f nl (u, u) is the vector of non linear forces and belongs to R N . The general system described by Eq. ( 1) represents a very wide class of mechanical problems. Indeed, the f nl term can be any nonlinear expression (e.g. cosine function, polynomial terms like u i u 3 j , non smooth function, etc.). This paper focuses on the particular study of NNMs, which will be considered as the periodic solutions of Eq. ( 1) in the particular case when no forcing neither damping is considered, or as an extension of the definition of LNMs which takes into account the nonlinear vector f nl . In this framework, introduced by Rosenberg, a NNM is a set of limit cycles -or periodic solutions -of the following system:

R(u(t)) = M ü + Ku + f nl (u, u) = 0 (2)
Unlike LNMs, NNMs do not decouple the equations, and there exists a frequency-amplitude dependency. However, they take into account non linear effects and can be useful to investigate modal interaction between widely spaced modes or modal bifurcations [START_REF] Renson | Nonlinear Normal Modes of Nonconservative Systems[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. Numerically, building a NNM branch requires a large number of descriptors and consequent computation times [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. This observation motivates the developments described in this work.

PGD and HBM combination to describe periodic solutions

Harmonic Balance Method

As NNMs are periodic solutions, a HBM framework is used. It does not require a time integration scheme as in shooting methods [START_REF] Sundararajan | Dynamics of forced nonlinear systems using shooting/arc-length continuation method -application to rotor systems[END_REF][START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF]. Detailed explanations about this method and its advantages are given in works [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF][START_REF] Urabe | Numerical computation of nonlinear forced oscillations by galerkin's procedure[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Sarrouy | Non-linear periodic and quasi-periodic vibrations in mechanical systems -on the use of the harmonic balance methods[END_REF]. Dealing with algebraic equations in the frequency domain has the disadvantages first to introduce a truncation order dependence and second to become costly as the number of dof increases and more harmonics are involved in the underlying physics of the system. A combination with a PGD reduction technique will permit to compute smaller systems without losing HBM advantages.

Some notations used later are here introduced. We are looking for a T -periodic signal u(t), with T = 2π/ω:

u(t) = a 0 √ 2 + H k=1 (a k cos(kωt) + b k sin(kωt)) (3)
where a k and b k are respectively the cosine and sine coefficients of the Fourier series. A compact way to write u(t) is then:

u(t) = U H h H (t, ω), where U H = [a 0 , a 1 , b 1 , . . . ] h H = [1/ √ 2, cos(ωt), sin(ωt), . . . ] T (4)
This formalism first shows computing u(t) is equivalent to compute both U H and ω, second allows to link the Fourier coefficient of u(t) to those of its temporal derivative u(t) via a matrix D:

D = diag 0, 0 1 -1 0 , . . . , H 0 1 -1 0 (5)
As HBM is a Galerkine method, this decomposition is injected in Eq. ( 2) and the equations are projected on each element of the Fourier basis h H using the inner product:

< f, g > T = 2 T T 0 f (t)g(t)dt (6)
This leads to a set of N (2H + 1) equations with N (2H + 1) + 1 unknowns, where N is the number of dofs.

As there is no external forcing in Eq. ( 2), ω is also an unknown, and the final problem can be written:

H(a 0 , a 1 , b 1 , . . . , ω) = 0 (7)
Whether a HBM is applied on the forced equation Eq. ( 1), ω is taken equal to the forcing frequency, and the system described by Eq. ( 7) is square. In the case of NNMs computation, an extra equation usually called phase condition [START_REF] Seydel | From Equilibrium to Chaos, Practical Bifurcation and Stability Analysis[END_REF] is needed. Indeed, no "initial time" exists for a given periodic solution u(t) so the phase condition allows to define the Fourier coefficients uniquely. Here are given some possibilities of equations: to lock the amplitude or the velocity of a given dof at t = 0 in Eq. (3) [START_REF] Seydel | From Equilibrium to Chaos, Practical Bifurcation and Stability Analysis[END_REF], or to nullify one of the Fourier coefficient components [START_REF] Sarrouy | Stochastic study of a non-linear self-excited system with friction[END_REF]. Considering a conservative system, please note that nullifying one sine coefficient implies that all sine coefficients are null, hence only (H+1) a k coefficients must be computed.

Considering Eq. ( 7) again, it is possible to analytically establish the expression of the projection of the linear part H l (ω) u H of Eq. ( 2), where u H = {a T 0 , . . . , a T k , b T k , . . . } T . Hence, Eq. ( 7) is rewritten H l (ω) u H + H nl (u H , ω) = 0, with:

H l (ω) = diag(Λ 0 , Λ 1 , . . . , Λ H ) Λ 0 = K, Λ k = diag(K -(kω) 2 M , K -(kω) 2 M ) (8)
The term H nl (u H , ω) corresponds to the projections of nonlinear forces f nl (u, u) onto the Fourier basis. Evaluating this nonlinear contribution is time consuming as analytical expressions can scarcely be derived. The most commun approach to deal with this term is to use the Alternating Frequency/Time (AFT) method introduced by Cameron and Griffin [START_REF] Cameron | An alternating frequency/time domain method for calculating the steadystate response of nonlinear dynamic systems[END_REF]. Basically, u(t k ) and u(t k ) are evaluated for specific t k values in [0, T ] using an Inverse Fast Fourier Transform (IFFT) based on a k and b k coefficients. Then, f nl (u(t k ), u(t k )) is computed for each time t k ; finally, a Fast Fourier Transform (FFT) is used to evaluate the projections, that is the coefficients of the Fourier series of f nl (u, u).

PGD combined to HBM

For the purpose of this work, the PGD process is divided in three main steps, which are separating variables (space and time here), obtaining as many subproblems as there are variables (time and space problems here), and run an alternated directions fixed point loop which solves each subproblem with the other variables fixed, for a given loop iteration. Here both subproblems are coupled with HBM features introduced in Sec. 3.1 to produce a global PGD/HBM algorithm. Many variants of the PGD method exist and we chose to describe the one usually called optimized PGD (oPGD) [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving timedependent partial differential equations[END_REF] in the following theoretical developments. A lighter adaptation of the oPGD from a computational point of view and called progressive PGD (pPGD) is detailed in Sec. 4.1.3.

Given a positive integer m << N , we look for a space-time separated solution u(t) of Eq. ( 2) such as:

u(t) ≈ m j=1 p j q j (t) ⇔ u(t) ≈ P q(t) with P = [p 1 , . . . , p m ] ∈ M N,m (9)
Assuming this notation, the m PGD modes are defined by their PGD mode shapes p j and their time dependences q j (t), assembled in a vector q(t) ∈ R m . In this two-term product, the amplitude information is chosen to be given by q j (t) only so the p j PGD mode shapes are normalized to 1.

This decomposition will to split the problem into two smaller problems (cf. Sec. 3.2.1 and Sec. 3.2.2) whose sizes depend on m. As NNMs are periodic solutions of Eq. ( 2), the temporal part q(t) will be decomposed using HBM:

q(t) = Q H h H (t, ω), where Q H = [a 0 , a 1 , b 1 , . . . ] h H = [1/ √ 2, cos(ωt), sin(ωt), . . . ] T (10)
Computing a NNM point by a combined PGD/HBM approach is then equivalent to compute a {P , Q H , ω} set.

Temporal problem T m

The objective of this section is to compute the temporal part -i.e Q H and ω -knowing the PGD mode shapes -i.e. the matrix P . A weak formulation is written on an oscillation period I T = [0, 2π/ω], with the test function u (t) = P q (t):

∀t ∈ I T ∀q (t)

I T q T (t)P T R(P q(t))dt = 0 (11)
The final temporal problem is established from the weak formulation Eq. ( 11):

∀t ∈ I T M r q(t) + K r q(t) + f nlr (P q(t), P q(t)) = 0 (12)
where

{M r = P T M P , K r = P T KP } ∈ M m (R) 2 and f nlr = P T f nl ∈ R m .
This system of m second order nonlinear ODEs is solved by HBM. Using the same notations as in Sec. 3.1, the temporal problem T m is defined as a square algebraic system with m(2H + 1) + 1 unknowns in the following manner:

T m (Q H , ω|P ) = 0 ⇔ H l (ω)q H + H nl (q H , ω) = 0 c(q H , ω) = 0 (13)
where q H = {a T 0 , . . . , a T k , b T k , . . . } T are the Fourier coefficients put into a vectorial form, H l is the linear contribution, H nl is the vector of Fourier coefficients of the nonlinear contribution of Eq. ( 12), and c(q H , ω) ∈ R is an arbitrary constraint equation which allows to get a square system. c(q H , ω) will be specified in Sec. 4.1 as a continuation criterion.

Spatial problem S m

The objective of this section is to compute the spatial part -i.e the matrix P -knowing the temporal parti.e. Q H and ω. This time, a weak formulation is obtained with the test function u = m k=1 p k q k , leading to a set of m systems with N equations:

∀k ∈ [[1; m]], m j=1 I T q k qj dt M + I T q k q j dt K p j + I T q k f nl (P q(t)) dt = 0 (14)
Finally, these equations can be condensed in a N × m algebraic system denoted S m :

S m (P |Q H , ω) = 0 ⇔ S l p + S nl (p) = 0 (15)
where p = {p T 1 , . . . , p T m } T ∈ R N ×m contains the columns of P , S l is the linear contribution with respect to p, and S nl (p) = [ I T q k f nl (P q(t)) dt] 1≤k≤m the nonlinear one. S nl is computed using the AFT approach, and S l can be established analitically:

S l = I 2 ⊗ M + I 0 ⊗ K , with I k = πω k-1 Q H D k T Q T H (16)
The spatial subproblem S m is more expensive to solve than T m as its size depends on N >> m.

PGD/HBM fixed point algorithm

The last step of the PGD process consists in integrating both temporal and spatial subproblems into a fixed point algorithm. The full fixed point loop is detailed in Algorithm 1. The subproblems resolution lies at lines 6 and 7. Two main issues have to be adressed: the choice of an error measure ε and the way in which potential additional PGD modes are initialized.

The error measure ε, defined at lines 1 and 11, compares the norm of the residue vector R(u(t)) -where u(t) = P Q H h H (t, ω) -to ||Ku(t)|| over a period I T . Unlike in Grolet and Thouverez work [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] in which k < k max = 3 was the only stopping criterion of the loop, here the convergence is checked by the use of this physical criterion, introduced at line 4. Moreover the fixed point loop is broken when k > k max or if the current iteration does not reduce ε enough, that is, when ε is greater than r fp ε prev where r fp < 1 and ε prev denotes the error measured at the end of the previous iteration.

If the convergence criterion ε ≤ ε max is not fulfilled at the end of the fixed point loop, the choice is made to add a new PGD mode (from line 13): a new column p j is added to P and the associated row is added to the Q H matrix. In previous works [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF], random values were used to fill these new vectors, and the order in which subproblems are dealt with in the algorithm did not seem to be important. Here, the initial observation is that it is difficult to initialize Q H with new temporal information. Hence, the temporal problem is computed first into the fixed point loop, with a new row filled with zeros. The new column p j is initialized with the shape of the (n 0 + m)-th LNM, with n 0 being the index of the n 0 -th NNM being computed. Indeed, this spatial data is easily computable and it seems relevant to add new shapes from the LNMs as they form a basis. Moreover, avoiding random initializations makes the algorithm more robust and results reproducible.

NNM continuation using PGD/HBM modes enrichment

In this section, a full NNM computation algorithm combining PGD/HBM and continuation schemes is detailed. Near the n 0 -th LNM, first NNM points are searched using a single PGD mode (m = 1) which initially has the shape of the n 0 -th LNM. When the mechanical energy is growing, new PGD modes are included to

k ← 0, ε prev ← 2 ε; 4 while k ≤ k max and ε ≥ ε max and ε ≤ r fp ε prev do / * FP loop * / 5 k ← k + 1, ε prev ← ε; 6 Q H , ω ← solutions of T m (Q H , ω|P ) = 0; 7 P ← solution of S m (P |Q H , ω) = 0; 8 for j ← 1 to m do / * Mode normalization * / 9 p j ← p j /||p j ||; end ε ← max t∈I T (||R(u(t))||/||Ku(t)||); end if ε ≥ ε max then / * Next mode initialization * / P ← [P , φ n 0 +m ], Q H ← Q H 0 1,2H+1 ; m ← m + 1; end end Return m, P , Q H , ω and ε;
the decomposition only when it is necessary along the path. The hightlights of the method are a dimension varying continuation scheme which keeps the number of NNM descriptors as small as possible and with an on the fly enrichment of the PGD modes set as the system behavior becomes more complex.

Before giving some details about the continuation scheme, we will define y n as the the n-th point on the current path. As explained in Sec. 3.2, a point of the NNM is fully determined by its {P , Q H , ω} set. Hence we will note:

y n ≡ {P n , Q H n , ω n } (size: N m + (2H + 1) m + 1) (17)
A distance d between two consecutive points y n-1 and y n is also defined. It is weighted by a set of positive constants {α P , α Q , α ω } which are arbitrary continuation control parameters, such as:

d(y n-1 , y n ) = α 2 P ||P n -P n-1 || 2 + α 2 Q ||Q H n -Q H n-1 || 2 + α 2 ω |ω n -ω n-1 | 2 (18)
The prediction and correction equations used in Sec. 4.1 both rely on this choice of distance.

Prediction/Correction continuation process

Although more complex continuation procedures exist [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF], the algorithm presented in this work can be used with a choice of simple predictor and correctors. The next point sought, denoted y , is computed from a predicted point y pred obtained with a secant method whose the arclength is set by ∆s. This predictor is chosen because of its simplicity. The correction step relies on the PGD/HBM solver described previously in Sec. 3.2 with a specific choice of constraint equation c(q H , ω) = 0, linked to a classical correction method (e.g. arclength, pseudo-arclength, etc.). The global algorithm is depicted in Algorithm 2. k max , ε max and r fp ; -LNMs: Φ = [φ n 0 , . . . , φ n ]; -Constraint equations: c 2 (q H , ω) and c p (q H , ω). Data: Circular frequency of n 0 -th mode ω n 0 .

1 y 1 ← {φ n 0 , 0 1,2H+1 , ω n 0 }; // First point 2 c(q H , ω) ← c 2 (q H , ω), m max ← m max tot ;
// Second point 3 Compute y 2 using Algorithm 1 with y 1 as initial value;

4 c(q H , ω) ← c p (q H , ω); //

Prediction step features

The prediction arclength ∆s is managed in a classical way, as follows: if the error ε measured after a correction step is too large (see Algorithm 1 in Sec. 3.2.3), ∆s is decreased by multiplying it by r div s < 1 (see line 17) leading to a closer prediction point y pred . On the contrary, if ε < r ε ε max with r ε < 1, ∆s is multiplied by r mult s > 1 (see lines 13 and 14). Else ∆s is kept constant for the next prediction step.

A secant prediction needs two solution points computed beforehand. The first point of the n 0 -th NNM is already known: y 1 ≡ {φ n 0 , 0 1,2H+1 , ω n 0 }. An amplitude condition imposed on the j-th dof at t = 0 is chosen to compute y 2 : u j (t = 0) = u j0 with u j0 = 0. This particular constraint, denoted c 2 (q H , ω) in Algorithm 2, is injected into Algorithm 1 to find y 2 starting from y 1 .

The reader could remark that during the prediction step, y n can have a greater size than y n-1 because of the possibility to get different problem dimensions m along the path. The distance given by Eq. ( 18) need summations of matrices that do not have the same sizes, and the choice of padding y n-1 missing components with zeros instead of truncating y n is made. Hence, all computed data are taken into account to reach a better prediction.

Correction step features

The startpoint of the correction step is y pred . It is injected in the PGD/HBM solver described in Algorithm 1, in which the constraint equation involved in the temporal problem T m is set to c p (q H , ω), as denoted in Algorithm 2. c p (q H , ω) defines the correction method and can be freely chosen among the classical possibilities [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. Here an arclength correction is used: c p (q H , ω ) = d(y n , y ) 2 -∆s 2 . This condition involves P which is the last value of the spatial modes returned after the spatial problem S m when processing the current temporal problem T m .

A remarkable point of this continuation scheme is the ability to increase the size of the problem when a new point y n+1 is computed from y n . At line 7 in Algorithm 2, the highest difference of size between y n and y n+1 is defined as ∆m. Consequently, two consecutive points can have different sizes, which has to be taken into account when c p (q H , ω) is evaluated. Indeed, an arclength condition requires d(y n , y ). This time, a truncation is applied for the current point y so the arclength condition is evaluated in a space which has the size of y n . Thus, no geometrical condition is imposed on the new component introduced in y .

Progressive PGD -a lighter variant

As said in Sec. 3.2.2, considering the oPGD method, S m is more expensive to solve than T m as its size is N × m. This observation motivates the use of a variant of the oPGD approach developed before, in which S m is solved as few times as possible, and in which only a few number of S m equations is solved. In order to reach these goals, severals choices are made:

• During the correction step, the spatial problem is only solved when a new PGD mode is added, i.e.

when the condition at line 13 in Algorithm 1 is fulfilled. Hence, most of the fixed point loops only solve the temporal problem. The matrix P is put as the initialization of T m . This is consistent with the fact that in practice, the nature of many points close together in the NNM branch is not sufficiently different to require the introduction of a new PGD mode shape.

• In addition to that, in order to deal with the case when S m has to be solved, a new lighter spatial problem is defined. The PGD mode shapes contained in P are definitively locked once computed, and only the new p j added at line 13 in Algorithm 1 is updated by the spatial solver. This constitutes a gain in algorithmic complexity since the new size of the spatial problem is only N . The procedure leading to the analytical expression of the spatial problem is analogous to the Galerkine process described in Sec. 3.2.2, except that this time the projection is only made on the new PGD mode shape p j . The expression of this N -sized system can be found in [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF].

This variant is usually called progressive PGD (pPGD). It allows to improve the numerical performance of the algorithm described before. However, its main disavantage is the lack of flexibility of the PGD mode shapes, which can finally generate NNM branches of higher dimension m max because more spatial information is needed when the system behavior becomes more complex.

Numerical application on two beam models

The main branch of the first NNM of a beam model with a cubic nonlinearity described in Sec. 5.1 is built using the full PGD/HBM continuation algorithm. The main results proper to the presented approach are given in Sec. 5.2. This numerical study has been performed for the presented oPGD/HBM method and its pPGD/HBM variant, but for the sake of clarity only the pPGD results are presented on the figures of this section. Finally in Sec. 5.3 some numerical considerations are made for the presented example and another similar beam model with a localized contact. Hence, two different classes of conservative nonlinearities are adressed in this paper.

Description of the beam models

The considered system is an Euler-Bernouilli bending beam with one end clamped and a spring at its free end; it is spatially discretized by FEM, as illustrated in Fig. 1. The beam properties are: Young modulus E = 210 GPa, density ρ = 7800 kg.m -3 , area of the square section S = 9 10 -4 m 2 , second area moment I = 6.71 10 -8 m 4 for an in-plane bending and length L = 1 m. The linear stiffness of the spring is k = 4 10 4 N.m -1 . The first eigenvalues of the underlying linear system consisting of the beam plus the linear spring are {ω 1 , ω 2 , ω 3 , ω 4 , ω 5 } = {217, 999, 2769, 5420, 8960} rad.s -1 . There are N e = 20 elements with two dof per node: a transverse displacement u and a rotation θ, so there are N = 40 dof.

Figure 1: Euler-Bernouilli cantilever beam (a) with a transverse cubic spring at dof 1 (b) with a transverse gap between dof 1 and the linear spring.

• On the left in Fig. 1, a polynomial non linearity is chosen by introducing a cubic stiffness k nl to the spring. The nonlinear contribution is

f nl (t) = [k nl u 3 1 (t) 0 N -1 ]
T with k nl = 9.2 10 5 N.m -3 . Otherwise, we consider an underlying linear system made of the mass and stiffness matrices of the beam plus the linear stiffness of the spring -i.e. K 11 = K 11beam + k. H = 50 harmonics are considered on an exploratory basis, in particular because the algorithmic complexity stays reasonable. Indeed, about twenty harmonics would be enough to deal with the involved physics.

• On the right in Fig. 1, a one-sided contact problem is also described: a gap g = 2 cm is introduced between the free end of the beam and the linear spring. The underlying linear system is the Euler-Bernouilli cantilever beam. The unilateral contact is modelled as follows:

f nl (t) = [k(u 1 (t) -g) 0 N -1 ] T , if u 1 (t) > g 0 N , otherwise (19)
This unsmoothed contact is quite hard as k/ 3EI L 3 = 94.56% with 3EI L 3 the equivalent stiffness of the cantilever beam taken alone. A contact problem between two identical beams would be similar from the stiffness point of view. Consistently with this class of problem, the arbitrary value H = 70 is chosen.

The results obtained for this second model will not be detailed in this paper. However, some numerical values will be given in Sec. 5.3 to extend some conclusions to this class of nonlinearities.

Parameters required by Algorithm 2 are given in Table 1. As both problems are conservative, a cosine basis is sufficient to describe the harmonic behavior of the system, so only (H + 1) Fourier coefficients have to be computed instead of (2H + 1).

Bending beam with a cubic spring at its free end

A convenient way to present the NNM branches is a Frequency-Energy plot (FEP) [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF], E being the mechanical energy of the NNM. The cubic spring case NNM1 (tangent to φ 1 at almost null amplitudes) is

∆m m max tot r ε r mult s r div s ∆s k max ε max r fp y 2 : u 1 (t = 0) {α P , α q , α ω } 1 10 0.1 1 0.5 1.5 10 10 -3 0.8 10 -3 m {10 -6 , 10 -1 , 10 -4 } (5) (0.5) (1.2) (1)
(5 10 -4 m) ({0, 10 -6 , 1})

Table 1: Algorithm 2 required parameters for the cubic spring case. The contact case parameters are given in parentheses on the second line, only when they differ from the previous ones.

displayed in a FEP given in Fig. 2. N pt is the index of converged solution points through the 131 points computed during this pPGD/HBM continuation. Only 7 PGD modes are required by the pPGD algorithm (cf. Fig. 3) to build this branch which reaches about E max = 2.25 10 5 J. The same branch can be computed with only 6 PGD modes using oPGD. Indeed the shapes of the PGD modes are always the same through a pPGD/HBM continuation process because the method involves blocking the shapes once they are integrated into the P matrix. This lack of flexibility can lead to the introduction of more PGD modes than with oPGD. Both variants obviously give identical FEP.

The interpretation of the modal interactions drawn in Fig. 2 (interaction tongue and final internal resonance) are deliberately limited in order to focus on the aspects specific to the numerical method developed. It is recalled that the NNM framework used in this paper is the Rosenberg's one [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF]. In particular, all dofs reach their maximum at the same time, hence the NNM shapes at null velocity plotted in Fig. 2 gives interesting physical data. The more energy grows the more the nonlinear spring takes importance against the bending inertia of the beam, generating localized curvature inversions near the free end. Qualitatively similar results are obtained by Kerschen for another cantilever beam with cubic spring in [START_REF] Kerschen | Computation of Nonlinear Normal Modes through Shooting and Pseudo-Arclength Computation[END_REF].

On Fig. 3, the error measure ε along the path is depicted along with the number of PGD modes. The process for adding PGD modes is clearly identified: ε reaches the limit value ε max then a new PGD mode is embedded in the calculations leading to a decrease of the error. Fig. 2 shows that the PGD modes are gradually added along the branch as expected with a smooth cubic nonlinearity. Please note that the error for the first point N pt = 1 is not displayed as it is null and cannot appear in a logarithmic scale. All these interpretations are retrieved in oPGD results. An analysis of the seven PGD mode shapes p j is provided on Fig. 4. An interesting feature proper to both pPGD and oPGD approaches is the presence of PGD modes which involve the participation of several LNMs. For instance such "combined mode" is here the PGD mode 3, which gathers contributions of the first 6 LNMs. With such an approach, high order LNMs can be taken into account in a reduced range of PGD mode shapes: here 5 PGD modes are sufficient to describe the system behavior until the interaction tongue and they involve the contributions of at least the first ten LNMs. Analoguously only 4 PGD modes are required to build the same area using oPGD, with the contributions of the first 6 LNMs. The temporal information associated with each PGD mode shape is plotted on Fig. 5. The significant variations of the temporal part through the NNM branch are expected as the previously described spatial part does not change through the pPGD/HBM continuation and the temporal part contains all the amplitude information. Only odd harmonics are plotted as all the even ones are null because of the cubic nature of the nonlinearity. As a global trend high order harmonics have less influence: at the higest energy point, all the amplitudes |a k | are sub-millimetric from k = 23, for all the PGD modes.

The PGD mode 1 has the simplest behavior as it stays mainly led by the fundamental cosine |a 1 | all along the continuation. Other PGD modes have a richer harmonic behavior with shared influences of the cosine coefficients. For example at the end of the branch, the fifteen first odd harmonics are predominant in PGD modes 2-4 and the highest harmonics are more represented by PGD modes 5-7. PGD mode shapes are combination of LNMs, as shown on Fig. 4, and this partly explains why the harmonics spectrum is complex on Fig. 5. All these qualitative interpretations also highlight the growing complexification of the beam behavior and the growing influence of high harmonics as energy grows. Here also similar results and interpretations about the temporal information are found with oPGD. Finally the main objectives of the method are met: a NNM branch is built with a few variables, the number of PGD modes is kept as small as possible by incrementing m only when the convergence criterion ε ≤ ε max is not fulfilled, and all physical data belonging to the case under study is retrieved. In the next section some details about the algorithmic performance are detailed. The presented oPGD algorithm is compared with its pPGD variant. Not only information about the the cubic spring case is given but also about the contact model described in Sec. 5.1.

Numerical considerations

Few descriptors are stored to describe the NNM branch so the main objective of reduced model approaches is reached. Table 2 shows this improvement for both PGD/HBM algorithms.

Please note that the first NNM of the contact problem has been obtained with the parameters given in Table 1 until a mechanical energy E ≈ 6.3 10 4 J. Succintly, 8 PGD modes and 115 points were required to compute the whole branch and setting ∆m = 5 allowed to directly add 3 PGD modes in a single solution point (that is when the beam first hits the spring). This shows that describing sharp or gradual changes of the system behavior by the PGD/HBM continuation scheme with on the fly enrichment is feasible. Obviously calculating a solution point via PGD/HBM is relevant if and only if m × (N + H + 1) < N × (H + 1). An easy proof gives the necessary condition for the PGD/HBM method to be more interesting than a classical HBM approach: m < H + 1. In other words, the PGD/HBM method is useful for problems which require a large number of harmonics.

In both examples the saving of variables with respect to a classical HBM can be assessed with the following ratios: 1 -m×(N +H+1) N ×(H+1) for oPGD, and 1 -m×(H+1) N ×(H+1) or 1 -m×(H+1)+N ∆m N ×(H+1)

for a pPGD solution point. The oPGD reduction is very efficient, with 78% and 81% less variables than HBM for the cubic spring and the contact cases respectively. The pPGD reduction goes further with about 88% less variables than HBM for both examples.

We focus now on the underlying trust-region solver and its iteration data used to solve the subproblems. This algorithm requires function and abscissa tolerances -set on 10 -12 here -and allows to set a maximum number of iterations N ite max . The values N ite max = 100 and N ite max = 25 are set respectively for the cubic spring and the contact cases but this limit is reached for only very few points (maxima reached when a PGD mode is about to be added or at the first turning point of the interaction tongue). In Table . 3 are given mean and maximum numbers of solver iterations for temporal and spatial problems through the continuation, for both examples and both PGD approaches. The oPGD iteration data is obtained with the same values than with pPGD (cf. Table 1) except the following ones to improve the continuation: α P = 10 -6 and y 2 : u 1 (t = 0) = 10 -3 m for the cubic spring case, and r mult s = 1.2 for the contact case. In a oPGD framework the spatial problem is heavier as expected. The pPGD allows to compute less spatial subproblems so many solution points include a zero iteration number for the spatial problem in the mean and standard deviation calculation. The cost of the space problem is thus well bypassed by this variant. Otherwise N ite max could be an interesting control parameter for bigger problems as one can keep it very small and save computation time by "converting" convergence issues into new PGD modes or a reduced arclength ∆s.

Mean

Conclusion

A full PGD/HBM-based continuation technique with on the fly enrichment is presented in this paper. It is applied to two examples of beam models representing different classes of nonlinearies. These examples show that a highly reduced description of the NNMs branches can numerically be obtained. It is to be noted that spatial initializations rely on the LNMs of the system, which gives to the PGD mode shapes a physical meaning that a classical PGD approach provide a priori. As many parameters are let free (cf. Table 1), this algorithm has a certain flexibility.

The presented approach can be used to compute bigger structures with different kinds of nonlinearities. In order to exploit the physical potential of NNMs, future prospects could include the computation of damped NNMs or forced responses of a structure.

Algorithm 2 :

 2 Continuation with on the fly PGD enrichment algorithm Parameters: Values for ∆m, m max tot , r ε , r mult s , r div s and ∆s; Values for Algorithm 1:

Figure 2 :

 2 Figure 2: (a) Main branch of NNM1 in a FEP and (b) its mode shape at null velocity points for four continuation points. Squares: solution points where a PGD mode is added.

Figure 3 :

 3 Figure 3: (a) Number of PGD modes m and (b) error measure ε against point index N pt . Squares: solution points where a PGD mode is added.

Figure 4 :

 4 Figure 4: Analysis of PGD spatial modes. (a) PGD mode shapes. Normalization is made at the maximum deflection point. (b) Participation factors of LNMs in each PGD mode. In each bar group, LNMs participations are numbered in ascending order until the tenth linear mode. p 1 is not displayed as it is LNM1 shape.

Figure 5 :

 5 Figure 5: Amplitudes of (cosine) Fourier coefficients of PGD modes. Crosses: coefficients when the PGD mode is first introduced; Circles: coefficients at the end of the branch (N pt = 131).

  Algorithm 1: PGD/HBM global algorithm Parameters: Values for m max , k max , ε max and r fp ; Definition of c(q H , ω) used in T m ; LNMs: Φ = [φ n 0 , . . . , φ n ]. Data: Initial values for m, P , Q H and ω. 1 ε ← max t∈I T (||R(u(t))||/||Ku(t)||); 2 while m ≤ m max and ε > ε max do 3

  Continuation loop 5 p ← 2; 6 while m ≤ m max tot do 7 m max ← m + ∆m;

8

Compute y using secant method with y p , y p-1 and ∆s; // Prediction 9 Compute y using Algorithm 1 with y as initial value; // Correction if ε < ε max then y p+1 ← y ; p ← p + 1; if ε < r ε ε max then ∆s ← r mult s ∆s; end else ∆s ← r div s ∆s; end end

Table 2 :

 2 Comparison between the number of descriptors needed by pPGD/HBM, oPGD/HBM and classical HBM for both examples. The reference for the N m pt values are the ones provided by pPGD/HBM computation. Only m × (N + H + 1) descriptors by solution point are required by the oPGD/HBM continuation whereas a classical HBM continuation needs N (H + 1) + 1 ones. In the case of the pPGD implementation, a single spatial problem of size N × ∆m is solved for a whole NNM segment with m PGD modes: the total number of required descriptors is further reduced to reach Σ m {N m pt × (m(H + 1) + 1) + N ∆m}, where N m pt is the number of solution points in the NNM segment with m PGD modes.

		pPGD/HBM	oPGD/HBM	HBM
		Σ m {N m pt × (m(H + 1) + 1) + N ∆m} N pt m(N + H + 1) N pt (N (H + 1) + 1)
	Cubic Spring case	33266	58695	267371
	Contact case	39161	60606	326715

Table 3 :

 3 Mean, standard deviation and maximum of iteration numbers for T m and S m problems, through the continuation.

					Standard	Maximum
				iteration number deviation iteration number
	Cubic spring case	oPGD pPGD	T m S m T m S m	6.23 16 6.56 1.29	4.89 10.92 17.96 6.23	22 54 206 46
	Contact case	oPGD pPGD	T m S m T m S m	2.18 9.77 3.56 0.35	1.53 8.68 5.52 1.92	16 49 42 14