Guillaume Girier
email: ggirier@bcamath.org

Mathieu Desroches

Serafim Rodrigues

From integrator to resonator neurons: A multiple-timescale scenario

Keywords: neuronal dynamics, excitability, multiple timescales, integrator neuron, resonator neuron, folded singularities, canards

Neuronal excitability manifests itself through a number of key markers of the dynamics and it allows to classify neurons into different groups with identifiable voltage responses to input currents. In particular, two main types of excitability can be defined based on experimental observations, and their underlying mathematical models can be distinguished through separate bifurcation scenarios. Related to these two main types of excitable neural membranes, and associated models, is the distinction between integrator and resonator neurons. One important difference between integrator and resonator neurons, and their associated model representations, is the presence in resonators, as opposed to integrators, of subthreshold oscillations following spikes. Switches between one neural category and the other can be observed and/or created experimentally, and reproduced in models mostly through changes of the bifurcation structure. In the present work, we propose a new scenario of switch between integrator and resonator neurons based upon multiple-timescale dynamics and the possibility to force an integrator neuron with a specific time-dependent slowly-varying current. The key dynamical object organising this switch is a so-called folded-saddle singularity. We also showcase the reverse switch via a folded-node singularity and propose an experimental protocol to test our theoretical predictions.

Introduction

Excitable systems, in particular neurons, can be classified according to the various criteria, one of them being the existence of sub-threshold oscillations [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF][START_REF] Prescott | Pyramidal neurons switch from integrators in vitro to resonators under in vivo-like conditions[END_REF]. This feature allows to distinguish between two types of neurons: integrator and resonator . Integrator neurons are defined by : a) the ability to get excited under high frequency pulses, b) the existence of a precise threshold, and c) the fact that they do not have a sub-threshold oscillations. They belong to what is usually referred to as type-I neuronal excitability, which means that their firing frequency starts from 0 at the transition between the stationary and the periodic regime. In contrast, resonator are neurons that : a) respond only to input with well-defined frequencies (they "resonate" with these special frequencies), b) do not have a well-defined threshold, and c) have sub-threshold oscillations. Therefore they display type-II neuronal excitability, which means that their firing frequency is bounded away from 0 at the transition between the stationary and the periodic regimes. In an experimental framework, neurons behave either as integrators or as resonators, and this feature is observed neither simultaneously nor in the same physiological conditions.

From a dynamical systems standpoint, the underlying models of these two types of neurons differ by their bifurcation structure upon variation of an applied current I as main parameter. Namely, in integrator-type models a saddle-node on invariant circle (SNIC) bifurcation organises the transition from rest to spiking, and their excitability threshold is defined by the stable manifold of the saddle equilibrium that disappears through the SNIC bifurcation. In contrast, in resonator-type neuron models this transition occurs via a Hopf bifurcation (HB), which is often subcritical and followed, in parameter space, by a saddle-node bifurcation of limit cycles [START_REF] Rinzel | Analysis of neural excitability and oscillations[END_REF]. The threshold is not well-defined however it can be approximated by a family of so-called canard cycles [START_REF] Benoît | Chasse au canard[END_REF][START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF][START_REF] Desroches | Inflection, canards and excitability threshold in neuronal models[END_REF][START_REF] Wechselberger | Canard Theory and Excitability[END_REF].

In in vitro experiments, it is currently possible to make an integrator-type neuron behave like a resonator neuron by means of pharmacological intervention. Indeed, the notion of behavioural switch between integrator and resonator neuron has long been described in the experimental and computational neuroscience literature, however using different approaches. At the experimental level, the environment of the neuron can be controlled in order to obtain this change of excitable behaviour. In particular, this has been achieved pharmacologically, to control the opening of ion channels [START_REF] Prescott | Pyramidal neurons switch from integrators in vitro to resonators under in vivo-like conditions[END_REF][START_REF] Kirst | Fundamental structure and modulation of neuronal excitability: Synaptic control of coding, resonance, and network synchronization[END_REF], by current injection [START_REF] Fellous | Frequency dependence of spike timing reliability in cortical pyramidal cells and interneurons[END_REF], using an electric field [START_REF] Yi | Exploring how extracellular electric field modulates neuron activity through dynamical analysis of a two-compartment neuron model[END_REF] or even by means of an excitation laser in neuromorphic experiments [START_REF] Feldmann | All-optical spiking neurosynaptic networks with self-learning capabilities[END_REF].

In mathematical and computational studies, these methods have also been demonstrated as viable [START_REF] Prescott | Pyramidal neurons switch from integrators in vitro to resonators under in vivo-like conditions[END_REF][START_REF] Feldmann | All-optical spiking neurosynaptic networks with self-learning capabilities[END_REF][START_REF] Roach | Acetylcholine mediates dynamic switching between information coding schemes in neuronal networks[END_REF][START_REF] Al | M-current induced bogdanov-takens bifurcation and switching of neuron excitability class[END_REF][START_REF] Kirst | Fundamental structure and modulation of neuronal excitability: Synaptic control of coding, resonance, and network synchronization[END_REF][START_REF] Yi | Exploring how extracellular electric field modulates neuron activity through dynamical analysis of a two-compartment neuron model[END_REF][START_REF] Guantes | Variability in noise-driven integrator neurons[END_REF], together with other approaches: to name a few, by adding terms taking into account new neuronal structures [START_REF] Kirst | From integrator to resonator: The effect of dendritic geometry on neuronal excitability[END_REF][START_REF] Zhao | Transitions between classes of neuronal excitability and bifurcations induced by autapse[END_REF], by varying some of the model's parameters [START_REF] Hutcheon | Resonance, oscillation and the intrinsic frequency preferences of neurons[END_REF][START_REF] Macherey | A dualprocess integrator-resonator model of the electrically stimulated human auditory nerve[END_REF], or by varying the input forcing frequency [START_REF] Fellous | Frequency dependence of spike timing reliability in cortical pyramidal cells and interneurons[END_REF][START_REF] Muresan | Resonance or integration? self-sustained dynamics and excitability of neural microcircuits[END_REF]. However, all these methods have in common that they result in changes of the system's bifurcation structure in order to allow this transition from integrator to resonator behaviour. In the present case, we want to keep the same bifurcation structure, that is, a SNIC bifurcation associated with integrator-type behaviour, and act differently upon the system so that it can be made to display the characteristic subthreshold oscillations of a resonator.

The main objective of the present work is to demonstrate mathematically that an integrator neuron -namely, a type-I neuron model-can be made to behave like a resonator neuron once an adequate slowly-varying current is applied to it with real-time feedback from the membrane potential; hence, we aim to obtain subthreshold oscillations in an integrator neuron model. Crucially, we want to achieve this apparent excitability switch without modifying the underlying bifurcation structure of the model.

To do so, we will exploit the multiple-timescale structure of the slowly-forced integrator model and show that subthreshold oscillations are possible in a specific parameter range, no matter which integrator model we are starting from, provided it has a SNIC bifurcation upon constant applied current and provided we apply to it a specific slowly-varying timedependent current.

In a nutshell, we will show that the forcing requires to consider two additional slow variables and that the extended (minimally 4D) model possesses a so-called folded-saddle singularity [START_REF] Desroches | Mixed-Mode Oscillations with Multiple Time Scales[END_REF][START_REF] Desroches | Spike-adding in parabolic bursters: The role of folded-saddle canards[END_REF]. It was recently discovered [START_REF] Mitry | Folded saddles and faux canards[END_REF] that subthreshold oscillations can appear near a folded saddle provided a certain algebraic condition is satisfied in the singular limit, that is, when the (explicit) timescale separation parameter ε tends to 0; see also [START_REF] Desroches | Canards, folded nodes, and mixed-mode oscillations in piecewise-linear slow-fast systems[END_REF]. It turns out that this condition cannot be obtained in a slowly-forced integrator system if the forcing is too simple, that is, harmonic; see Section 2. As we will show in Section 3, one needs a feedback term from the voltage in the forcing equation in order to obtain the subthreshold oscillations, which suggests in the context of real neurons, an autaptic behaviour, and can be tested experimentally using a dynamic-clamp protocol.

We will showcase our strategy with a simple biophysical example of 2D type-I neuron model, namely the I Na /I K model proposed by Izhikevich in [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF]; however our approach will work with any type-I neuron model. Noteworthy, it does not require to alter the underlying bifurcation structure of the model, which is customary in studies reporting a switch from integrator to resonator. [START_REF] Prescott | Pyramidal neurons switch from integrators in vitro to resonators under in vivo-like conditions[END_REF][START_REF] Yi | Exploring how extracellular electric field modulates neuron activity through dynamical analysis of a two-compartment neuron model[END_REF][START_REF] Guantes | Variability in noise-driven integrator neurons[END_REF] The excitability switch that we propose here is purely due to timescale separation between the model and the forcing. We also showcase the reverse scenario, namely a switch from resonator to integrator. This requires to have a folded node instead of a folded saddle, therefore a different slow forcing structure, however here again we obtain the switch from one neuronal type to the other by staying within the same bifurcation scenario, only playing with the slow-fast structure of the model and of its folded singularity. The folded-node scenario indeed induces a particular geometry for the trajectories passing near such a folded singularity. Namely, they make transient small-amplitude oscillations, which correspond to subthreshold oscillations in the neuronal context [START_REF] Desroches | Mixed-Mode Oscillations with Multiple Time Scales[END_REF][START_REF] Wechselberger | Canard Theory and Excitability[END_REF] and their number can be controlled, e.g., by varying initial conditions. Indeed, families of initial conditions giving rise to the same number of oscillations form so-called rotation sectors in phase space. Furthermore, by controlling the trajectory to flow into the first rotation sector, one can suppress these subthreshold oscillations and hence, turn the behaviour into an integrator, hence obtaining the reverse switch. However, we will argue that the folded-saddle scenario is more appropriate for this switch between integrator and resonator neuron in order to obtain a behaviour as close as possible to experiments.

