Ste ´phanie Boule ˆtreau

Julien Cucherousset

Se ´bastien Ville ´ger

Re ´mi Masson

Fre ´de ´ric Santoul

Colossal Aggregations of Giant Alien Freshwater Fish as a Potential Biogeochemical Hotspot

The ubiquity and fascinating nature of animal aggregations are widely recognised. We report here consistent and previously undocumented occurences of aggregations of a giant alien freshwater fish, the Wels catfish (Silurus glanis). Aggregative groups were on average composed of 25 (610 SD, ranging from 15 to 44) adults with estimated average total biomass of 651 kg (386 -1132) and biomass density of 23 kg m 22 (14 -40). Aggregations always occurred within the same location. No foraging, reproductive or anti-predator behaviour were observed during the aggregations. A mass-balance model estimated that these colossal aggregations of an alien species can locally release, through excretion only, up to 70 mg P m 22 h 21 and 400 mg N m 22 h 21 , potentially representing the highest biogeochemical hotspots reported in freshwater ecosystems and another unexpected ecological effect of alien species.

Introduction

The establishment of vast groups of animals (insect swarms, fish schools, mammal herds or bird flocks) is a phenomenon that has always fascinated humans and scientists. Animal aggregation is ubiquitous, occurring in virtually all taxa, and is driven by a tradeoff whereby group members obtain benefits (protection against predators, information to optimize migration route, foraging or mate choice) that are counterbalanced by costs (e.g. intra-specific competition) [START_REF] Krause | Living in groups[END_REF][START_REF] Parrish | Complexity, Pattern, and Evolutionary trade-Offs in Animal Aggregation[END_REF]. One of the most fascinating aspects in the establishment of ephemeral animal aggregations is the occurrence of high local density and biomass. For instance, large numbers of Pacific salmon aggregate in coastal rivers to reproduce, with densities reaching 10 ind. m 22 in spawning grounds [START_REF] Schindler | Pacific salmon and the ecology of coastal ecosystems[END_REF]. Here, we report the consistent occurrence of previously undocumented and colossal aggregations of a giant alien freshwater predator (Wels catfish, Silurus glanis). This species is the world's third largest and Europe's largest freshwater fish, originates from Eastern Europe and has been introduced westward [START_REF] Stone | The last of the Leviathans[END_REF][START_REF] Copp | Voracious invader or benign feline? A review of the environmental biology of European catfish Silurus glanis in its native and introduces ranges[END_REF].

Results

Monospecific circularly-moving aggregations (Figure 1 and Movie S1) were consistently observed at the same location during 17 snorkelling surveys performed in the Rho ˆne River, France (details in Methods). The number of individuals in the aggregations was estimated to average 25 (610 SD) adults ranging from 15 to 44 (Figure 2) with estimated body size ranging from 120 to 210 cm and body weight ranging from 12 to 65 kg. Observed aggregations of 15 to 44 individuals represented an estimated total biomass of 386 to 1132 kg, corresponding to an estimated density of 0.5-1.6 ind. m 22 and an estimated biomass density of 14-40 kg m 22 (Table 1). Using a mass-balance model, we estimated that the colossal aggregations of giant alien catfish locally release, through excretion only, 21-74 mg P m 22 h 21 and 132-419 mg N m 22 h 21 (Table 1).

Discussion

Evidence suggests that the mechanism responsible for the observed aggregations were not associated with schooling behavior, reproduction, foraging or anti-predator behavior. Indeed, individuals were active, always swimming, but were not all pointing in the same direction as in polarized shoals. No synchronous movements were observed and, contrary to schooling fish that maintain a minimum distance between conspecifics, individuals were swimming while rubbing against each other. Therefore, the observed aggregations do not represent fish schools sensu stricto [START_REF] Pitcher | Heuristic definitions of fish shoaling behaviour[END_REF]. It is unlikely that they were linked to reproduction since no mating behaviours were observed and groups occurred throughout the year at temperatures below the spawning threshold [START_REF] Copp | Voracious invader or benign feline? A review of the environmental biology of European catfish Silurus glanis in its native and introduces ranges[END_REF]. In addition, foraging behaviour was unlikely since no prey was captured and no foraging behaviour was displayed. Additionally, all individuals were large enough to be released from any predation risk.

Because of the very large size of these alien individuals (approximately five times heavier than native fish species), their aggregation can potentially lead to important functional consequences in recipient ecosystems. For instance, defecation and excretion from dense aggregations of fish that rest over coral reefs can provide important quantity of nitrogen and phosphorus that subsequently increase the growth rate of corals [START_REF] Meyer | Fish Shools: An Asset to Corals[END_REF]. In some cases fish can translocate nutrients within the ecosystem by feeding in one location while defecating in another [START_REF] Meyer | Fish Shools: An Asset to Corals[END_REF][START_REF] Krone | Defecation behaviour of the Lined Bristletooth Surgeonfish Ctenochaetus striatus (Acanthuridae)[END_REF]. Additionnally, heterogenous spatial distribution of fish can also create biogeochemical hotspots, i.e. places where nutrient release by animals exceeds the need of primary producers [START_REF] Mcintyre | Fish distributions and nutrient cycling in streams: can fish create biogeochemical hotspots?[END_REF].

Here, the observed aggregations of alien fish potentially represent the highest biogeochemical hotspots ever reported for freshwater ecosystems, as our estimates correspond to 83-286 fold and 17-56 fold the maximal fish excretion values for P and N, respectively, reported in the literature (e.g. [START_REF] Mcintyre | Fish distributions and nutrient cycling in streams: can fish create biogeochemical hotspots?[END_REF][START_REF] Schaus | Nitrogen and phosphorus excretion by detritivorous gizzard shad in a reservoir ecosystem[END_REF]). Therefore, these aggregations can potentially have strong implications on ecosystem functioning since these fish may translocate nutrients from their feeding areas, concentrate locally these nutrients in the aggregation area and subsequently affect primary production and nutrient cycling. Therefore, this phenomenon represents another example of unexpected potential ecological impacts of alien species (e.g. [START_REF] Cucherousset | Ecological impacts of non-native freshwater fishes[END_REF]).

Methods

Observations

Wels catfish (Silurus glanis) were observed by snorkelling in a stretch (mean length = 400 m, mean width = 200 m) of the Rho ˆne River located downstream of Lyon (France). The studied stretch has a mean discharge of 150 m 3 s 21 , ranging annually from 50 to 200 m 3 s 21 . Seventeen surveys (2-hour long) were performed in the early afternoon by the same person (RM) at the same place, when permitted by water visibility. At this place, current is relatively slow, with a mean depth of 5 m and a substrate composed of gravels and pebbles.

Biomass and density estimates