This article is organised as follows. In Section 2, we present the I N a /I K model and analyse numerically its integrator structure. Then, we apply to it a first slow harmonic forcing and show that it is sufficient to create a folded-saddle singularity but insufficient to obtain subthreshold oscillations, which require a more elaborate forcing. This is why, Section 3, we adapt the forced current in order to obtain the singular-limit algebraic condition giving rise to subthreshold oscillations in the full model, and hence the resonator behaviour. In Section 4, we present the reverse scenario, whereby a resonator neuron can behave like an integrator, and show that this is due to another type of folded singularity, namely a folded-node singularity. Finally, we conclude in Section 5 and propose a strategy to verify experimentally our theoretical predictions. [START_REF] Amir | Burst discharge in primary sensory neurons: Triggered by subthreshold oscillations, maintained by depolarizing afterpotentials[END_REF] Izhikevich's I Na /I K model.

Under constant external current

We consider a two-dimensional conductance-based neuron model with minimal components for excitability. This model was proposed by Izhikevich in [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF] and it was referred to as the I Na /I K model, since it only assumes basic persistent sodium with instantaneous activation, potassium and leak currents; this is the name we shall use throughout this article. The model's equations are:

CV = I -g L (V -E L ) -g N a m ∞ (V )(V -E N a ) -g K n(V -E K ), n = n ∞ (V ) -n τ n (V ) , (1)
with steady-state functions:

x ∞ (V ) = (1 + exp((V x,1/2 -V )/k x )) -1 , x = {m, n}.
For simplicity, we take the time constant τ n (V ) to be independent of V and we shall fix its value to 1. System (1) is based on a simplification of the two-dimensional reduction of the Hodgkin-Huxley model proposed by Krinskii & Kokoz [START_REF] Vi | Analysis of the equations of excitable membranes. i. reduction of the hodgkins-huxley equations to a 2d order system[END_REF], independently by Rinzel [START_REF] Rinzel | On repetitive activity in nerve[END_REF], and further studied, e.g., by Moehlis [START_REF] Moehlis | Canards for a reduction of the hodgkin-huxley equations[END_REF]. Thus, the variable V represents the membrane potential of the neuron, and n, the activation of potassium channels. The constants g x Parameter values are:

C = 1, E L = -80, E Na = 60, E K = -90, g L = 8, g Na = 20, g K = 10, V m,1/2 = -20, K m = 15, V n,1/2 = -25, K n = 5, τ n (V ) = 1.
(where x corresponds to L, Na, or K) are the maximal conductances of the ionic currents considered, E x are the Nernst potentials of the ionic species and C is the capacitance of the neural membrane; I denotes an externally applied current.

The bifurcation diagram of system (1) with respect to I, shown in Fig. 1, is typical of a neuron with type-I excitability [START_REF] Ermentrout | Type I membranes, phase resetting curves, and synchrony[END_REF]. Namely, a family of low-voltage equilibria (rest states of the neuron) destabilise and give way to a family of stable limit cycles (spiking states of the neuron) via a SNIC bifurcation, which occurs at an input current value I ≈ 4.51. The SNIC bifurcation being a homoclinic-type bifurcation, the emerging stable cycle has a very large period (tending to infinity at the bifurcation), hence a very small frequency, which is a key hallmark of type-I excitability. At a much higher value of the input current, the branch of stable cycles disappears via a supercritical Hopf bifurcation at I ≈ 200. Therefore, system (1) is considered to be in an integrator regime here.

Applying a slow sinusoidal external current

We now consider a periodic forcing to system (1) in the form of a slow externally-applied sinusoidal current. This can be done by replacing the constant term I in the V -equation of ( 1) by a time-dependent function I(t) = I 0 + sin(εt), with ε > 0 a small constant. However, to further analyse the resulting periodically forced system using geometric singular perturbation theoretic tools [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], it is more appropriate to write it in autonomous form and obtain the slow sinusoidal forcing I(t) as the solution of a harmonic oscillator. Hence, we consider system (1) forced by the following slow differential equations:

I = -εJ, J = ε(I -I 0 ). (2)
We are therefore considering the following extended 4D system:

CV = I -g L (V -E L ) -g Na m ∞ (V )(V -E Na ) -g K n(V -E K ), n = n ∞ (V ) -n τ n (V ) , I = -εJ, J = ε(I -I 0 ), (3)
where the prime denotes differentiation with respect to the fast time τ . Hence, ( 3) is a slow-fast dynamical system with two fast variables V and n, and two slow variables I and J. As customary in multiple-timescale dynamics, we can rescale time by a factor ε and introduce the slow time t = ετ , which brings the system in a different time parametrisation that will be helpful when styding its slow singular (ε = 0) limit, namely:

εC V = I -g L (V -E L ) -g Na m ∞ (V )(V -E Na ) -g K n(V -E K ), ε ṅ = n ∞ (V ) -n τ n (V ) , İ = -J, J = I -I 0 , (4)
where the over-dot denotes differentiation with respect to t.

As long as ε = 0, the systems ( 3) and ( 4) are equivalent, they have the same phase portraits, but the solution trajectories are parameterized differently. Furthermore, their respective singular limits are different and highlight different aspects of the original system's dynamics: the fast components in the case of system (3), and the slow components for (4).

The fast singular limit (i.e., system (3) with ε = 0) corresponds to the original integrator I Na /I K model [START_REF] Al | M-current induced bogdanov-takens bifurcation and switching of neuron excitability class[END_REF], which is logical since the extended system (3) was obtained by slowly forcing this integrator model. Recall that its bifurcation structure, shown on Fig. 1, is characterised by the presence of a SNIC bifurcation. The S-shaped curve of equilibria of this fast subsystem (1) is called critical manifold of the full system and we label it as S 0 in Fig. 1 and subsequent figures. Its algebraic expression is given by

S 0 := I -g L (V -E L ) -g Na m ∞ (V )(V -E Na ) -g K n ∞ (V )(V -E K ) = 0 (5)
Hence, the critical manifold, which in the present case is a surface in R 4 (since it has the dimension of the slow variables), can be seen as a graph over V and its equation can be written 2 shows a trajectory of the slowly forced system (3) superimposed onto the critical manifold S 0 and displayed only in the vicinity of the fast subsystem SNIC bifurcation point, for specific initial conditions of the forcing I(0) = I 0 and J(0). The trajectory (in black, with arrows representing the direction of motion) follows the stable branch of equilibria of the fast subsystem (red branch of the parabolicshaped bifurcation curve), which is the behaviour predicted by slow-fast theory. However, near the SNIC point -also labelled FS for reasons related to the slow subsystem (see below) -the trajectory turns around the point, instead of being repelled away, before flowing backwards as the forcing changes direction. This behaviour is counter-intuitive, it has to do with certain types of canards, and it is best explained by considering the other singular limit, namely the slow limit, of the forced integrator system (3), as we do next. Similar to the fast subsystem, we take the ε = 0 limit of system ( 4), which provides a good approximation of the slow dynamics of the forced integrator system. However, the slow singular limit is very different from the fast one and it is given by the following

S 0 = {I = f (V )}. Figure
0 = I -g L (V -E L ) -g Na m ∞ (V )(V -E Na ) -g K n(V -E K ), 0 = n ∞ (V ) -n τ n (V ) , İ = -J, J = I -I 0 . (6)
The first two equations of (6) express the fact that in the slow singular limit, the dynamics is constrained to evolve only on S 0 . In passing, this shows the importance of the critical manifold in both singular limits. The other two equations are the slow differential equations written in the slow-time parametrisation. The resulting system ( 6) is complicated to study as such, in particular because of the algebraic constraint which hides the limiting fast dynamics. However, one can rescue it by differentiating this algebraic constraint with respect to time which, after projecting the dynamics onto the (V, J) plane (because the slow singular dynamics is essentially 2D) and rearranging terms, yields the following version of the slow subsystem:

f V (V ) V = -J, J = I -I 0 , (7)
with I = f (V ) due to the constraint to evolve on S 0 and where f V denotes the derivative of f with respect to V , namely,

f V (V ) = g L + g Na (m ∞,V (V )(V -E Na ) + m ∞ (V )) + . . . . . . g K (n ∞,V (V )(V -E K ) + n ∞ (V )), (8)
with:

x ∞,V (V ) = dx ∞ (V ) dV = exp (V x,1/2 -V )/k x k x x ∞ (V ) 2 .
System [START_REF] Desroches | Spike-adding in a canonical three-time-scale model: superslow explosion and folded-saddle canards[END_REF] is the more practical form of the slow subsystem or reduced system (RS). This limiting system is singular along the zero set of f V , which geometrically corresponds 9), B: RS system (1). In both panels, the orange lines are the components of the V -nullcline, while the green curve is the J-nullcline. A: the DRS has 3 equilibria, 2 centers (blue dots) and 1 saddle (red dot); the red curves are the stable and unstable manifolds of the saddle, they coincide and form a double homoclinic connection, each loop surrounding one of the centers. B: there are 2 center equilibria (blue dots) and one folded saddle (red dot). The red curves are the singular true and singular faux canards, they coincide to form a double folded homoclinic connection [START_REF] Desroches | Spike-adding in parabolic bursters: The role of folded-saddle canards[END_REF][START_REF] Desroches | Spike-adding in a canonical three-time-scale model: superslow explosion and folded-saddle canards[END_REF].

Parameter values are as in Fig. 1.

to the fold set F := {f V (V ) = 0} of the critical manifold. The critical manifold of system ( 4) is a cubic surface and its fold set F has two connected components; see, e.g., Fig. 4 for an illustration of the lower fold curve F of S 0 . Noteworthy, the fold curve locally separates each of the two attracting sheets of S 0 , along which f V (V ) > 0, from the repelling sheet defined by f V (V ) < 0. The attractiveness and repulsiveness of S 0 is inherited from the stability and instability of equilibria of the fast subsystem (1), given that its set of equilibria precisely corresponds to S 0 . For a similar reason, the fold set F corresponds to the set of saddle-node bifurcations of the fast subsystem, hence the lower fold curve F corresponds to a family of SNIC bifurcation points of the fast subsystem. Namely, the SNIC point shown in Fig. 1 and detected when varying parameter I in (1) does not depend on J, hence we obtain a line of such point in the forced system (4).

In order to understand the flow of the RS near the fold curve F , one classical approach is to desingularise system (7) by rescaling time by a factor f V (V ), which brings forth the so-called desingularized reduced system (DRS)

V = -J, J = f V (V )(f (V ) -I 0 ), (9)
where the prime denotes differentiation with respect to the new time, i.e., after desingularization. System ( 9) is now defined everywhere on R 2 including along the fold set F.

Two points are worth noting about the DRS. First, the change of time by a factor f V (V ) artificially creates in (9) the possibility for equilibria on the fold set F, in particular on the lower fold curve F . Such equilibria satisfy the algebraic conditions

J = 0, f V (V ) = 0,
whereas the other equilibria of (9) satisfy the algebraic conditions and these are also equilibria of the RS system [START_REF] Desroches | Spike-adding in a canonical three-time-scale model: superslow explosion and folded-saddle canards[END_REF]. Second, the same change of time, because of the V -dependent factor, makes the orientation of trajectories of the DRS be opposite to that of the RS whenever f V (V ) < 0, that is, along the repelling (middle) sheet of S 0 . Hence, the DRS is a standard planar dynamical system, which can for instance have a saddle equilibrium on F , and the RS has the same geometrical orbits but the reversal of orientation along the repelling sheet of S 0 implies that the saddle equilibrium of the DRS is not an equilibrium anymore in the RS. Rather, it is a special point called folded saddle, which two special trajectories reach in finite time and cross. One trajectory crosses it from the attracting side of S 0 upwards, continuing along the repelling side, and it corresponds to the stable manifold of the saddle equilibrium of the DRS. The other trajectory crosses it in the opposite direction and it corresponds to the unstable manifold of the saddle of the DRS. Both trajectories are related to canards in that they cross from one side to the other of the critical manifold via a folded singularity. In the folded-saddle case, they are called true singular canard and faux singular canard, respectively [START_REF] Desroches | Mixed-Mode Oscillations with Multiple Time Scales[END_REF].

J = 0, f (V ) = I 0 ,
Figure 3 shows the phase portrait of the DRS of system (4) on panel A, and of the corresponding RS's phase portrait on panel B. A saddle equilibrium of the DRS is located on F , therefore it corresponds to a folded saddle in the RS. The DRS has two other equilibria, both of center type, they are still (true) equilibria of the RS and, hence, they will also influence the dynamics of the full system. Figure 4 shows the RS's phase portrait in a 3D (I, J, V ) projection, where the critical manifold S 0 is indeed a surface and F a curve; in fact, F is a straight line here since it does not depend on J. The figure illustrates well the geometry of such problems and the role of folded-saddle singularities in shaping the dynamics of slowly periodically-forced type-I neuron models. Indeed, such systems effectively correspond to parabolic bursters and folded saddles organise the appearance of spikes in such bursters along solution branches in parameter space; see [START_REF] Desroches | Spike-adding in parabolic bursters: The role of folded-saddle canards[END_REF] for details.

On Fig. 4, a spiking solution of system ( 4) is shown on top of S 0 and it clearly appears that, as it reaches the (lower) fold curve F of S 0 and comes close to the folded-saddle point FS, the trajectory follows the true singular canard, then it makes a spike and, as the voltage is going down back to baseline, the trajectory follows the faux singular canard. Therefore in this context of type-I membrane model with slow periodic forcing, the spike-adding threshold is organised by folded-saddle canards [START_REF] Desroches | Spike-adding in parabolic bursters: The role of folded-saddle canards[END_REF]. Therefore, with a slow harmonic forcing, an integrator neuron like system (1) displays a folded-saddle singularity and still behaves as an integrator. In particular, it cannot generate subthreshold oscillations. This is essentially due to the eigenvalue ratio of the saddle equilibrium of the DRS, as we will see next.

3 Integrator neuron with resonator behaviour

The eigenvalue ratio of the DRS's saddle equilibrium

We have just seen that the slowly forced integrator neuron model has trajectories with no subthreshold oscillations. Our aim is now to create these subthreshold oscillations in order to effectively obtain a switch from integrator to resonator behaviour.

Recent results by Mitry and Wechselberger [START_REF] Mitry | Folded saddles and faux canards[END_REF], also confirmed in the context of piecewise-linear slow-fast systems in [START_REF] Desroches | Canards, folded nodes, and mixed-mode oscillations in piecewise-linear slow-fast systems[END_REF], show that it is possible to obtain subthreshold oscillations near a folded saddle. More precisely, near the faux canard a folded saddle, which is the perturbation for ε > 0 small enough of the singular faux canard described in the previous section and shown in Fig. 4. As proven in [START_REF] Mitry | Folded saddles and faux canards[END_REF], the algebraic condition to obtain these subthreshold small-amplitude oscillations is on the ratio µ of the eigenvalues of the saddle equilibrium of the DRS, the ratio being of the unstable eigenvalue over the stable one. Necessarily µ is negative, however if it is strictly contained between -1 and 0, then such oscillations appear around the folded saddle's faux canard; see already Fig. 5 for an illustration.

As explained in the previous section, spiking trajectories of the forced system (4) follow the singular faux canard -hence, they also stay close to its ε-perturbation, the faux canard -right after making a spike, as the voltage goes down towards baseline. Hence, provided we can obtain subthreshold oscillations near the faux canard, then these will adequately resemble those obtained in resonator neurons, which will provide us with our objective of turning an integrator neuron into a resonator one purely based on a slow-fast effect. Now, a rapid glance at the DRS (9) makes us conclude that its Jacobian matrix has zero trace, whatever the equilibrium solution around which one linearises. Hence, the saddle equilibrium of the DRS, which corresponds to the folded-saddle of system (4), is a neutral saddle, implying that the ratio of its eigenvalues is necessarily equal to -1. This will always be the case when slowly forcing an integrator model if the forcing is harmonic. Therefore, to obtain subthreshold oscillations one needs to consider a more elaborate forcing, namely one that includes a feedback term from the voltage.

Slow forcing with feedback from the voltage

We now consider a more general forcing, with a feedback term in V in the I equation in order to obtain a non-zero trace in the Jacobian matrix of the new DRS evaluated at the saddle equilibrium of interest. We will keep the J equation only dependent on I as we simply need one non-zero diagonal element in the Jacobian matrix evaluated at this saddle equilibrium in order to ensure that its eigenvalue ratio will be different than -1. For simplicity, we will keep the dependence in V in the I equation linear and show that it suffices to obtain the expected behaviour both at the level of the eigenvalue ratio and in the full system's solutions. Specifically, we define the new slow forcing (written in fast time) as

I = ε (-J + αV ) , J = ε (I -I 0 ), (10)
which then yields the new DRS (after rescaling to the slow time)

V = -J + αV, J = f V (V )(f (V ) -I 0 ). (11)
Note that the voltage equation of the DRS is obtained from the slow differential equation for I in the forcing system, and that I must be kept equal to f (V ) in the slow singular limit. For these reasons, it would suffice that the I equation of the forcing depends on I and not on V , on top of its J dependence, in order to obtain a non-zero trace in the Jacobian matrix of the DRS corresponding to this new forcing. This would amount to replacing V by I in the first equation of [START_REF] Fellous | Frequency dependence of spike timing reliability in cortical pyramidal cells and interneurons[END_REF]. As a consequence, this in particular would avoid to use a feedback term in V in the full system [START_REF] Feldmann | All-optical spiking neurosynaptic networks with self-learning capabilities[END_REF]. However, it turns out that the full dynamics would not be able to exploit within its spiking regime the eigenvalue ratio of the slow singular limit, as additional unwanted equilibria would arise.