Estimated aggregations size averaged 25 adults Wels catfish and ranged from 15 to 44 adults among the 17 surveys. The groups were mainly composed of 120-cm long individuals, followed by some 210-cm long individuals and few 170-cm long individuals (a 60-cm long individual, observed at a few occasions, was not used for the subsequent estimation of aggregation biomass, density and excretion). Estimated body-length distribution was computed as followed: 66% 120 cm TL, 14% 170 cm TL and 20% 210 cm TL. Total biomass of the aggregation was estimated assuming a sex ratio of 1:1 and using the length-weight relationships for females and males Wels catfish (Table 2). Total biomass density of the aggregation was estimated assuming that aggregations had a circular-like shape with a maximal estimated diameter of 6 m.

Excretion estimates

Nitrogen and phosphorus excretion rates were calculated from a mass-balance model [START_REF] Nakashima | The role of fishes in the regulation of phosphorus availability in lakes[END_REF] commonly used in fish bioenergetics modelling [START_REF] Hanson | Fish Bioenergetics 3.0[END_REF] as follows:

U P ~CP {G P {F P where U P is the mass of P excreted, i.e., lost in urine (g), C P is the mass of P consumed (g), G P is the mass of P allocated to growth (g) and F P is the mass of P lost in faeces (g).

Figure 1 .Figure 2 .

 12 Figure 1. Wels catfish (Silurus glanis) aggregations. Bottom-right panel illustrates the size of Wels catfish in the aggregation. doi:10.1371/journal.pone.0025732.g001

Table 1 .

 1 Mass-balance input variables and estimates.Density, total length, biomass and excretion rates were calculated for average (n = 25 individuals), minimal (n = 15 individuals) and maximal (n = 44 individuals) aggregation sizes. Excretion rates were also calculated for two diet compositions. doi:10.1371/journal.pone.0025732.t001

	Mass-balance model inputs Unit	Average Min	Max
	Density	ind. m 22	0.9	0.5	1.6
	Total length (TL)	cm	145	120	210
	Total biomass (B aggregation )	kg	651	386	1132
	Biomass density	kg m 22	23	14	40
	Mass-balance model
	estimates
	Diet: 100% fish
	P excretion rate	mmol P m 22 h 21 1366	810	2376
	P excretion rate	mg P m 22 h 21	42	25	74
	N excretion rate	mmol N m 22 h 21 17006	10201 29924
	N excretion rate	mg N m 22 h 21	241	143	419
	Diet: 80% fish & 20% crayfish
	P excretion rate	mmol P m 22 h 21 1167	692	2029
	P excretion rate	mg P m 22 h 21	36	21	63
	N excretion rate	mmol N m 22 h 21 15847	9396	27560
	N excretion rate	mg N m 22 h 21	222	132	386

Acknowledgments

We thank J. R. Britton, R. Ce ´re ´ghino, J. D. Olden, M. J. Vanni and two anonymous reviewers for their valuable comments on earlier versions of the manuscript.

The authors have no support or funding to report.

Table 2. Mass-balance model parameters.

Mass-balance model parameters Unit Value or relation References

Fish P concentration ([P] fish ) % of wet mass 0.5 [START_REF] Penszak | Phosphorus, nitrogen, and carbon cycling by fish populations in two small lowland rivers in Poland[END_REF][START_REF] Davis | Concentrations of selected elements and ash in bluegill (Lepomis macrochirus) and certain other freshwater fish[END_REF] Crayfish P concentration ([P] crayfish ) % of wet mass 0.16 [START_REF] Hanson | Fish Bioenergetics 3.0[END_REF][START_REF] Davis | Concentrations of selected elements and ash in bluegill (Lepomis macrochirus) and certain other freshwater fish[END_REF] Fish N concentration ([N] fish ) % of wet mass 2.54 [START_REF] Penszak | Phosphorus, nitrogen, and carbon cycling by fish populations in two small lowland rivers in Poland[END_REF][START_REF] Davis | Concentrations of selected elements and ash in bluegill (Lepomis macrochirus) and certain other freshwater fish[END_REF] Crayfish N concentration ([N] crayfish ) % of wet mass 1.6 [START_REF] Hanson | Fish Bioenergetics 3.0[END_REF][START_REF] Davis | Concentrations of selected elements and ash in bluegill (Lepomis macrochirus) and certain other freshwater fish[END_REF] Daily ration (DR) % of wet mass day 21 1.32 [START_REF] Copp | Voracious invader or benign feline? A review of the environmental biology of European catfish Silurus glanis in its native and introduces ranges[END_REF] Specific growth rate (SGR) % day 21 0.066 * [START_REF] Britton | Investigating the catch returns and growth rate of Wels catfish, Silurus glanis, using mark-recapture[END_REF] Assimilation efficiency of P (AE P ) % 72 [START_REF] Nakashima | Natural sources and requirements of phosphorus for fishes[END_REF] Assimilation efficiency of N (AE N ) % 80 [START_REF] Brett | Bioenergetics and growth[END_REF] Female weight (W female ) g 0.0038 TL 3.1295 [START_REF] Alp | Age and growth of the European catfish (Silurus glanis) in a Turkish Reservoir and comparison with introduced populations[END_REF] Male weight (W male ) g 0.0104 TL 2.9133 [START_REF] Alp | Age and growth of the European catfish (Silurus glanis) in a Turkish Reservoir and comparison with introduced populations[END_REF] *The specific growth rate is the mean of specific growth rates of the biggest individuals. doi:10.1371/journal.pone.0025732.t002

Excreted P (or N) was estimated as the difference between the P (N) consumed and the P (N) lost in faeces and allocated to growth. Faecal loss can be accounted for as a direct proportion of consumption [START_REF] Nakashima | Natural sources and requirements of phosphorus for fishes[END_REF] by applying a gross assimilation efficiency AE P (AE N ) for a given prey type, simplifying equation 1 to:

The mass of P (N) consumed was calculated as the product of the mass of prey consumed and the concentration of P (N) in prey tissue ([P] prey and [N] prey ). The mass of prey was deduced from the average daily ration (DR) of Wels catfish of size ranges 51-70 cm, 71-100 cm and 40-160 cm from June, July, August and September (reviewed in [START_REF] Syva ¨ranta J | Dietary breadth and trophic position of introduced European catfish Silurus glanis in the River Tarn (Garonne River basin), Southwest France[END_REF]). Two diet compositions were simulated, i.e. 100% fish or 80% fish and 20% crayfish [START_REF] Syva ¨ranta J | Dietary breadth and trophic position of introduced European catfish Silurus glanis in the River Tarn (Garonne River basin), Southwest France[END_REF][START_REF] Alp | Age and growth of the European catfish (Silurus glanis) in a Turkish Reservoir and comparison with introduced populations[END_REF].

The mass of P allocated to growth was calculated as the product of the growth rate and the concentration of P (N) in predator tissue ([P] pred and [N] pred ). The growth rate was estimated from the product of the specific growth rate (SGR). The model becomes:


Axel Hutt

Jérémie Lefebvre

Periodic external input tunes the stability of delayed nonlinear systems: from the slaving principle to center manifolds

Keywords: slaving principle, center manifold theorem, delay equations, bifurcation theory

The work illustrates a recent analysis technique that demonstrates that external periodic input affects the stability of the timeaveraged nonlinear dynamics of a delayed system. At first, the article introduces the fundamental elements of delayed differential equations and then applies these to a nonlinear delayed problem close to a transcritical bifurcation. We observe a shift of stability in the system induced by the fast periodic driving.

Introduction

Self-interactions are one of the fundamental components of complex systems. The consideration of these reentrant contributions, oftentimes exhibiting some form of latency or delay, play an important role in numerous areas of research (lasers, machining, chemistry, control), especially in models of biological systems. As such, over the last decades, retarded dynamical systems have been used to successfully describe physiological systems like the eye light pupil reflex [START_REF] Longtin | Modelling autonomous oscillations in the human pupil light reflex using nonlinear delay-differential equations[END_REF], blood circulation [START_REF] Glass | From Clocks to Chaos: The Rhythms of Life[END_REF] or postural motor control [START_REF] Boulet | Stochastic two delay-differential model of delayed visual feedback effects on postural dynamics[END_REF]. The influence of time lags is also ubiquitous in neuroscience, where various feedback loops have been exposed throughout neural circuitry. It has indeed been shown that delays play many important functional roles in neural systems and constitute one of the main mechanism underlying network synchronization and spatio-temporal activity patterns [START_REF] Hutt | Pattern formation in intracortical neuronal fields[END_REF][START_REF] Atay | Neural fields with distributed transmission speeds and constant feedback delays[END_REF][START_REF] Coombes | Bumps, breathers, and waves in a neural network with spike frequency adaptation[END_REF][START_REF] Franovic | Spontaneous formation of synchronization clusters in homogenous neuronal ensembles induced by noise and interaction delays[END_REF]. Over the past years, the question as to how spatio-temporal forcing interacts with retarded dynamics has received a vivid interest. The exact function of delays in the integration of temporally fluctuating temporal signals is still unknown, mainly because of the lack of tools the dynamical systems theory provides for the study of the stability of non-homogeneous and/or nonautonomous problems.

The center manifold theorem and the slaving principle have proven to be powerful tools in this task, and since then they have been successfully applied to many non-driven delayed problems [START_REF] Campbell | Delay Differential Equations: Recent Advances and New Directions, chapter Calculating Center Manifolds for Delay Differential Equations Using Maple[END_REF][START_REF] Redmond | Bifurcation analysis of a class of firstorder nonlinear delay-differential equations with reflectional symmetry[END_REF][START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF][START_REF] Schanz | Synergetic system analysis for the delay-induced hopf bifurcation in the wright equation[END_REF][START_REF] Haken | Synergetics[END_REF][START_REF] Schoener | The slaving principle for Stratonovich stochastic differential equations[END_REF]. The slaving principle may be seen as the physical equivalent, i.e. the representation in nature, of the major statement of the center manifold theorem. In this context, the center manifold theorem represents a promising candidate in the approach of the non-autonomous cases. However, the question whether the center manifold theorem may be extended to non-autonomous, i.e. forced, delayed problems is currently left poorly addressed. They are indeed indications that center manifolds do exist in infinite dimensional non-autonomous dynamical systems [START_REF] Chicone | Center manifolds for infinite dimensional nonautonomous differential equations[END_REF]. Furthermore, stochastic center manifold theory has been established in non-delayed noisy systems [START_REF] Boxler | A stochastic version of center manifold theory[END_REF][START_REF] Xu | On the low-dimensional modelling of Stratonovich stochastic differential equations[END_REF][START_REF] Hutt | Additive global noise delays turing bifurcations[END_REF], and a similar approach has been used in the numerical analysis of non-linear ODEs, subject of time-dependent forcing, providing accurate results. However, as the need to analyse forced system arises, it is still unfortunately unclear how to apply and compute center manifolds for non-autonomous delayed feedback systems and to expose the underlying mechanism by which the unstable modes govern the dynamics. A solution to this problem would greatly enhance the possibilities of theoretical analysis of delayed systems, of prime importance in mathematical neuroscience and related fields.

In the following, we present a method that allows to compute center manifolds of delayed system with time-dependent driving. We first outline the fundamental elements of delay-differential equations and the corresponding center manifold reduction in delayed systems, summarizing the detailed discussions of [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF][START_REF] Hale | Introduction to functional differential equations[END_REF][START_REF] Quesmi | A maple program for computing a terms of a center manifold, and element of bifurcations for a class of retarded functional differential equations with hopf singularity[END_REF][START_REF] Campbell | Analytical and symbolically-assisted investigations of hopf bifurcations in delay-differential equations[END_REF][START_REF] Faria | Normal forms for retarded functional differential equations with parameters and applications to hopf bifurcation[END_REF]. We then present results for a periodically driven delayed model with quadratic non-linearity, and show that it describes accuratly the dynamics near a transcritical instability.

General treatment of autonomous DDEs

Consider the general autonomous scalar delay differential equation,

ẋ(t) = f (x(t), x(t -τ )) = L({x(t), x(t -τ )}) + F ({x(t), x(t -τ )}), (1)
where L is a linear function with L(0)=0, and where F is a non linear and sufficiently smooth function, satisfying F (0) = DF (0) = 0. Both linear and nonlinear parts of this system may contain delayed components. In the following, we will consider the control parameter ε and investigate the stability of Eq.( 1) unfolded around the point ε = 0. It will prove convenient to work with the augmented system

ẋ(t) = f (x(t), x(t -τ ), ε) = L(x(t), x(t -τ ), ε) + F (x(t), x(t -τ )ε). (2) dε dt = 0.
This allows to take immediately into account the role played by the parameters in subsequent derivations.

The embedding

In order to consider solutions x(t) of Eq. ( 1) for t ≥ 0, one needs a complete description of initial value problem, corresponding to the retarded dynamics into the interval [-τ, 0]. This criterion implies that the map from the interval [-τ, 0] into R is not injective. Consequently, the system of Eq. ( 1) has ill-defined initial conditions, and uniqueness of solutions is not guaranteed [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF]. It is therefore imperative to consider Eq. ( 1) in an appropriate phase space which would ensure its self-consistency. To take into account the continuous dependence of the flow x(t) on the retarded dynamics, we introduce the parameter θ with -τ ≤ θ ≤ 0 and the new variable z t (θ) ∈ ℜ 2 , so that

z t (θ) ≡ (x(t + θ), ε) T . (3)
Based on this definition, an appropriate phase space can be shown to be the Banach space of continuous maps C ≡ C([-τ, 0], R × R) [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF][START_REF] Faria | Normal forms for retarded functional differential equations with parameters and applications to hopf bifurcation[END_REF]. Reformulating Eq.( 1) to take into account the continuous dependence of the flow on θ, we formally obtain

[ dx(t + θ) dt , dε dt ] T = d dt z t (θ) ≡ lim δt→0 (z t (θ + δt)) -z t (θ) δt , (4)
where we will distinguish instantaneous and retarded dynamics. For -τ ≤ θ < 0, expanding around δt = 0 to first order yields [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF]

d dt z t (θ) = lim δt→0 z t (θ) + ∂zt(θ) ∂(θ) δt -z t (θ) δt = ∂z t (θ) ∂θ . (5)
On the other hand, if θ = 0, we simply get

d dt z t (0) = ( dx(t) dt , dε dt ) T = (f (x(t), x(t -τ )), 0) T = (L[z t ] + F [z t ], 0) T . (6)
It is important to note that L[z t ] and F [z t ] are functionals. For instance, L[z t ] may be written as [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF] L

[z t ] = 0 -τ z t (θ)ω(θ)z t (θ)dθ = 0 -τ dη[θ]z t (θ), (7)
with the density function ω(θ). For the scalar DDE and a single parameter, ω(θ) is a 2x2-matrix. We may conveniently summarize these specific cases by writing the system of Eq. ( 1) in C as a infinite dimensional ODE,

d dt z t (θ) = A(z t (θ)) + X o F [z t (θ)], (8)
where the infinitesimal generator A is defined by [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF] A(z t (θ))

≡ ∂z t (θ) ∂θ + X o (L[z t ] - ∂z t (θ) ∂θ | θ=0 ]. (9)
The connection function X o (θ) allows the simultaneous treatment of the cases τ ≤ θ < 0 and θ = 0. It is defined ad hoc by

X o (θ) = 0 -τ ≤ θ < 0 I θ = 0.
Hence, the dynamics of (x(t + θ), ε) T = z t (θ) is governed by the infinite dimensional ODE in Eq. ( 8) appropriately defined in C and parameterized by θ.

Spectral analysis

Let us investigate the spectrum σ(A) of the linear operator A taken from Eq.( 9).

To this end, we consider the linearized problem of Eq. ( 8)

d dt z t (θ) = A(z t (θ)). ( 10
)
Substituting the ansatz z t (θ) = φ(θ)e λt at θ = 0 yields the eigenvalue problem,

L[z t ] -λz t (0) = (L[z t ] -λI)φ(0) ≡ ∆(λ)φ(0) = 0. (11)
Here φ(θ) is the eigenvector associated with the Lyapunov exponent λ. The spectrum of A is defined by σ(A) = {λ ∈ C|∆(λ) = 0}. For -τ ≤ θ < 0, we obtain an ODE on θ that defines the eigenvectors

λφ(θ)e λt = φ ′ (θ)e λt , (12)
where φ ′ (θ) denotes the derivative of φ with respect to θ.

A solution of Eq.( 12) is φ(θ) = φ(0)e λθ . Hence, for all the eigenvalues λ i ∈ σ(A), one finds an associated eigenvector φ i (θ), which constitutes a basis Φ(θ) = [φ 1 (θ), φ 2 (θ), ...]. This basis spans C, and as such, we can choose to write any state vector v ∈ C in terms of the eigenbasis of A.

However, it is possible that the basis Φ is neither orthogonal or normalized. Thus, consider the adjoint basis Ψ † (s) = [ψ † 1 (s), ψ † 2 (s), ...] where ψ † i (s) are eigenvectors of the adjoint linear problem

d dt y t (s) = A † (y t )) = -A(y t ) = -∂yt(s) ∂s 0 > s ≥ τ -L[y t ] s = 0. ( 13
)
Note that the adjoint problem of Eq.( 10) is found to be the foward problem, with t → -t [START_REF] Wischert | Delay-induced instabilities in nonlinear feedback systems[END_REF][START_REF] Hale | Introduction to functional differential equations[END_REF]. To ensure the bi-orthonormality of the eigenbases of A and A † , we further normalize the adjoint basis by

Ψ (θ) = (Ψ † (θ), Φ(θ)) -1 Ψ † (θ), (14)
where we introduce the bilinear form operator ( , ) in C, defined by

(a(θ), b(θ)) ≡ a(0)b(0) - θ 0 dθ ′ 0 -ξ dξ[dη(θ ′ )]a(ξ -θ ′ )b(ξ), (15)
with the measure dη ′ (θ) defined in Eq. ( 7). This bi-linear form is playing the role of the dot product in the space of functions. Normalizing the eigenbases provides

(Ψ, Φ)(θ) = I.
Considering the previous results, an arbitrary state vector v = v(θ) can be expressed in terms of the eigenbasis of (10) as

v(θ) = Φ(θ)(Ψ (θ), v(θ)).

Phase space decomposition and subspace dynamics

A well chosen decomposition of the spectrum can play a fundamental role characterizing instabilities of Eq. ( 8). Indeed, in the vicinity of an instability, we might assume, without loss of generality, that a finite number of Lyapunov exponents cross the imaginary axis while changing a control parameter and all other exponents are bounded to the left-hand plane. These bifurcating, or unstable exponents near the transition point introduce a very slow time scale, while the stable components relax much faster to their steady state. As a consequence, after a sufficiently long time, the dynamics of the system is essentially determined by the slow unstable modes: this is the essence of the slaving principle [START_REF] Haken | Synergetics[END_REF]. Indeed, one can choose σ U (A) ≡ {λ ∈ C|Re(λ) = 0} which leads to C = U + S, where S = U . The space U is the eigenspace spanned by the eigenvectors associated with unstable Lyapunov exponents. These eigenvectors constitute a basis of U , namely Φ U (θ) ⊂ Φ(θ). This implies that there exists a complementary subspace S, spanned by Φ S , associated with stable Lyapunov exponents i.e. σ S (A) ≡ {λ ∈ C|Re(λ) < 0}. Here, we label U the subspace spanned by the unstable eigenmodes i.e. for which the eigenvalue have a zero real part, which is analogous to the center subspace, but only the terminology differ. Following the discussion in the previous section, one can project the state vector z t (θ) governed by [START_REF] Campbell | Delay Differential Equations: Recent Advances and New Directions, chapter Calculating Center Manifolds for Delay Differential Equations Using Maple[END_REF] with respect to the unstable basis Φ U (θ),

U t (θ) = Φ U (θ)(Ψ U (θ), z t (θ)) = Φ U (θ)(u(t), ε) T , (16)
where (u(t), ε) T = (Ψ U (θ), z t (θ)) is a vector containing the expansion amplitudes of z t (θ) with respect to the unstable eigenbasis Φ U (θ). According to this phase space decomposition, the state vector z t (θ) may be separated into two disjoint elements, its stable and unstable components in S and U , respectively. Consequently we can write z t (θ) = U t (θ) + s t (θ), where

s t (θ) = z t (θ) -Φ U (θ)(Ψ U (θ), z t (θ)) ≡ Φ S (θ)(Ψ S (θ), z t (θ)). ( 17
)
Following this separation, we can project the dynamics of Eq.( 3) in S and U , by applying the projector [START_REF] Xu | On the low-dimensional modelling of Stratonovich stochastic differential equations[END_REF] to Eq.( 8), to obtain

Φ U (θ)(Ψ U (θ), d dt z t (θ)) = Φ U (θ)(Ψ U (θ), A(z t (θ))) + Φ U (θ)(Ψ U (θ), X o F [z t ]).
(18) Using the linearity of the bilinear form Eq.( 15), the left hand side is just

Φ U (θ)(Ψ U (θ), d dt z t (θ)) = d dt Φ U (θ)(Ψ U (θ), z t (θ)) (19)
≡ d dt U t (θ). ( 20
)
The projection of the linear operator A is simply

Φ U (θ)(Ψ U (θ), A(z t (θ))) = Φ U (θ)(Ψ U (θ), A(U t (θ))) + Φ U (θ)(Ψ U (θ), A(s t (θ))).
Then, we find with Eq.( 11),

Φ U (θ)(Ψ U (θ), A(U t (θ))) = Φ U (θ)(Ψ U (θ), Λ U U t (θ))), (21)
where Λ U is a diagonal matrix with entries being the elements of σ U (A).