Therefore we keep the new slow forcing system [START_REF] Feldmann | All-optical spiking neurosynaptic networks with self-learning capabilities[END_REF] to obtain resonator behaviour. This new forcing may appear for now as the result of some ad hoc reverse engineering process, however we shall propose in the discussion section an explanation for its form and a possible experimental implementation of it.

We can now verify that the new DRS [START_REF] Fellous | Frequency dependence of spike timing reliability in cortical pyramidal cells and interneurons[END_REF] does possess a saddle equilibrium on F and that one can take a value of parameter α so that its eigenvalue ratio is strictly between -1 and 0. At an equilibrium (V * , J * ) located on the lower fold curve F of S 0 , the Jacobian matrix of (11) reads

J = α -1 f V V (V * )(f (V * ) -I 0 )) 0 ,
where f V V (V ) is the second derivative of f with respect to V . Given that α contributes to the trace of J and not to its determinant, it is clear that we still have a saddle equilibrium on F , and hence a folded saddle in the full system with the new forcing. Therefore, the eigenvalue ratio µ is given by

µ = α + α 2 -4f V V (V * )(f (V * ) -I 0 ) α -α 2 -4f V V (V * )(f (V * ) -I 0 ) . (12)
We then verify numerically that, for α negative and sufficiently large in absolute value, µ is indeed strictly between -1 and 0. For instance, by fixing α = -4 one can observe the expect subthreshold oscillations when simulating the full system, as illustrated in Fig. 5 where we also show the critical manifold S 0 (blue surface), its lower fold curve F (dotted line), the folded-saddle singularity, labelled FS (red dot), and the two singular canards (red curves traced on S 0 ). Any values of α such that µ in expression ( 12) is between -1 and 0 will work equally well. The plotted trajectory is entirely subthreshold, however it comes close to the folded saddle and then turns back; as the voltage is going down towards baseline, it oscillates around the singular faux canard. Hence we have obtained a key feature of a resonator neuron by simply modifying the slow forcing received by an a priori integrator neuron without modifying its bifurcation structure obtained with constant forcing. Resonator models have also other features which we can as well recover here.

In Figure 6, we highlight this resonator effect in the time series of the full system with the new forcing, obtained by taking an ensemble of initial conditions for the forcing, namely varying I 0 . What we observe is a transition in the voltage response from no spike, to one spike and then two spikes, depending on the value of I 0 . Every trajectory with at least one spike has clear subthreshold oscillations after the spike, or after the second spike for trajectories that have two spikes. What is more, the values of I 0 at which we observe a first spike in the voltage response, and then a second spike, are quite specific. This confirms that the forced system resonates with specific inputs. Consequently, we have obtained the main features of a resonator neuron. Again, the main mechanism behind this switch of behaviour is purely due to the slow forcing received by the integrator model, and based upon multiple-timescale dynamical phenomena.

4 A multiple-timescale scenario for the reverse switch: from resonator to integrator So far, we have mostly focused on the switch from an integrator to a resonator while retaining the characteristics of the integrator in absence of forcing. It is also possible to obtain the reverse switch, that is, to have a resonator neuron that behaves like an integrator. To do so, we keep system (1) but now consider yet another slow forcing, namely:

I = ε (-βJ + α(V -V 0 )) , J = ε (I -I 0 ), ( 13
)
where β is a new parameter that allows us to regulate the applied current; β was implicitly equal to 1 in the previous slow forcing [START_REF] Feldmann | All-optical spiking neurosynaptic networks with self-learning capabilities[END_REF]. We also add V 0 as new parameter as it helps controlling the amplitude of the feedback term from the voltage; V 0 was implicitly equal to 0 in the previous slow forcing. The new parameter β also enables to change the type of folded singularity by changing the topological type of the DRS equilibrium located on the fold curve F . In particular, varying β may turn the saddle into a node, hence giving a folded node in the full system. The non-singular (ε > 0 small) dynamics near a folded node is well known to produce small oscillations, which in the neuronal context correspond to subthreshold oscillations [START_REF] Desroches | Mixed-Mode Oscillations with Multiple Time Scales[END_REF]. Also known is the fact that, near a folded node, the phase space is locally partitioned into rotation sectors in which trajectories make a fixed number of subthreshold oscillations. Hence, it suffices to take initial conditions into the first rotation sector in order to obtain that the dynamics of this resonator system appears to behave like an integrator; Fig. 7 illustrates this effect. Therefore, it is interesting to showcase this effect of a resonator that behaves like an integrator without changing anything to its structure. However, the observed behaviour has a notable discrepancy with standard resonator neurons, whereby the subthreshold oscillations occur before the spike and not after; see Fig. 7 A. Yet, this type of scenario can be related to experimental observations since there are neurons with this kind of electrical behaviour recorded in vitro; see [START_REF] Amir | Burst discharge in primary sensory neurons: Triggered by subthreshold oscillations, maintained by depolarizing afterpotentials[END_REF] for an example in dorsal root ganglion neurons of rats.

Conclusion

In this work, we have studied a novel mechanism for excitability switch, from integrator to resonator and vice versa, within the framework of multiple-timescale dynamical systems. As a proof-of-concept, we considered the simple yet biophysical 2D example of the I N a /I K model by Izhikevich, but it would work just as well in any type-I model. Starting from a parameter set in which the model behaves as an integrator in the classical sense, that is, where the spiking regime occurs through a SNIC bifurcation, we first applied a slow periodic forcing in order to create a multiple-timescale structure in the forced model, and to show that the resulting 4D model possesses a folded-saddle singularity. Following standard geometric singular perturbation theory (GSPT), we derived the desingularized reduced system (DRS) and showed that it had a saddle equilibrium on the fold curve of the critical manifold of the original system. This slow periodic forcing preserved the integrator behaviour. However, modifying the forcing by including a feedback term from the membrane potential, we managed to affect the eigenvalues of the saddle equilibrium of the DRS. This, according to recent results from GSPT [START_REF] Mitry | Folded saddles and faux canards[END_REF], has the effect to allow for small-amplitude oscillations in trajectories that pass near the folded saddle. In turn, this provided us with a mechanism to obtain subthreshold oscillations in a neuron model that was a priori an integrator, hence making it behave like a resonator. The feedback term in V in the slow current forcing can be seen as a simple form of autaptic connection, which is physiologically plausible.

Therefore, we can see these results as theoretical predictions that we would like to verify experimentally. From the standpoint of electrophysiological measurements on real neurons, the possibility to apply a current that depends in real time on the readout potential from the neuron can be obtained through a dynamic-clamp protocol [START_REF] Ori | Dynamic clamp constructed phase diagram for the hodgkin and huxley model of excitability[END_REF][START_REF] Sharp | Dynamic clamp: computergenerated conductances in real neurons[END_REF]. The experimental setup allows to inject currents to a cell that depends upon the measured voltage in real time. It can be used in conjunction with pharmacological blockade of e.g. one ionic channel of the cell and injecting back the corresponding current as a result of a computer simulation using the measured voltage. We are current working on validating our theoretical prediction using dynamic-clamp experiments.

From the modeling's point of view, adding this term αV in the slow differential equation of the forcing current was a way to modulate the determinant of the Jacobian matrix of the DRS, and therefore obtain subthreshold oscillations after a spike as well as modulate their number. Thus, it is possible to keep the bifurcation structure and the characteristics of an integrator neuron model while forcing it to display the specific behaviour of a resonator neuron.

This difference with the classical scenario of integrator neuron models can potentially have interesting fallouts in the study of information transmission between neurons. Indeed, by definition, an integrator neuron integrates the message and returns to the rest electrical potential without subthreshold oscillation. The fact that an integrator neuron can resonate like a resonator neuron shows that the transmission of information is more complex than a simple integration of a message given upstream. This is an interesting avenue for follow-up research on this topic.

We have also shown that the inverse scenario of switch from a resonator towards an integrator is also possible, by using a folded-node scenario in place of a folded-saddle one. However, the subthreshold oscillations obtained with a folded-node scenario occur before the spike rather than after, which is uncommon for resonator neurons. This may be associated with certain types of neurons (as reported in e.g. [START_REF] Amir | Burst discharge in primary sensory neurons: Triggered by subthreshold oscillations, maintained by depolarizing afterpotentials[END_REF]) and to the dynamical phenomenon of mixed-mode oscillations (MMOs) [START_REF] Desroches | Mixed-Mode Oscillations with Multiple Time Scales[END_REF]. The most important aspect of this numerical experiment is that the resonance phenomenon does occur and it can be related to known properties of folded-node canards, in response to specific forcing inputs [START_REF] Wechselberger | Existence and bifurcation of canards in \mathbbrˆ3 in the case of a folded node[END_REF]. Yet, one can control the trajectories in such a scenario so that no subthreshold oscillation occurs as the solution flows past the folded node, which effectively makes the resonator behave like an integrator. We plan to verify experimentally this theoretical and computational prediction, first in neuromorphic analog circuits and then, using dynamic clamp, in real neurons.

Figure 1 :

 1 Figure 1: Bifurcation diagram of the system (1) with respect to parameter I. Red (resp. black) segments of the S-shaped curve of equilibria denote stable (resp. unstable) branches. As the applied current I is increased, oscillations (spikes) appear through a SNIC bifurcation (red dot) and then disappear through a supercritical Hopf bifurcation (blue dot denoted HB). Parameter values are: C = 1, E L = -80, E Na = 60, E K = -90, g L = 8, g Na = 20, g K = 10, V m,1/2 = -20, K m = 15, V n,1/2 = -25, K n = 5, τ n (V ) = 1.

Figure 2 :

 2 Figure 2: Folded-saddle canard trajectory (black curve with arrows) from (3) superimposed on the fast subsystem's bifurcation diagram already shown in Fig. 1. The red dot corresponds to a SNIC bifurcation point of the fast subsystem, as well as the folded-saddle singularity (FS) found in the slow subsystem; see the end of section 2.2. Black curve: unstable stationary points, red curve: stable stationary points. Parameter values are as in Fig. 1 except: I 0 = 4, ε = 0.001.

Figure 3 :

 3 Figure3: Phase portraits in the (J, V ) plane. A: DRS system (9), B: RS system (1). In both panels, the orange lines are the components of the V -nullcline, while the green curve is the J-nullcline. A: the DRS has 3 equilibria, 2 centers (blue dots) and 1 saddle (red dot); the red curves are the stable and unstable manifolds of the saddle, they coincide and form a double homoclinic connection, each loop surrounding one of the centers. B: there are 2 center equilibria (blue dots) and one folded saddle (red dot). The red curves are the singular true and singular faux canards, they coincide to form a double folded homoclinic connection[START_REF] Desroches | Spike-adding in parabolic bursters: The role of folded-saddle canards[END_REF][START_REF] Desroches | Spike-adding in a canonical three-time-scale model: superslow explosion and folded-saddle canards[END_REF]. Parameter values are as in Fig.1.

Figure 4 :

 4 Figure 4: Phase portrait of system (3) projected onto the (I, J, V ) space. Also shown is the critical manifold S 0 (blue surface), the lower fold curve F (dotted line), the singular true and faux canards (red curves), the folded-saddle singularity (red dot), labelled FS, and a trajectory making one spike without subthreshold oscillations (black curve). Parameter values are as in Fig. 1 except: V (0) = -76, n(0) = 3 • 10 -05 , I(0) = I 0 = -5.48 and J(0) = 10.

Figure 5 :

 5 Figure 5: Phase portrait of system (4) projected onto the (I, J, V ) space. Also shown at the critical manifold S 0 (blue surface), the lower fold curve F (dotted line), the singular true and faux canards (red curves), the folded-saddle singularity (red dot), labelled FS, and a trajectory making subthreshold oscillations (black curve). Parameter values are as in Fig. 1 except for the new parameter α: α = -4.

Figure 6 :

 6 Figure 6: Integrator neuron (1) acting as a resonator with the slow forced current (10). Parameter values are as in Fig. 5 except : ε = 0.02, α = -8, V (0) = -113.11, n(0) = 3 • 10 -05 , J 0 = 546, and I(0) = I 0 ranging between -4 and -7.

Figure 7 :

 7 Figure 7: A: Resonator neuron model (1) with the slow forced current (13) for J(0) = -370. The folded node is labelled FN. B: Same model acting as an integrator for J(0) = -380, I(0) = -6.0. Parameter values are as in Fig. 5 except for the new parameter β: β = -1, and V (0) = -65.933.

AcknowledgmentsSR was supported by Ikerbasque (The Basque Foundation for Science); GG and SR acknowledges support from the Basque Government through the BERC 2022-2025 program and by the Ministry of Science and Innovation: BCAM Severo Ochoa accreditation CEX2021-001142-S / MICIN / AEI / 10.13039/501100011033 and through project RTI2018-093860-B-C21 funded by (AEI/FEDER, UE) and acronym "MathNEURO". MD and SR acknowledge the support of Inria via the Associated Team "NeuroTransSF".


Dr Andrew Nichols French
email: afrench@uswcl.ars.ag.gov

PhD; M C Frederic Jacob

Martha C Anderson

W P Kustas

W Timmermans

A Gieske

Z Su

H Su

M F Mccabe

F Li

PhD; W P Anderson

PhD; W Timmermans A Kustas

; Z Gieske

Su Su

F Mccabe

; J Prueger

; N Brunsell

B Su

J Prueger

N Brunsell

Surface energy fluxes with the Advanced Spaceborne Thermal Emission and Reflection radiometer

Keywords: Thermal infrared, ASTER, high spatial resolution, multispectral TIR, surface energy balance modeling, spatial variability Thermal infrared, ASTER, high spatial resolution, multispectral TIR, surface energy balance modeling, spatial variability

HAL is

only satellite sensor collecting multispectral thermal-infrared images, a capability allowing unprecedented surface temperature estimation accuracy for a variety of surface cover types. Availability of ASTER data to study surface energy fluxes allows direct comparisons against ground measurements and facilitates detection of modeling limitations, both possible because of ASTER's higher spatial resolution. The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is a special remote-sensing tool intended to help estimate the surface energy balance, a critical attribute for monitoring land surface climate, hydrological processes and vegetation health. Using ASTER observations with spatial resolutions ranging from 15 to 90 m, in combination with a model simulating energy transfer, detailed estimates of surface energy fluxes may potentially be retrieved over large portions of the Earth's landmass. Knowledge of energy fluxes is important, providing meteorological boundary conditions for convection within the atmospheric boundary layer, and evapotranspiration (ET) estimates for hydrological studies, including vegetation transpiration estimates for crop growth assessment. The required flux retrieval accuracy varies by application, but is typically ∼50 W m -2 , as suggested by Seguin et al. (1999).

The essential idea for energy flux estimation of ET is recognition that water vapor mass transport can equivalently be represented as energy transport.

Due to the large amount of heat represented by the evaporation process, water mass transport can be accurately estimated from known energy transport.

Since energy transport of water vapor is directly connected to thermal and water vapor gradients, one can estimate ET primarily from knowledge of four quantities: surface temperature, vegetation density, near surface air temperature and near surface humidity. The first two are readily estimated from remote-sensing, while the last two are obtainable from ground observations or possibly from remote-sensing contextual data.

Estimation of land surface ET is of course not new, and in recent years accurate point-based observations have been used to monitor ET variability at hourly time scales. Examples include experiments at Monsoon '90 [START_REF] Moran | Use of ground-based remotely sensed data for surface energy balance evaluation of a semiarid rangeland[END_REF], FIFE [START_REF] Kanemasu | Surface fluxes measurements in FIFE: an overview[END_REF]Shuttleworth et al., 1989), HAPEX-SAHEL [START_REF] Goutorbe | HAPEX-Sahel: a large scale study of land-atmosphere interactions in the semi-arid tropics[END_REF], EFEDA [START_REF] Pelgrum | An inter-comparison of techniques to determine the area-averaged latent heat flux from individual in situ observations : a remote sensing approach using the European Field Experiment in a Desertification-threatened Area data[END_REF], BOREAS (Sellers et al., 1997), NOPEX [START_REF] Halldin | Boreal forests and climate[END_REF], JORNEX [START_REF] Havstad | Jornada experimental range: a unique arid land location for experiments to validate satellite systems[END_REF] and ReSeDA (Olioso et al., 2002b). However, it has also been recognized that spatially variability of ET is large and that even the most advanced ways of measuring ET with tower-mounted Eddy Covariance (E-C) systems are often not representative of landscape scale ET [START_REF] Moran | Use of ground-based remotely sensed data for surface energy balance evaluation of a semiarid rangeland[END_REF][START_REF] Pelgrum | An inter-comparison of techniques to determine the area-averaged latent heat flux from individual in situ observations : a remote sensing approach using the European Field Experiment in a Desertification-threatened Area data[END_REF]. Consequently, incorporation of remote-sensing data is required, since this is the only way to simultaneously observe surface properties over large surface areas.