Computing the term Φ U (θ)(Ψ U (θ), A(s t (θ))) uses the biorthonormality of the stable and unstable eigenbases. In the same spirit as in the case of the stable mode projection in equation Eq.( 16), the stable component s t (θ) ∈ S of the state vector z t (θ) may be written with respect to the eigenbasis of S,

S t (θ) = Φ S (θ)(Ψ S (θ), z t (θ)). ( 22
)
Using this fact, along with (Ψ U , Φ S ) = 0 yields

Φ U (θ)(Ψ U (θ), A(s t (θ))) = 0. ( 23
)
Thus, grouping the projection over the elements in U and S in Eq.( 21) and Eq.( 23) gives the projected linear component of Eq.( 8)

Φ U (θ)(Ψ U (θ), A(z t (θ))) = Φ U (θ)(Ψ U (θ), Λ U U t (θ))) = Λ U U t (θ). ( 24
)
The projection over the non-linear component X o F [z t ] of Eq.( 8) is computed from Eq.( 15), and reads

Φ U (θ)(Ψ U (θ), X o F [z t ]) = Φ U (θ)Ψ U (0)F [z t ]. (25)
Thus, combining Eq. ( 24) and Eq. ( 25) for z t (θ) = Φ U (θ)(u(t), ε) T + s t (θ), yields the dynamics of the unstable modes of Eq. ( 3)