ET monitoring with remote sensing has yet to be proven as a reliable, robust approach. Recent experiments show promising results, with systematic agreement between estimates from ground and aircraft measurements over several test areas [START_REF] French | Surface energy fluxes over El Reno, Oklahoma, using high-resolution remotely sensed data[END_REF][START_REF] Kustas | Using satellite remote sensing to extrapolate evapotranspiration estimates in time and space over a semiarid rangeland basin[END_REF][START_REF] Moran | Use of ground-based remotely sensed data for surface energy balance evaluation of a semiarid rangeland[END_REF][START_REF] Norman | A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature[END_REF]Kustas et al., 1996;[START_REF] Zhan | An intercomparison study on models of sensible heat flux over partial canopy surfaces with remotely sensed surface temperature[END_REF][START_REF] Anderson | A twosource time-integrated model for estimating surface fluxes from thermal infrared satellite observations[END_REF][START_REF] Norman | Surface flux estimation using radiometric temperature: A dual-temperature-difference method to minimize measurement errors[END_REF][START_REF] Hasager | Land cover, surface temperature and leaf area index maps from satellites used for the aggregation of momentum and temperature roughnesses[END_REF][START_REF] Wassenaar | Estimation of evapotranspiration over heterogeneous pixels[END_REF]Jacob et al., 2002a;Olioso et al., 2002a). On the other hand, previous criticisms (e.g., [START_REF] Hall | Satellite remote sensing of surface energy balance: sucess, failures and unresolved issues in FIFE[END_REF] suggest that thermal infrared (TIR) satellite observations, particularly, are not sufficiently accurate to constrain flux models. Although such critiques have been refuted [START_REF] Kustas | Reply to comments about the basic equations of dual-source vegetation-atmosphere transfer models[END_REF]Bastiaanssen et al., 1998a,b), further examples are clearly needed before remote sensing can be used operationally for this purpose.

Accurate spatially distributed estimates of surface energy fluxes require physically based energy flux models driven by accurate remote-sensing observations. Such ideal models, which do not depend on empirical calibration, are unrealizable since complete consideration of processes at the local scale requires large amounts of unavailable information. Some models (e.g., SiB2, Sellers et al., 1986) estimate ET from the Penman-Monteith equation [START_REF] Penman | Natural evaporation from open water, bare soil and grass[END_REF][START_REF] Monteith | Principles of Environmental Physics[END_REF], which considers both the radiative flux and the water mass flux, and can use remote sensing to derive net radiation (R n ). This widely used approach, however, does not adequately model heterogeneous surfaces or water-stressed vegetation. Other approaches, such as the simplified method of [START_REF] Jackson | Wheat canopy temperature: a practical tool for evaluating water requirements[END_REF], can model some of these conditions by incorporating surface temperatures observed from thermal remote sensing, but do so in an empirical way, meaning that local calibration is required.

Recent versions of remote-sensing energy flux models have progressed to include treatment for a wide range of conditions including heterogeneous surfaces and variable meteorological constraints, thereby reducing problems encountered previously. These newer models are similar in some respects to pre-decessors such as Penman-Monteith models, but represent an advance to the state of the art because they incorporate higher spectral and spatial resolution remote-sensing data across the visible, near-infrared and thermal wavelengths.

Five examples are: the Two Source Energy Balance model (TSEB, Norman et al., 1995), DisAlexi [START_REF] Norman | Remote sensing of surface energy fluxes at 10 1 -m pixel resolutions[END_REF], the Surface Energy Balance Algorithm for Land model (SEBAL, Bastiaanssen et al., 1998a), the Surface Energy Balance System (SEBS, Su, 2002), and the NDVI/Temperature triangle method [START_REF] Gillies | A verification of the triangle method for obstaining surface soil water content and energy fluxes from remote measurements of the Normalized Difference Vegetation Index NDVI and surface radiant temperature[END_REF]. Although these models rely upon different assumptions and interpret remote-sensing data in different ways, they have all been tested and verified under local conditions. This means that any one of them ought to return accurate surface flux estimates given accurate input observations. However the question remains whether or not these models are generally applicable and can do what previous models could not, namely reliably estimate spatially distributed surface energy fluxes.

One way to answer this question is to assemble multispectral remote-sensing observations of high quality and high spatial resolution over various landscapes and to compare model results with each other and with independently measured ground flux measurements. In this respect, ASTER should be wellsuited to the task since it can accurately retrieve vegetation density, spectral reflectances over visible and near-infrared (VNIR) wavelengths (Jacob et al., 2002b), and surface temperatures potentially accurate to ∼0.5 • C [START_REF] Hook | Land surface temperature measured by ASTER-first results[END_REF]. In recent work by [START_REF] Jacob | Comparison of land surface emissivity and radiometric temperature from MODIS and ASTER sensors[END_REF], land surface temperature retrieval agreement between ASTER and MODIS sensors was ∼0.9 • C , indicating consistent data quality.

In this paper we demonstrate estimation of surface energy fluxes using ASTER data by examining first results from a model intercomparison study based on observations over a 2002 experimental study in Iowa, USA. The approach is to first discuss basic observational requirements for landscape-scale surface energy balance modeling and how ASTER data meet those requirements. Next we briefly describe implementations of the TSEB and SEBAL models. Third, these flux models are compared with each other and ground-level measurements.

Observational Requirements

Regardless of spatial scale, estimation of surface energy fluxes [START_REF] Brutsaert | Evaporation into the atmosphere: theory, history and applications[END_REF][START_REF] Mcnaughton | An evaluation of the Priestley-Taylor equation and the complementary relationship using results from a mixedlayer model of the convective boundary layer[END_REF]) can be represented by the balance of turbulent fluxes against available energy:

H + LE = R n -G ( 1 )
where H, sensible heat and LE, latent heat, are turbulence terms, G is conducted soil heat and R n is net radiation. These terms, all in W m -2 , respectively represent energy conducted away from the surface and energy delivered to the Earth's surface. Neglected are storage and photosynthesis terms. These latter two terms are usually thought unimportant for short-term estimates over non-forested areas. However, recent work by [START_REF] Meyers | An assessment of storage terms in the surface energy balance of maize and soybean[END_REF] suggests that storage in corn fields can be significant. Spatial implementation of Eq.1 requires an underlying energy balance model and boundary condition constraints from remote sensing data observations. Terms in Eq.1 constrained by remote sensing depend upon spatial and spectral resolutions. , 1996).

For this paper we examine ASTER's suitability by evaluating results from an experimental site in central Iowa known as the Soil Moisture Atmosphere

Coupling Experiment 2002 (SMACEX). The site consisted of intensive soil, vegetation and meteorological observations over the Walnut Creek Watershed, an area just south of Ames, Iowa, USA [START_REF] Kustas | Effects of remote sensing pixel resolution on modeled energy flux variability of croplands in Iowa[END_REF]. Two observational factors we consider are spatial and spectral resolution.

Spatial Resolution

For surface energy balance modeling the chief requirement for spatial resolution is the ability to distinguish land cover types that have distinctly different heat flux properties and resolve spatial variability of the considered processes within a land-class type. Over agricultural regions, these cover types include grazinglands, managed crops, riparian zones, water bodies and fields of bare soil. Each of these types differs significantly in surface roughness, potential heat capacity, moisture content, and spectral reflectance. Spatial variability of energy fluxes at local scales is strongly controlled by soil and micro-meteorological conditions.

ASTER images are potentially able to make these distinctions over all three band groups, with 15-m resolution for VNIR bands 1,2 and 3; 30-m resolution for SWIR bands 4-9; and 90-m resolution for TIR bands 10-14 [START_REF] Yamaguchi | Overview of Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER)[END_REF]. Considering the context of currently available images (e.g., Ikonos, Quickbird, TM7 and MODIS), ASTER data may be considered high-to moderate-resolution because they match TM7's capabilities in visible and near infrared wavelengths and have a significantly higher resolution than MODIS' best resolution (250 m).

An example of ASTER's spatial capability can be demonstrated using SMACEX imagery from 1 July 2002. Using 15-m NDVI images (Fig. 1), land-use patterns and highways are clearly defined at the landscape scale. NDVI was derived from Level 2 ASTER processing of band 2 and 3N reflectances [START_REF] Abrams | The Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER): Data products for the high spatial resolution imager on NASA's Terra platform[END_REF]).

In the Iowa example above, the adequacy of remote-sensing resolution was determined qualitatively, where imagery was checked against expected land cover patterns at the section and quarter section levels. Knowing that adequate resolution must be about half the distance of the dominant land-surface scale to avoid dominance of mixed pixels (e.g., [START_REF] Woodcock | The factor of scale in remote sensing[END_REF], minimal resolution over SMACEX is ∼1/8 mile (400 m). A more quantitative way to evaluate land-surface scales and required resolution, is to generate histograms at successively coarser resolutions and assess where the distributional pattern changes. Using ASTER 15-m NDVI images, landscape scale is identified from a series of histograms generated from image reflectance aggregations. In the SMACEX NDVI example (Fig. 3), source red and near-infrared A critical spatial issue for energy balance modeling is the ability to distinguish land-use patterns in the TIR. Surface temperature, rather than vegetation density, is the most influential parameter for instantaneous energy fluxes. ASTER surface temperature data collected over SMACEX with a 90-m spatial resolution resolves dominant land-use patterns discerned with VNIR data (Fig. 2). Surface temperatures, derived from the Temperature Emissivity Separation (TES) algorithm [START_REF] Gillespie | A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images[END_REF] and atmospherically corrected with radiosonde data and MODTRAN [START_REF] Berk | MODTRAN cloud and multiple scattering upgrade with application to AVIRIS[END_REF], range between 25 and 40 (Fig. 4).

ASTER's multispectral TIR imaging capability is its most distinguishing characteristic and is important for retrieving surface temperatures, potentially accurate to 0.5 • C . Reasons for this potential are its low NE∆T (noise-equivalent change in temperature, a measure of the signal-to-noise ratio for thermal detectors), <0.3 K, and its multiband placement within the TIR window Five TIR bands allow better estimation of spectral variability of surface brightness temperatures and surface emissivities [START_REF] Gillespie | A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images[END_REF]

than otherwise 3 Energy Flux Modeling

To show how ASTER can be used to model the surface energy fluxes, we combine ASTER imagery over the SMACEX site with implementations of two well-known energy balance models, the Two Source Energy Balance approach (TSEB, [START_REF] Norman | A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature[END_REF], and the Surface Energy Balance Algorithm for Land approach (SEBAL Bastiaanssen et al., 1998a). By showing results from two models, rather than one, model performance can be seen for the same landscape, and allows identification of model strengths and weaknesses.

TSEB

TSEB is a remote-sensing two-source approach based upon land surface separability into two distinct, but linked components: soil surfaces and vegetation canopies. By combining remote observations of surface temperature, vegetation data, land cover class, and near-surface meteorological data, surface resistance characteristics can be modeled in a more physically meaningful way (Kustas andNorman, 2000, 1999) than a more typical one-source approach based on the kB -1 parameter [START_REF] Garratt | Momentum, heat and water vapour transfer to and from natural and artificial surfaces[END_REF]. TSEB derives energy flux estimates from modeling the land as a resistance network, between energy sources from soil, vegetation, and the overlying atmosphere. Two main variants of TSEB exist, one applicable strictly at local scales (as described in [START_REF] Norman | A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature[END_REF], while the other, known as DisAlexi [START_REF] Anderson | A twosource time-integrated model for estimating surface fluxes from thermal infrared satellite observations[END_REF][START_REF] Mecikalski | Estimating fluxes on continental scales using remotely sensed data in an atmospheric-land exchange model[END_REF] is also useful at regional scales since it also models energy exchange at the atmospheric boundary layer. The model outputs used here are from the DisAlexi implementation of TSEB, but we note that results from the local implementation are similar due to the local scale of SMACEX. R n , computed from a multiple scattering model [START_REF] Monteith | Principles of Environmental Physics[END_REF][START_REF] Goudriaan | Crop micrometeorology: a simulation study[END_REF] is partitioned for each of the considered spectral ranges: visible, near-infrared and TIR wavelengths [START_REF] Norman | A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature[END_REF]. This partitioning allows discrimination between energy available at the soil surface and energy delivered to the vegetation canopy. TSEB has three key assumptions: turbulent fluxes are constant within the near surface layer (Monin-Obukhov similarity is used for stability correction), radiometric temperature can be repartitioned into soil and vegetation components, and Priestley-Taylor transpiration [START_REF] Priestley | On the assessment of surface heat flux and evaporation using large-scale parameters[END_REF] applies for unstressed vegetation. More complete details are described in [START_REF] Norman | A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature[END_REF] and in [START_REF] French | Surface energy fluxes over El Reno, Oklahoma, using high-resolution remotely sensed data[END_REF].

SEBAL

The SEBAL model, in contrast to TSEB, is a one-source modeling approach which does not discriminate soil and vegetation components. This model is based upon the estimation of surface energy fluxes with minimal ancillary information by considering spatial variability within the images themselves.

This design principal has great practical value since instantaneous evapotranspiration estimates can potentially be retrieved worldwide without the need to collect ground-based data. SEBAL does this by contextual separation of the land-surface images into 'wet' and 'dry' areas. Distinction is achieved using the 'evaporative controlled' and 'radiative controlled' branches of a surface temperature-versus-albedo diagram. If both areas are included within the images, the diagram commonly shows a characteristic parabolic pattern. The vertex marks the maximum observed radiometric temperature and separates the lower albedo wet areas from the higher-albedo dry areas.

An aggregation scheme of the dry area properties is used along with air temperature at an estimated blending height in order to compute a regional layer aerodynamic resistance, and then wind speed at a standard reference level

(2 m). Blending height is the height above the ground where the influence of local-scale surface heterogeneity upon atmospheric turbulence is relatively unimportant (e.g., [START_REF] Mahrt | Surface heterogeneity and vertical structure of the boundary layer[END_REF]. Air temperature is assumed to be a linear function of surface temperature. The slope and offset of this function are derived considering both wet and dry areas. Over the driest areas, LE flux is zero, allowing retrieval of H as a residual of R n -G. Over the wettest areas, H is negligible, with surface temperature equal to near surface air temperature.

Key assumptions for SEBAL are the simultaneous presence of wet and dry areas within the remote-sensing scene and the ability to parameterize landsurface resistance with a factor called kB -1 [START_REF] Garratt | Momentum, heat and water vapour transfer to and from natural and artificial surfaces[END_REF]. This factor is a means for obtaining thermal heat roughness from momentum roughness. Further details of SEBAL can be found in Bastiaanssen et al. (1998a), Su et al. (1999),and Jacob et al. (2002a).

Model Validation Results

To are from E-C flux observations, nominally 2 m above the ground in soybean fields, and 4 m above the ground in corn fields [START_REF] Kustas | Effects of remote sensing pixel resolution on modeled energy flux variability of croplands in Iowa[END_REF]. In rough terms, E-C measurements are represented by no more than ∼4 ASTER pixels in the upwind direction. Energy budgets at all stations were closed by the Bowen ratio approach (Twine et al., 2000), then adjusted for wind speed, wind direction and fetch effects using a flux footprint model from [START_REF] Schuepp | Footprint prediction of scalar fluxes from analytical solutions of the diffusion equation[END_REF].

Comparisons

Discussion

In light of the excellent retrieval capabilities of ASTER and the reasonably good agreement between modeled estimates over the SMACEX study site we can say that ASTER is performing consistently. ASTER data provide one of the few avenues to evaluate and test remote sensing-based surface energy flux models at higher resolutions. Spatially varying surface temperatures closely correspond to observed variations in vegetation densities and surface fluxes.

Spatial resolution of ASTER resolves physically distinct land-cover types that in the Iowa study have a dominant length scale of 250 m. Without the 15-to 90-m resolution capabilities, flux modeling would have necessitated more empiricism, and therefore more calibration than currently used in the TSEB and SEBAL models. Furthermore, without the higher-resolution images, model comparisons against ground observation sites would not be very meaningful.

Energy flux model intercomparison shows functional relationships exist between the SEBAL and TSEB, indicating both are responding similarly to ASTER temperature, NDVI and albedo input data. However, in general terms, SEBAL predicted fluxes are more sensitive to input changes than are TSEB predicted fluxes. For example, over low H regions, SEBAL and TSEB commonly agree within 20 W m -2 , but over high H regions, SEBAL predicts values up to 200 W m -2 greater than TSEB values. Spatial representation of this sensitivity difference is also strongly expressed with LE estimates (Fig. 5).

Distributional patterns are similar for both models, but in the SEBAL case, estimated LE areal values are much less than TSEB estimates.

General results from the tabular data shows cross-comparison results seen previously. Agreement between R n ground observations, and model estimates from TSEB and SEBAL are acceptable with mean deviations of 31 W m -2

or less, indicating realistic model estimation of both short-and long-wave radiation components. We note, however, that poor agreement would have been surprising because modeled estimates had the benefit of calibrated spatial albedo data and measured incoming solar radiation. Agreement between G observations and modeled estimates was also good, with average deviation of 15 W m -2 . This good agreement indicates that fractional estimation of net radiation based upon NDVI-derived fractional cover works well for the SMACEX area at satellite overpass time.

Agreement for turbulent fluxes H and LE was more problematic to assess.

In some cases, (e.g., dense corn field sites 6 and 151), SEBAL LE model results agreed very closely with ground observations (< 10 W m -2 ), while TSEB results showed differences on the order of 50 W m -2 . But in other cases (e.g., sparse soybean field sites 3 and 23), TSEB model results were in much better agreement with E-C observations. TSEB differences for H estimates were typically within 35 W m -2 of E-C values, while SEBAL estimates were highly discrepant (up to 150 W m -2 ).

Because these flux estimate discrepancies arise from identical input data, their causes have to be explained by TSEB/SEBAL model differences and assumptions. One important model difference lies with model representation of land and near surface air temperatures. In TSEB, H was determined by the gradient between observed land surface radiometric temperatures and an observed air temperature. At SMACEX a constant near surface air temperature was used with TSEB, meaning that modeled H values were directly related to observed radiometric surface temperatures. Based on the 7 W m -2 average deviation between E-C and TSEB H fluxes, this gradient estimation approach seems accurate within E-C measurement error. By contrast, the SEBAL model, returned comparatively less accurate H estimates, with an average deviation of 89 W m -2 . SEBAL does not use near surface air temperatures because of its temperature normalization using wet and dry areas. Though known to work well elsewhere, the application of this approach did not work well at SMACEX due to an inability to fully distinguish wet and dry areas. As noted previously, SEBAL anticipates a parabolic relationship between surface temperature and albedo, thus distinguishing evaporative and radiative controlled domains. At SMACEX this parabolic relationship did not exist (Fig. 7), despite known existence of sparsely vegetated and riparian zones. This means successful retrieval by SEBAL of meteorological and resistance parameters is questionable because no distinction between minimum and maximum temperatures within the radiative and evaporative domains is possible for this 1 July scene.

Current validation results at SMACEX suggest more accurate flux estimation from TSEB than from SEBAL (Table 1). As noted, typical H estimates from TSEB were within 7 W m -2 of ground observations, while typical LE estimates from TSEB deviated significantly less from observations than did SEBAL estimates. Most of the disagreements between TSEB and SEBAL estimates occurred over sparsely vegetated sites, suggesting that the soil-vegetation differentiation accommodated by TSEB is a significant model benefit.

Nevertheless, some caution is required with these results. especially considering the preliminary nature of surface flux observations. Issues such as energy budget closure and flux footprint estimates need more thorough investigation.

Such investigation is underway, which will also consider flux estimates from the SEBS and the NDVI/Temperature triangle approaches.

  Surface energy flux retrieval from ASTER is demonstrated using data collected over an experimental site in central Iowa, USA, in the framework of the Soil Moisture Atmosphere Coupling Experiment (SMACEX). This experiment took place during the summer of 2002 in a study of heterogeneous agricultural croplands. Two different flux estimation approaches, designed to account for the spatial variability, are considered: the Two-Source Energy Balance model (TSEB) and the Surface Energy Balance Algorithm for Land model (SEBAL). ASTER data are shown to have spatial and spectral resolution sufficient to derive surface variables required as inputs for physically based energy balance modeling. Comparison of flux model 1 Introduction

  Using VNIR to TIR detectors, constraints on three of four terms are directly possible: H, R n and G. H is constrained by surface temperatures derived from TIR data, while both R n and G fluxes are constrained by both VNIR reflectances and emitted TIR radiances. The remaining component, LE can be constrained indirectly by residuals. The useful role of VNIR-TIR remote-sensing in energy flux modeling is therefore determined by a sensor's ability to accurately determine the three flux terms from observed reflected and emitted radiances. Without sufficient accuracy, errors from estimates of H, R n and G accumulate in LE flux estimates and could overwhelm results by 100's of W m -2 (Kustas and Norman

  reflectance data were aggregated by spatial averaging to 90 m (ASTER TIR nominal resolution) and 240 m (close to 250 m, the best MODIS nominal resolution), then combined to yield NDVI. Since NDVI is closely coupled to energy balance model results, knowing the minimum resolution needed to represent relatively homogeneous areas of vegetation, rather than as clumped mixtures is important. Histograms shown in Fig. 3 indicate that the minimum required resolution is somewhere between 90-m and 240-m resolutions due to a significant distributional change. At 15-to 90-m spatial resolutions, vegetation cover is bimodally distributed. At 240-m resolution, the bimodal distribution is mostly lost, meaning that distinctive sparse vegetation values are mixed with dense vegetation ones.

4

  Intercomparison Results at SMACEX 2002 Results from the ASTER-based model intercomparison between TSEB and SEBAL model estimates are summarized as cross-plot comparisons on a pixelby-pixel basis (Fig. 6). These plots help identify how model estimates relate to each other by flux component and represent frequency counts of the components R n (upper left), G (upper right), H (lower left) and LE (lower right).The plotted distributions represent spatially comparative SEBAL vs. TSEB fluxes for corn and soybean fields over an area approximately 13 km x 8 km over the Walnut Creek watershed. This subset corresponds to ∼13,000 ASTER pixels. The pixel frequency counts in Fig.6were binned according to range of observed fluxes with highest occurrence frequencies coded dark gray.In aggregate, R n estimates by SEBAL and TSEB show moderate systematic agreement (R 2 = 0.68), but show a bias of 13 W m -2 . For soil heat flux G, TSEB and SEBAL calculations are strongly correlated (R 2 = 0.95), moderate bias of 11 W m -2 , but have differing sensitivities. Where G values are less than ∼100 W m -2 , SEBAL and TSEB values agree within standard error of mean values of 8 W m -2 , but with G greater than 100 W m -2 , SEBAL estimates are on the order of 50 W m -2 greater than TSEB estimates. The relationships between turbulent fluxes H and LE are similarly systematic, but with greater differences between low flux values and high flux values. For H estimates, the two models showed good correlation (R 2 = 0.80), and bias of 33 W m -2 . However, SEBAL model sensitivity was greater than that of TSEB, where the crossplot slope is 1.46, yielding flux estimates differing by over 100 W m -2 . The LE comparison results complement the H results, also with good correlation (R 2 = 0.89), overall bias of 8 W m -2 , and greater SEBAL model sensitivity (crossplot slope of 1.25).

  distributed surface energy fluxes over agricultural areas requires high-quality remote sensing, a need mostly fulfilled by ASTER. With spatial resolution capabilities ranging between 15 and 90 m, ASTER can detect and discriminate variations in surface temperature, emissivity, vegetation densities and albedo corresponding to distinct land use types and thus reduce problems from mixed-pixels. In this study we have combined consistent ASTER observations with physically based surface energy flux models to retrieve reasonable estimates of instantaneous surface energy fluxes. Using model intercomparison results showed systematic agreement between all flux components, indicating that the two models tested, TSEB and SEBAL, operate similarly when provided identical remote-sensing inputs. These first results from the multi-model intercomparison study also showed significantly different model sensitivities, indicating a need for further, in-depth analysis of each of model. Such investigations are possible because of multispectral, high spatial resolution over VNIR-TIR bands observed by ASTER. 8 Acknowledgments This work would not have been possible without support from the U.S. Water Conservation Lab, NASA EOS Grant 03-OES-02, and the Department of Research, Development, and International Relations of the Purpan Graduate School of Agriculture. We especially thank John Eaton for Octave (www.octave.org) and the many contributors to the R project (www.r-project.org) for providing superb open-source software. ASTER NDVI over SMACEX site on 1 July 2002. NDVI values range from 0.3 (white) to 0.9 (black), respectively sparse vegetation and thick vegetation. E-C tower locations are indicated by white circles. Nested squares at upper right indicate sizes of 250-m and 1-km pixels. This scene and subsequent ones are projected to UTM zone 15 coordinates. The white area at the left denotes the edge of the ASTER scene. 33 2 ASTER surface temperature over SMACEX site on 1 July 2002. Temperatures range from 25 • C (dark gray) to 40 • C (white). Numbered points indicate field campaign E-C station locations. 34 3 NDVI histograms of SMACEX at 15-, 90-and 240-m resolutions. Bi-modal distributions for 15-m and 90-m images reflect the known land cover patterns of relatively sparse soybean fields and dense corn fields.

Fig. 1 .Fig. 2 .

 12 Fig. 1. ASTER NDVI over SMACEX site on 1 July 2002. NDVI values range from 0.3 (white) to 0.9 (black), respectively sparse vegetation and thick vegetation. E-C tower locations are indicated by white circles. Nested squares at upper right indicate sizes of 250-m and 1-km pixels. This scene and subsequent ones are projected to UTM zone 15 coordinates. The white area at the left denotes the edge of the ASTER scene.

Fig. 4 .Fig. 5 .

 45 Fig.4. Spectral properties within longer wavelengths of the TIR window. Lines represent, from top to bottom: example dry silt-loam soil emissivity (Emis.), from central Oklahoma, example atmospheric transmissivity (Trans.), from a SMACEX radiosonde, filter response (Resp.) functions for ASTER bands 13 and 14 (solid) and TM7 band 6 (dashed). Response functions have been scaled by 1 3 to aid display.

Fig. 7 .

 7 Fig. 7. Temperature vs. albedo distribution plot over SMACEX, 1 July 2002 showing mostly constant albedo values for temperatures between 30 and 40 • C . This distribution falls within the radiative controlled domain.

  05

  Figures

  evaluate model flux estimates independently of each other, comparisons with ground measured fluxes are needed. Using remote-sensing derived model fluxes and footprint-weighted ground based E-C measurements, flux compo-

	nents are compared for 8 ground locations within the SMACEX experimental
	area (Table 1, see locations in Fig. 2). The table is divided into four groups,
	one for each flux component (H, LE, G and R n ). Within each group are 1/2
	hourly averaged ground flux values (O) and TSEB (T) and SEBAL (S) flux
	values. The leftmost column contains the ground station identifier, appended
	by crop cover (C=corn, S=soybean). The bottom row contains mean flux es-
	timate deviations from E-C measurements when considering all 8 stations.

Ground observations for R n are from tower-mounted net radiometers. G values are from soil heat flux plates. Preliminary turbulent fluxes, H and LE,

 Boundary-Layer Meteorol. 50, 355-373. Seguin, B., Becker, F., Phulpin, T., Gu, X., Guyot, G., Kerr, Y., King, C., Lagouarde, J., Ottlé, C., Stoll, M., Tabbagh, T., Vidal, A., 1999. IRSUTE:A minisatellite project for land surface heat flux estimation from field to regional scale. Remote Sens. Environ. 68, 357-369. Sellers, P., Hall, F., Kelly, R., Black, A., Baldocchi, D., Berry, J., Ryan, M., Ranson, K., Crill, P., Lettenmaier, D., Margolis, H., Cihlar, J., Newcomer, J., Fitzjarrald, D., Jarvis, P., Gower, S., Halliwell, D., Williams, D., Goodison, B., Wickland, D., Guertin, F., 1997. BOREAS in 1997