U ∋ du(t) dt = Λ U u(t) + Ψ U (0)F [Φ U u(t) + s t (θ)] (26)
U ∋ dε dt = 0.
Using the same approach, we apply the operator (I -Φ U (θ)(Ψ U (θ), • ) to the state vector z t (θ) of Eq. ( 8) to obtain the complementary dynamics of the stable modes in

S S ∋ d dt s t (θ) = A(s t (θ)) + (X o -Φ U (θ)Ψ U (0))F [Φ U u(t) + s t ]. ( 27
)
2.4 The time-independent center manifold reduction

The manipulations described above can be seen as a procedure first identifying the stable and unstable manifolds, and secondly writing down the dynamics of Eq.( 1) in U and S explicitly. The key idea behind this projection is first that one can reduce the dynamics of the infinite dimensional system in Eq. ( 8) onto the finite dimensional center eigenspace in Eq.( 26) by applying the center manifold theorem, and second to get rid of delays. It provides a useful analysis platform on which one can investigate dynamic instabilities using the standard tools of linear analysis of ODEs.

Bifurcations are characterized by unstable Lyapunov exponents (i.e. exhibiting a zero-real part), or equivalently by a non-empty unstable eigenspace. The precise point in parameter space where σ U (A) = ∅ is called the instability threshold, and can be quantitatively described by max ε Re(λ) = 0|λ ∈ σ(A). When this critical point is reached, the center manifold theorem applies, and the stable modes in S are slaved by the dynamics of the unstable modes in U [START_REF] Haken | Synergetics[END_REF]. Then

s t (θ) = h(θ, u(t), ε), (28)
holds true. Consequently, the unstable modes are further described by the order parameter equation (OPE)

du dt = Λ U u + Ψ U (0)F [Φu + h(u)]. (29)
While in the vicinity of an instability and given that the functional h(u) is known, this system captures the dynamics of Eq.( 1) entirely. Delayed components are not present anymore, and the dimensionality of this representation is finite, making the OPE very useful in the treatment of non linear DDEs. Most importantly, it is possible to reconstruct the flow x(t) of the original delayed system of Eq.( 1) solely from the unstable modes by

x r (t) = Φ U (0)u(t) + h(0, u(t), ε). ( 30
)
Although the center manifold theorem ensures that the functional h(θ, u(t), ε) exists, it is typically not unique and usually challenging to compute explicitly. Such a derivation is often realized using algebraic manipulation softwares, and other methods [START_REF] Quesmi | A maple program for computing a terms of a center manifold, and element of bifurcations for a class of retarded functional differential equations with hopf singularity[END_REF][START_REF] Campbell | Analytical and symbolically-assisted investigations of hopf bifurcations in delay-differential equations[END_REF].

The reconstructed flow of x r (t) converges to the original flow x(t) as the manifold h(θ, u) gets closer to its exact form. Since ds/dt = D u h(u)(du(t)/dt) from Eq.( 28), the center manifold satisfies the implicit relationship

D u h(θ, u, ε)[Λu+Ψ U (0)F [Φu+h(u)] = A(h(u))+(X o -Φ U (θ)Ψ U (0)) F [Φ U u + h(u)] , (31)
taking into account Eqs. ( 26) and [START_REF] Yeganefar | Input-to-state stability of time-delay systems : A link with exponential stability[END_REF]. Here D u denotes the partial derivative with respect to u.

A typical ansatz to compute h(θ, u, ε) is a polynomial expansion in powers of u and the control parameter ε. Then sorting the terms by orders of ε and u yields a set of first order linear differential equations in each of the polynomial coefficients for -τ ≤ θ < 0. The initial conditions of these are fixed by solving Eq.( 31) for θ = 0. The dimensionality of the differential equations in each the coefficients is the same as the codimension of the bifurcation considered. With this ansatz, one can compute h up to any desired accuracy, by computing higher orders coefficients in the expansion and , hence, make x r (t) as close as desired to x(t). However, for most applications an expansion to low order is sufficient.

The time-dependent center manifold reduction

Now consider the non-autonomous DDE

ẋ(t) = f (x(t), x(t -τ ), t) = L({x(t), x(t -τ )}) + F ({x(t), x(t -τ )}) + I(t). (32)
Equation ( 32) is a non-autonomous delay-differential equation and it is a challenging problem to find conditions for its stability. Close to a bifurcation point, the analysis of such DDEs has attracted increasing attention in the last years, e.g. considering more general [START_REF] Hale | Introduction to functional differential equations[END_REF][START_REF] Asl | Analysis of a system of linear delay differential equations[END_REF], deterministic [START_REF] Li | Forced phase-locked response of a nonlinear system with time delay after hopf bifurcation[END_REF][START_REF] Xu | Effects of time delayed position feedback on a van der pol 2013duffing oscillator[END_REF] or stochastic forces [START_REF] Frank | Stationary solutions of linear stochastic delay differential equations: Applications to biological systems[END_REF][START_REF] Guillouzic | Small delay approximation of stochastic delay differential equation[END_REF]. The approach discussed in the previous sections cannot, formally, be used since the origin is no longer a fixed point and the eigenbases definition and associated phase space decomposition are not valid anymore. In addition, by virtue of the new time scales introduced by the external input, it is more difficult to identify separate time scales which is necessary in the center manifold technique. However, when the driving I(t) is small and fast compared to the relatively slow unstable modes, one may consider the fixed point of the autonomous system of Eq. ( 1) for the analysis of Eq. ( 32). This step is reasonable since previous studies on nonlinear delayed systems have shown that such driven delayed systems are stable under certain conditions in the sense of Input-to-state Stability [START_REF] Yeganefar | Input-to-state stability of time-delay systems : A link with exponential stability[END_REF]. One can then use the spectrum and subspace eigenbases from the autonomous system to project the dynamics of the non-autonmous system of Eq. ( 32) and subsequently obtain the non-stationary version of Eqs. ( 26)-( 27)

U ∋ du(t) dt = Λ U u(t) + Ψ U (0)F (Φ U (θ)u(t) + s t (θ)) + Ψ U (0)I(t) (33) S ∋ d dt s t (θ) = As t (θ) + (X o -Φ U (θ)Ψ U (0))[F (Φ U (θ)u(t) + s t (θ)) + I(t)].
This result has been demonstrated formally in the case of linear non-autonomous delayed system [START_REF] Hale | Introduction to functional differential equations[END_REF]. Hence the approximation is reasonable since the amplitude of both the stable and unstable modes in the vicinity of an instability can be taken arbitrarily small by adjusting the control parameter, making the non-linear component F small enough. In other words, close to the origin, the nonlinear dynamics of the system is close to its linear dynamics.

To reduce the dimensionality of the system from infinity to few modes, we proceed with the asumption that the center manifold theorem still applies close to the instability and that the functional h exhibits an explicit time dependence now. This assumption has been considered successfully in non-delayed systems for quasi-periodic input [START_REF] Cox | Center manifolds of forced dynamical systems[END_REF]. Moreover the existence of time-dependent center manifolds in non-delayed systems has been proven for stochastic inputs I(t) [START_REF] Boxler | A stochastic version of center manifold theory[END_REF]. Accordingly we choose

s t (θ) = h(θ, u(t), ε, t) . (34)
As in the case of autonomous systems, the functional h has yet to be at least approximated, to be of any use in the analysis of Eq.( 32). In particular, now the manifold h(θ, u, ε, t) satisfies

D u h(θ, u, ε, t)[Λu + Ψ U (0)F [Φ U u(t) + h(u, ε, t)] + Ψ U (0)I(t)] + ∂h(θ, u, ε, t) ∂t = A(h(u, ε, t)) + (X o -Φ U (θ)Ψ U (0)) (F [Φ U u(t) + h(u, ε, t)] + I(t)) . (35)
As an ansatz, we add a time-dependent correction h t to the expansion used in the autonomous case, such that the time-dependence in the center manifold takes the form of a fast additive perturbation

h(θ, u, t, t ′ ) = h n (θ, ε, u) + h t (θ, t) + O(m > n; u, ε, t), (36)
where O(m; u, ε, t) denotes terms of order of magnitude m in u, ε and the timedependent contribution. The ansatz of Eq. ( 36) assumes that the autonomous center manifold h n of order n and the correction term h t have similar order of magnitude. In the following we assume O(h t ) = 2 and the order of h n may be n = 2 or n = 3. This ansatz implies time-dependent corrections that are small compared to the amplitude of the unstable modes. This ansatz is analogous to the one used by [START_REF] Cox | Center manifolds of forced dynamical systems[END_REF] for forced non-linear ODEs, which proved to accurately reproduce the dynamics for various types of driving.

Moreover, Eq. ( 36) assumes that, to second order, the center manifold has a separable form in time t and modes u which facilitates the resolution of the resulting ODE system for -τ ≤ θ < 0. Indeed, the substitution of this ansatz in Eq. ( 35) up to quadratic order, i.e. n = 2 leads to the same set of differential equations as the autonomous problem, with the notable exception of an additional slow equation in h t (θ, t) of order 2. The terms in h t are decoupled from the autonomous contribution h n (θ, u) and obey

∂h t (θ, t) ∂t = ∂h t (θ, t) ∂θ -Φ(θ)Ψ (0)I(t) . ( 37
)
Here we chose O(I(t)) = 2 Equation ( 37) is a linear first order non-homogeneous partial differential equation of the time correction coefficient h t , which may be solved using the method of characteristics, given that I(t) and the entries of the bases Φ(θ) and Ψ (θ) are smooth enough. To solve Eq. ( 37), we have to distinguish the two cases:

for t + θ ≤ 0, Eq. ( 37) is an intial value-problem with the history function g(t) , -τ ≤ t ≤ 0, i.e. h t (θ, 0) = g(θ). Then the method of characteristics leads to

h t (θ, t) = - t 0 Φ(t + θ -s)Ψ (0)I(s)ds + H(θ + t) , t + θ < 0, ( 38
)
with H(θ) = g(θ). We point out that this solutins holds for the time interval t ∈ [0; -θ] only. -for t + θ > 0, Eq. ( 37) is a boundary value-problem at θ = 0 and we find by the method of characteristics

h t (θ, t) = - -θ 0 Φ(-s)Ψ (0)I(t + θ + s)ds + H(θ + t) . (39)
Indeed, writing Eq.( 35) for θ = 0 yields

∂h t (0, t) ∂t = L[h t (t)] + (1 -Φ(0)Ψ (0)) I(t) . (40)
According to Eq.( 39), we may write h t (θ, t) = r(θ, t) + H(t + θ), where r(θ, t) = --θ 0 Φ(-s)Ψ (0)I(t + θ + s)ds. Then inserting this expression into (40) leads to

∂H(t) ∂t = L[H(t)] + L[r(t)] + (1 -Φ(0)Ψ (0)) I(t) . (41)
Recall that L[H(t)] is a functional of H(θ + t) and thus L[H(t)] may depend on H(t) and H(tτ ). Similarly L[r(t)] is a functional of r(θ, t). Given that the functional r(θ, t) is known and the linear term contains retarded terms, Eq.( 41) is a non-autonomous linear delay-differential equation in H which can be solved analytically [START_REF] Hale | Introduction to functional differential equations[END_REF][START_REF] Amann | Some basic remarks on eigenmode expansions of time-delay dynamics[END_REF].

The approach illustrates the hypothesis that the non-autonomous case can be analyzed, in lowest order approximation, by a fast and small time-dependent correction on the center manifold. Higher degrees of accuracy than the second order could be achieved by proceeding to higher order terms in both modes and time dependent components of h(θ, u, ε, t) in Eq.( 36). The time-mode separable form, i.e. the separation of the mode-dependent part h a and the time-dependent part h t , combined with the time-scale separation asumption allows to compute higher order terms in the modes expansion, while keeping the time-dependent component to second order. The subsequent section examines the hypotheses made and the corresponding results by applying the method to a specific example.

The asymmetrical transcritical bifurcation

To validate our approach, let us apply the procedure discussed in the previous sections to a non-autonomous delayed differential equation with quadratic nonlinearity near a codimension 1 instability. The system we consider reads

dx(t) dt = -x(t) -R 1 x(t -τ ) -R 2 x 2 (t -τ ) + I(t), (42)
with the augmented system

dx(t) dt = -x(t) + x(t -τ ) -εx(t -τ ) -R 2 x 2 (t -τ ) + I(t) dε dt = 0, (43)
where we introduced the control parameter ε ≡ R 1 + 1. Applying the steps described in the previous section we obtain the reduced dynamics [START_REF] Lefebvre | Reduced dynamics for delayed systems with harmonic or stochastic forcing[END_REF]

du(t) dt = 1 1 + τ F [u(t) + s t ] + 1 1 + τ I(t) dε(t) dt = 0 (44) d dt s t (θ) = A(s t ) + X o - 1 1 + τ (F [u(t) + s t ] + I(t)).
Now applying the center manifold theorem implies that the functional h depends on time explicitely and with the separation ansatz 36 h t obeys

∂h t ∂t = ∂h t ∂θ - 1 1 + τ I(t) . (45)
Then Eqs.( 38), (39) lead to the solutions

h t (θ, t) = - 1 1 + τ t 0 I(s)ds + H(θ + t) , t ≤ -θ (46) h t (θ, t) = - 1 1 + τ -θ 0 I(t + θ + s)ds + H(θ + t) , t > -θ,
with the initial condition H(θ) = g(θ). Substituting these results into Eq.( 41) yields the evolution equation of H(t), t ≥ 0 for θ = 0

dH dt = -H(t) + H(t -τ ) - 1 1 + τ τ 0 I(t -τ + s)ds + τ 1 + τ I(t) . (47)
We point out, that the last two terms may be viewed as an external driving and the linear terms are the same as the linear terms in the original system. By virtue of the spectrum of the linear operator the linear system in H(t) in ( 47) is marginally stable close to the stability threshold (the maximum Lyapunov exponent is close to zero). The other Lyapunov exponents have negative real parts and their contribution vanish for large times.

To verify these results, let us consider the periodic driving I(t) = I o sin(w o t), whose amplitude I o is small compared to the amplitude u(t) and whose oscillation period is short compared to the slow evolution of u(t). The autonomous components of the center manifold were found previously. The time-dependent correction is given by Eq.( 46)

h t (θ, t) = H(t + θ) + I o w o (1 + τ ) (cos(w o t) -1) , t ≤ -θ (48)
h t (θ, t) = H(t + θ) + I o w o (1 + τ ) (cos(w o t) -cos(w o (t + θ))) , t > θ .
Assuming the initial function g(t) = 0, -τ ≤ t ≤ 0, and for large times t → ∞, the solution of (47) reads Taking into account the time-dependence of the center manifold gives access to more than just the good reconstruction of the systems dynamics. In general, nonautonomous components play a major role in the stability of dynamical systems, especially in the vicinity of dynamic instabilities, i.e. in the presence of different time scales [START_REF] Hutt | Additive global noise delays turing bifurcations[END_REF][START_REF] Hutt | Additive noise-induced turing transitions in spatial systems with application to neural fields and the swift-hohenberg equation[END_REF]. Therefore, the study of time-corrected center manifolds yield details about input-induced bifurcations. Let us investigate the interaction of the input I(t) with the transcritical bifurcation studied.

H(t) = - I o τ (1 + τ ) 2
Based on the calculations in the previous section, we may write the order parameter equation as

du(t) dt = 1 1 + τ -ε(u + h(-τ, u, t)) -R 2 (u + h(-τ, u, t)) 2 + I(t) , (50)
where we omitted the trivial dynamics of the control parameter ε. Over a finite time interval, we might time-average both sides of Eq.( 50) to obtain

du(t) dt = 1 1 + τ ( I(t) -ε u -ε h(-τ, u, t)) -R 2 u 2 -2R 2 uh(-τ, u, t)) -R 2 h(-τ, u, t) 2 , ( 51
)
where the time averaging operator • is defined over an interval T by

• = 1 T t+T t ( • )dt.
Following the time scale separation issued by the center manifold theorem and the OPE, the input I(t) is considered fast, compared to the unstable mode u(t). Hence, for T sufficiently small, u(t) is approximately constant, and we may consequently write u(t) ≈ u(t). Thus, if the input is chosen such that I(t) = 0, Eq.( 51) may be written as

du(t) dt = - 1 1 + τ R 2 u 2 + ε R 2 + 2 h(-τ, u, t) u + ε R 2 h(-τ, u, t) + h(-τ, u, t) 2 . ( 52
)
Because of the separation of the time scales, the fixed points of the averaged equation correspond to the fixed points of Eq.( 50). Inserting the stable manifold h(-τ, u, t) with the terms in Eqs. ( 48) and (49) found in section 3, we gain the stationary states by setting du(t) dt = 0. In addition the focus to the solutions close to the origin, leads to the stationary states

u o = 1 2A 0 + 2A 1 h t -(B 0 + B 1 h t ) ± (B 0 + B 1 h t ) 2 -4(A 0 + A 1 h t ) h 2 t , , (53)
with functions A 0 = A 0 (ε), A 1 = A 1 (ε), B 0 = B 0 (ε), B 1 = B 1 (ε) and h t = h t (-τ, t). Since the external input can be viewed as a a linear superposition of oscillations according to Fourier theory and hence h t exhibits oscillations with the same frequencies, it is reasonable to assume that h t ≈ 0. This assumption implies that the time-average window T is smaller than the period of the slowest Fourier component in the input. Then we find

u o = 1 2A 0 -B 0 ± B 2 0 -4A 0 h 2 t . (54)
with h 2 t ≥ 0. If h 2 t = 0, then I(t) = 0 and u o = 0, -B 0 /A 0 ≈ -ε × const and the origin u o is the only stationary solution for ε = 0. In contrast, if h 2 t > 0 for I(t) = 0 and the origin is not a stationary solution of the dynamics. New equilibria are moved to ε min,1 > 0, ε min,2 < 0. These solutions are roots of the polynomial B 2 0 (ε) -4A 0 (ε) h 2 t = 0. These results demonstrate that the external input destroys stationary states, that existed without external input, and breaks the symmetry of the transcritical bifurcation: the external input changes the stability of the system.

To verify these analytical results, we choose the averaging interval as the period of one input cycle i.e. T = 2π/w o so that I(t) = 0. Considering the full terms in Eq. ( 52) and by setting du(t) dt = 0, we may find the stationary states numerically, see Fig. 1. We immediately see that whenever h 2 t > 0, the symmetry of transcritical bifurcation is broken, and we obtain a imperfect bifurcation scenario. This symmetry breaking replaces the intersecting branches by two disjoint saddle-node curves. In order for the order parameter equation to capture this particular bifurcation diagram, the precision of the center manifold h(θ, t, u) is very important. The time-corrected center manifold brings a considerable amount of accuracy to the OPE, not only by adjusting the phase but also the amplitude of the system's response. Fig. 1 shows how the OPE with time-dependent center manifold reproduces the bifurcation diagram of the original DDE with an improved accuracy compared to the same problem but without any time-dependency on the center manifold. 53) delimits the basins of attraction of the original system [START_REF] Bechhoefer | Feedback for physicists: A tutorial essay on control[END_REF] shown in tones of shaded gray for different input amplitudes. The input amplitudes have been set to Io = 0.05 (stationary solutions 1 and basin of attraction in light gray), Io = 0.1 (stationary solutions 2 and basin of attraction in gray) and Io = 0.3 (stationary solutions 2 and basin of attraction in dark gray). As the input amplitude increases i.e. Io > 0, the basin of attraction of the stable fixed points splits and exhibit a band of unstable initial conditions, indicating that the input induces an imperfect transcritical bifurcation. In this case, the stable and unstable branches do not meet at ε = 0 as expected and are replaced by two saddle nodes bifurcations. (b) Plot of the fixed point curves predicted by the averaging of the original system [START_REF] Bechhoefer | Feedback for physicists: A tutorial essay on control[END_REF] where a standard transcritical case is predicted. This result does not correspond to the dynamics of the original system. In contrast, the fixed point curves (53) of the averaged order parameter equation using time-dependent center manifold show saddle node bifurcations for different input amplitudes, delimiting the basins of attraction of the original system accurately for Io = 0.05 and Io = 0.1. For Io = 0.3, the input amplitudes becomes large compared to the unstable mode amplitude. Additional parameters are wo = 15, τ = 2.0, R2 = 1.5

Concluding Remarks

In this essay, we showed that the dynamics of a non-autonomous delayed feedback system could be captured by center manifold reduction. This is made possible by allowing an explicit time dependence of the manifold, taking the form of an additive time-dependent correction to the non-driven problem. We illustrated the approach by considering a scalar delay differential equation with quadratic non-linearity, driven by an additive time-periodic term, in the vicinity of a transcritical bifurcations. Numerical experiments are in good agreements with the analytical results. Higher degrees of accuracy could be reached by considering higher order terms in both time and mode dependent components of the center manifold.

It is still unclear how the initial conditions of the original DDE are mapped to those of the OPE. There appears to be discrepancy at t = 0 between the reconstructed and original flow of the system considered in this example, that we corrected manually to match the initial conditions of both the OPE and DDE. This deviations seems to be due to the yet unknown map from the interval [-τ, 0] to the initial value problem x(0) = x o , induced by the projection onto stable and unstable subspaces. This discrepancy usually decays numerically with the transients. A tentative solution to this problem would require a consideration of the individual stable modes initial value problems, an information that appears to be lost with the application of the center manifold theorem and following approximations.

The time-dependent correction considered here is appropriate for additive driving only, as no time-mode mixing is present at lowest order. We might consequently expect that mutiplicative time-dependent forcing would require a different working ansatz, which would not allow a separation between autonomous and non-autonomous problems, as the example detailed here shows. This case would invariably lead to more involving calculations.

  w o (cos(w o t -1)) + τ I o w 2 o (1 + τ ) 2 (sin(w o t)sin(w o (tτ ))sin(w o τ )) +R 1 (I o , w o ) sin(w o t) + R 2 (I o , w o ) cos(w o t). (49) with constants R 1 (I o , w o ) and R 2 (I o , w o ) depending on the stimulus frequency w o , the input strength I o and the delay τ .

Fig. 1 .

 1 Fig. 1. Bifurcation diagram of the averaged order parameter equation. (a) Comparison of the basins of atttaction of the original DDE for different Io with those predicted applying the proposed time-dependent center manifold reduction. The fixed point curves (solid and dashed lines) of the averaged OPE in Eq.(53) delimits the basins of attraction of the original system[START_REF] Bechhoefer | Feedback for physicists: A tutorial essay on control[END_REF] shown in tones of shaded gray for different input amplitudes. The input amplitudes have been set to Io = 0.05 (stationary solutions 1 and basin of attraction in light gray), Io = 0.1 (stationary solutions 2 and basin of attraction in gray) and Io = 0.3 (stationary solutions 2 and basin of attraction in dark gray). As the input amplitude increases i.e. Io > 0, the basin of attraction of the stable fixed points splits and exhibit a band of unstable initial conditions, indicating that the input induces an imperfect transcritical bifurcation. In this case, the stable and unstable branches do not meet at ε = 0 as expected and are replaced by two saddle nodes bifurcations. (b) Plot of the fixed point curves predicted by the averaging of the original system[START_REF] Bechhoefer | Feedback for physicists: A tutorial essay on control[END_REF] where a standard transcritical case is predicted. This result does not correspond to the dynamics of the original system. In contrast, the fixed point curves (53) of the averaged order parameter equation using time-dependent center manifold show saddle node bifurcations for different input amplitudes, delimiting the basins of attraction of the original system accurately for Io = 0.05 and Io = 0.1. For Io = 0.3, the input amplitudes becomes large compared to the unstable mode amplitude. Additional parameters are wo = 15, τ = 2.0, R2 = 1.5