Simple models including energy and spike constraints reproduce complex activity patterns and metabolic disruptions

In this work, we introduce new phenomenological neuronal models (eLIF and mAdExp) that account for energy supply and demand in the cell as well as the inactivation of spike generation how these interact with subthreshold and spiking dynamics. Including these constraints, the new models reproduce a broad range of biologically-relevant behaviors that are identified to be crucial in many neurological disorders, but were not captured by commonly used phenomenological models. Because of their low dimensionality eLIF and mAdExp open the possibility of future large-scale simulations for more realistic studies of brain circuits involved in neuronal disorders. The new models enable both more accurate modeling and the possibility to study energy-associated disorders over the whole time-course of disease progression instead of only comparing the initially healthy status with the final diseased state. These models, therefore, provide new theoretical and computational methods to assess the opportunities of early diagnostics and the potential of energy-centered approaches to improve therapies.


Introduction
Brain metabolism, even in its resting state, constitutes a major source of energy consumption in mammalian species. Indeed, cells -and especially excitable cells such as neurons -undergo constant ion fluxes both along and against the concentration and electric gradients. To move ions against these gradients, an active mechanism is required, which consumes energy in the form of ATP. In cells, this work is mostly associated with the sodium-potassium pump (Na/K pump or NKP) which moves 3 sodium ions out of the cell in exchange for 2 potassium ions moving in for every hydrolyzed ATP molecule, thus creating a net electric current [1]. As a result, Na/K pump is responsible for roughly 75% of the total energy consumption in neurons [2], which arguably makes it one of the most important players in the cell: its action makes the energy from the hydrolysis of ATP available to most other processes [3], allowing changes in the membrane potential, regulation of the volume, or transport of nutrients inside the cell. Thus the energy level, through the Na/K pump activity, modulates neuronal response and directly influences information processing [4].
Though the Na/K pump has been thoroughly researched in the past decades [3,1], surprisingly few neuronal models include the pump and its electrogenic properties [5,6,7] and even fewer account for its underlying energy substrate [8,9]. A probable reason for this fact comes from the significant focus of theoretical studies on cortical areas that generally display sparse activity. Such conditions put little or no metabolic stress on the neurons and thus limit the influence of the Na/K pump and energetic constraints on the dynamics. However, the story changes drastically when energy-intensive behaviors such as bursting or fast pacemaking dynamics are considered, or when studying neuronal disorders. Indeed, both situations can place neurons under significant metabolic stress and induce fluctuation in the metabolite and ion concentrations which, from NKP-driven coupling between metabolism and activity, can then lead to major changes in the neuronal dynamics.
Outside of neuroscience, the influence of Na/K pump and energy consumption on activity and disorders were investigated in the context of the cardiac electrophysiology [10,11,12]. However, awareness is now raising in the neuroscience community, including its most theoretically-oriented members, as an increasing number of publications start to stress the critical influence of mitochondria [13,14] and Na/K pump [15] and the intricate feedback loops between activity and energetics. Some well-known works on energetics in computational neuroscience include the energy budgets from [16] and [2], as well as studies related to the link between action potential shape and ATP consumption [17,18]. Yet, these studies deal with general budgets from the point of view of optimality theory and do not describe the local interactions between energy levels and spike initiation.
The interactions between energetics and neuronal activity are most visible in neuronal disorders such as epilepsy [19,20,21], Alzheimer [22], or Parkinson's disease [23,24]. It is therefore in the context of neuronal diseases that one can find the few studies that really focused on these interactions [9,8,25,26]. Unfortunately, because such studies are still scarce and the associated modeling frameworks are still limited, computational studies of neuronal disorders currently suffer from at least one of the following issues: a) they do not account for energetic constraints, b) the models do not reproduce important features of the relevant neuronal behaviors, or c) the size of the simulated networks is extremely small (notably due to the use of complex conductance-based models).
Here we present new models to help tackle these issues through theoretical descriptions of neuronal dynamics that a) account for energy levels and their influence on neuronal behavior, b) are able to reproduce most relevant neuronal dynamics in the context of disorders such as seizures or Parkinson's disease, and c) can be used in simulation of networks up to several million neurons.

Methods
In the following, we describe the implementation of the new models. We discuss the biological mechanisms that gave rise to the variables and equations in our models and list the associated properties that an energetic model should satisfy. Two major biological components considered in the models are the pumps that degrade ATP into ADP to maintain ion gradients (most notably the Na/K and calcium pumps), and ATP-gated potassium (K-ATP) channels that open or close depending on the ATP/ADP ratio. When ATP concentration or the ATP/ADP ratio decreases, the pump's effectiveness decreases, resulting in a rise of sodium concentration, thus increasing membrane potential and sometimes excitability. Conversely, a decrease in the ATP/ADP ratio tends to open K-ATP channels, allowing potassium to flow out of the cell and decrease membrane potential and excitability. These mechanisms directly or indirectly influence a neuron's excitability and its ability to generate action potentials. Depending on their relative importance, a neuron can, therefore, end up in a depolarized or hyperpolarized state when energy levels decrease.
Based on these main mechanisms, the models were design to meet several conditions that can be split into 1) behavioral requirements, associated to the type of responses and biological situations that the models can account for, and 2) practical constraints associated to the computational cost and theoretical complexity of the model. As energetic constraints are especially relevant for behaviors associated with diseased or hypoxic state, we designed our models so that they would be able to provide meaningful behaviors in such conditions.
Regarding behavioral requirements, we took care of reproducing the effects ATP/ADP changes on pumps and K-ATP channels so that the models could account for three major observations: • as mitochondrial health or metabolic resources decrease (e.g. during hypoxia), the excitability and resting potential of the neuron can increase [27,25], notably due to insufficient activity of the Na/K pump, • decrease in metabolic resources is also associated with an increase in calcium levels [27] due to insufficient activity of the pumps, • during seizures, or when submitted to excessive excitation, neurons undergo depolarization blocks characterized by "superthreshold" membrane potential without spike emission [28] that is caused by sodium-channel inactivation as the Na/K pump cannot move sodium out quickly enough.
In addition, the specific form of the equation was chosen to allow two specific behaviors to be switched on or off depending on the parameters used: • neuronal bistability, observed in several brain regions [29,30], is involved in important mechanisms such as up-and-down states and could also explain discontinuities in the progression of neurodegenerative diseases [31], • adaptation currents and bursting or rebound activities that are major players in neuronal disorders [32,33].
Our central goal is to develop models that do not only reproduce important behaviors, but also allow for large-scale event-based simulations. To achieve this, the computational cost and complexity of the models should be minimal. Thus, we decided to use hybrid models based on the integrate-and-fire paradigm.
We established that models including an adaptation current, such as the Quadratic Integrate-and-Fire and the AdExp neurons [34,35], were able to provide most of the required dynamics such as bursting and rebound activity [36,37]. The missing requirements -depolarization block and bistability -as well as the inclusion of metabolic resources would thus come from the addition of dynamic resource availability (broadly called energy in the following), as shown on Figure 1. This purpose of this variable is to represent the ATP/ADP ratio in biological neurons, though the phenomenological nature of the models implies that there are limits to this analogy.
For applications where bursting behavior and adaptation do not play an important role, a simple model that accounts only for energy dynamics is provided: the eLIF neuron. It introduces energy dynamics as an addition to the simpler leaky integrate-and-fire (LIF) model and enables us to analyze the consequences of these constrains in a more straightforward and visual manner. The behavior of this model can also be fully investigated analytically compared to the 3-dimensional system that arises in a second time when both energy and adaptation dynamics are considered. This second model, called mAdExp, is built upon the AdExp equations and cam reproduce all desired behaviors. Though analytical analysis of this model can prove challenging, most of its dynamics can be understood from the complementary analyses of the eLIF and AdExp models. Variables and interactions that must be present in the models to capture all relevant behaviors, the main molecules associated to each of the variables are also displayed. The type of interaction is marked on the arrow. For instance, w modulates (M) V as it influences the intrinsic dynamics of V but does not usually cause it directly. On the other hand, as changes in the membrane potential are the main cause of variations in w, V is said to drive (D) w. Eventually, all mechanisms consume (C) energy.

Introducing energy: the eLIF model
The first proposed model is a straightforward modification of the standard Leaky Integrate-and-Fire (LIF) model [38]. In order to provide an intuitive and analytically tractable implementation that would illustrate the consequences of energy dynamics and the constraints it places on spike-emission, we developed a twodimensional dynamical system describing the evolution of a) the membrane potential V of a point neuron and b) the available amount of energy that the neuron can access. To make the equations more readable and the parameters easy to interpret, the model is displayed using three equations; however, it can be easily simplified to a system of two equations only.
As in other standard integrate-and-fire models, the neuron possesses a leak potential E L , a membrane capacitance C m , and a leak conductance g L , the combination of the last two defining the membrane timescale τ m = C m /g L . Input from other neurons are represented by I syn while external input currents are associated to I e . When either of these inputs brings the neuron above its threshold potential V th , provided that there is enough energy ( > c ) a spike is emitted and the voltage is instantaneously reset to V r .
The available energy is introduced as a proxy for the ATP/ADP ratio in biological neurons. Its value varies with a typical timescale τ e and is regulated by an interplay of the energy production (which tries to maintain it close to the typical energy value defined by the energetic health α 0 ) and two consumption mechanisms. Phase space of the eLIF model in bistable parameter regime. V -nullcline is given by the blue line, -nullcline by the red curve. Fixed points (FPs) are shown by the circles (filled for stable and empty for unstable) and the cross marks the inflection point of the -nullcline. Dashed lines represent the shifts in the V -nullcline which lead to the disappearance of the unstable fixed point and of one of the stable fixed points (saddle-node bifurcation via the external current I e ). The super-threshold region, where spikes are elicited upon entrance, is marked by the light grey shading; the energy-limiting region ( < c ) is marked by the grey shading and overlaps with the super-threshold region in the dark grey area, where energy limitations prevent spiking though the neuron is above threshold.
The production term reflects mostly the oxidative phosphorylation performed in mitochondria [39] that enables a tight regulation of ATP levels in the cell.
The first consumption mechanism is associated with the fluctuations of the membrane potential and accounts for the ATP consumed by the Na/K pump to maintain ion homeostasis [3]. Since there is no available information about the functional form of the relationship between membrane potential and energy consumption, we have almost no constraints on the choice of the functional class. We selected a function allowing for a wide range of behaviors as observed in experiments while remaining as simple as possible: a 3rd order polynomial (see Figure 2, red line). Indeed, this is the simplest nonlinearity that can, depending on parameter values, either lead to a behavior that is qualitatively equivalent to a linear relationship or to the presence of a bistability, making it possible to asses the influence of bistable states on neurons' and network dynamics. The parameters defining the shape of the nullcline are: the flex potential E f (that corresponds to the inflection point, or flex, of the curve) and the energy-depletion potential E d , that is a potential at which -nullcline crosses the x-axis -E d thus corresponds to the lowest voltage-clamp potential that will lead to complete energy depletion and therefore neuronal death. The second source of energy consumption is the energetic cost δ of the spike generation mechanisms. Though the biological reason for this energy consumption is the same as the first term (ionic transfer by the Na/K pump), a separate term is necessary because of the reset mechanism of integrate-and-fire neurons: in such models, part of the spike duration is compressed into an instantaneous jump; δ thus accounts for the energy consumed during this compressed period. The normal energy level that the neuron is able to maintain depends on its "energetic health" described by the α parameter: a healthy neuron would have a value of α equal to one, while diseased neuron would see their α parameter decrease towards zero.
Contrary to most previous models, the leak potential is not constant, as it depends on the energy level of the neuron. The steady-state value E L of the membrane potential thus varies linearly, starting from E u when the energy is zero and decreasing as increases, crossing the potential E 0 for = 0 (see Figure 2) for details). Biologically, this account for the fact that a decrease in energy availability inhibits the function of the Na/K pump, leading to sodium accumulation inside the cell and thus to a depolarization. The behavior of the standard LIF is recovered when E u = E 0 and δ = 0.

Adaptation and bursting: mAdExp model
In order to model the whole range of biologically-relevant behaviors that can be observed in neuronal disorders such as epilepsy or Parkinson's disease, it is necessary to include a modulatory mechanism to account for cellular and spike-driven adaptation. This second dynamical system keeps the basic properties introduced in the eLIF model and extends them to accommodate the cellular adaptation and spike initiation mechanisms of the adaptive Exponential Integrate-and-Fire model (aEIF or AdExp) by [35]. This leads to a 3D model with three dynamical state variables which are the membrane potential V , the energy level (as for the eLIF model), and an adaptation current w: Compared to the eLIF implementation, the presence of the spike initiation mechanism through the exponential function removes the necessity of a hard threshold for spike prevention due to energy limitation: the ( − c )/ 0 factor suppresses the exponential divergence as soon as the amount of available energy goes below c .
The dynamics of the variable remains mostly unchanged except for the addition of a new consumption term associated with the adaptation current w: biologically γ −1 corresponds to the energetic cost of bringing back the potassium ions which exited the cell (through calcium-gated potassium channels) per pA unit of the adaptation current. The model thus clearly separates the contributions of the energy ( ) and of the calcium-gated adaptation (w).
Compared to the original AdExp model, the w dynamics includes an additional term, c c +2 I KAT P , to account for ATP-sensitive potassium channels that trigger potassium outflow when the ATP/ADP ratio becomes small, with a typical activation-threshold depending on the ADP/ATP ratio [40]. I KAT P is thus the maximum current at zero energy. Because of the numerous calcium exchangers in neuronal cells [41,42], the term responsible for the exponential decay of the adaptation current with timescale τ w is considered to be energy-independent. Thus, only E L and K-ATP induce energy-dependent changes in the adaptation current.

Numerical implementations
Implementations of the models are available for three major simulation platforms: NEST [1], through the NESTML language [44], BRIAN [45], and NEURON [46]. Models are available on ModelDB and on GitHub 1 , together with code to reproduce the figures. Networks were generated using NNGT 2.0 [4] and simulated using NEST 2.20 [1].

Fitting procedure
To reproduce experimental recordings, we could set some of the model parameters directly from the data. The rest had to be manually adjusted. The following parameters can be informed from the data: a) E L was obtained by measuring the median resting value b) the membrane timescale τ m was measured from the initial slope of the membrane dynamics in response to hyperpolarizing currents c) the sum g L + a was obtained through a linear regression from the difference between resting E L and steady-state E ss potentials in response to depolarizing currents as ∆V = E ss − E L = I a+g L . These properties were used to constrain the following parameters: C m , g L , a, E L , E 0 , E u . All other parameters were then manually adjusted to minimize the discrepancy between subthreshold dynamics, number and time of spikes. Further research would be necessary to find how to automate this procedure using a proper distance function in optimization toolboxes.

Results
The new eLIF and mAdExp models enable us to obtain a variety of new dynamics such as rebound spiking, depolarization block, cellular bistability and up-and-down states, as well as biologically relevant transitions from a healthy to a diseased state.
For hybrid models, most of the neuronal dynamics can be understood through two main concepts: a) fixed points (FPs), which are equilibrium states of the model, and b) bifurcations, which are sudden changes in the number or stability of the fixed points, and which make the neuron change its behavior, for instance from resting to spiking.
This section details the aforementioned behaviors and their mechanistic origins through the theory of dynamical systems, using fixed points and bifurcations.

Behaviors and bifurcations of the eLIF model
The eLIF model, like the integrate-and-fire (LIF) neuron, has only two dynamical states: quiescent or active (spiking). Due to the energetic constraints, the model has two possible quiescent states which are the "normal" resting state, with a membrane potential located below threshold, and a super-threshold state where depleted energy levels prevent spike emission. The finite energy resources also imply that, contrary to the LIF neuron, the active state can be transient, as the neuron transits from its resting state to a quiescent, super-threshold state through an active period.
In the language of dynamical systems, the quiescent states are associated to FPs inside the continuous region (if either V F P < V th or F P < c ), whereas the active state is associated to the absence of a stable FP that can be accessed continuously in the region of phase-space where the neuron lies -see Figure 3.
We will focus here on the situation that is most relevant for the study of neuronal disorders, i.e. the case where E u > E 0 , meaning that decrease in energy levels leads to increase in membrane potential. This situation leads to a neuronal behavior which is that of an integrator; another type of behavior, closer to that of a resonator, with dampened oscillations is also possible for E u < E 0 and is discussed in section 3.4 and in the Supplementary Information.
In this situation, due to the nonlinearity of the -nullcline, the biophysically acceptable domain for steady states ( ≥ 0 and V in a reasonable range of potential) can contain either zero, one, or three FPs. In the case of a single, necessarily stable FP, it corresponds to a standard neuron with a single resting state. For certain combinations of the neuronal parameters, the V -nullcline can intersect the -nullcline three times, leading to two stable FPs and one unstable point. This situation corresponds to a bistable cell, where two distinct resting states are possible: an up-state, characterized by lower energy levels and high membrane potential, and a down-state, associated to higher energy and hyperpolarized membrane potential. Responses of the bistable neuron to the different step-currents are illustrated in Figure 3. Depending on initial state and the input the neuron transitions between the up-and down-states. Finally, the situation without FPs in the biophysical domain is unsustainable and will lead to rapid neuronal death. Possible reasons for transitions between these states will be detailed in the following section.  Table S2 in Supplementary Information E.1 for detailed parameters.
We use the transitions in the number of FPs, called bifurcations, to predict the behavior of the neuron. The bifurcations can have two separate kinds of consequences, that can potentially happen simultaneously: a) a change in the steady-state behavior of the neuron such as the switch from a unistable to a bistable state or vice-versa, b) a transition from a quiescent to an active state.
Let us discuss these bifurcations in response to an external stimulation associated to an applied current I e . The consequence of I e is to shift the V -nullcline horizontally (towards more negative potentials if I e < 0, or towards more positive if I e > 0), which can lead to transition between the unistable and bistable states as one stable FP either splits into one stable and one unstable FP or, on the contrary, merges with the unstable FP and disappears. This type of transition is called a saddle-node bifurcation and occurs for: Depending on the value of I e , the neuron can thus display either a single or two stable FPs -see Figure 3 and Supplementary Information B.3 for the analytic derivation of the FPs.
As I e increases, the transition from three FPs to one FP can also lead the neuron to fire, either transiently if the remaining FP is located in the continuous region (if either V F P < V th or F P < c ) or continuously (if V F P ≥ V th and F P > c ).

Transition from health to disease
As energy availability decreases, either due to disease [25] or hypoxia [27], neurons often display a parallel increase in their resting membrane potential and excitability, which can lead to highly active periods before the neuron ends up in a highly depolarized yet completely non-responsive state also called depolarization block. Biologically, this low-energy state -(d) and below on Figure 4 -would be associated to deregulation of calcium levels and accumulation of oxidizing agents which eventually lead to cell death (occuring when α reaches zero in the model).
Due to the interaction between energy and membrane potential in the eLIF neuron, the model can reproduce this kind of dynamics through the evolution of one or more parameters. The most straightforward way to model this transition is through the α parameter which represents the energetic health of the neuron -see Figure 4. The progressive decrease in the value of α, from values close to 1 for a healthy neuron to values that tend towards zero for a diseased cell, leads to progressive changes in the membrane potential and excitability of the neuron. The typical behavior of the model, illustrated on Figure 4, consists of a slow increase of the resting membrane potential, and thus of the excitability, until the background noise or external input is sufficient to trigger spike emission from the neuron. Once that happens, the cell enters a highly active state in which it remains until the progressive decrease of α brings the target energy below c , at which point spike emission stops and the neurons enters a highly depolarized and non-responsive state.
When it comes to collective dynamics, a decrease in neuronal health can go unnoticed especially if the homeostatic regulation adjusts excitability of individual neurons. It happens, for instance, in excitatory and inhibitory networks displaying asynchronous-irregular (AI) activity with low firing rates -see Figure 5.
Without external input, the distribution of firing rates as well as the average properties of the activity (cross-correlations and coefficients of variation) can remain stable despite neuronal health decrease ( Figure 5, panels B and C). This happens because compensatory mechanisms enable neurons to maintain firing rate homeostasis by means of synaptic scaling and regulation of cell excitability, that we modeled numerically by a decrease of excitatory synaptic weights and an increase of V th -see Supplementary Information E.2 for more details.
However, the response to an external input can be drastically modified ( Figure 5, panel A), transitioning from an almost continuous tonic response (top), to an intermittent, bursty dynamics (bottom). This example demonstrates how, depending on the homeostatic capabilities of the brain region of study and the recording protocol, the effect of energetic constraints can be either masked or clearly visible in the neuronal responses. There are multiple ways in which the energetic health can influence the information processing capabilities. Using our models these mechanisms can be studied further in large recurrent networks.

Dynamics of the mAdExp model, biologically-relevant behaviors
Despite the multiple interesting features of the eLIF model, several important dynamics such as bursting or adaptation cannot be reproduced within the model. In order to recover all relevant behaviors, we added a spike-generation mechanism as well as an adaptation current to the eLIF model to obtain the mAdExp model (modified AdExp with energy dependency).
This 3-dimensional model is then able to provide all the features of the eLIF and AdExp models while bringing the dynamics closer to biological observations, especially in large-input or stress-inducing situations. Figure 6 shows several standard neuronal responses reproduced by the model, as well as how these responses evolve as the input intensity increases up to values where the neuron cannot sustain continuous activity.
Though the theoretical analysis of the model becomes more complex, "standard" resting states 2 for healthy neurons can be very well approximated by the fixed point of the eLIF model because the adaptation current is usually close to zero at rest. Furthermore, their response to low-intensity stimuli can be accurately predicted 2 "standard" meaning that V F P is several ∆ T smaller than V th − ∆ T ln Eu−E 0 . Similar responses to the lower (yellow) currents can be achieved by the original AdExp model. However, each of these dynamics now comes with an "energy-depleted" state for high input current (red), associated to a depolarization block (responses associated to red bars), that cannot be captured by AdExp model. In addition to these standard behaviors, dynamical repertoire of the mAdExp neuron also includes a different mechanism for post-inhibitory rebound spiking (IR), and can display post-excitatory rebound (ER) or intermittent spiking dynamics (IS). See Table S3 in Supplementary Information E.1 for detailed parameters. by the AdExp model with the same common parameters and the corresponding E L value 3 . Most healthy neurons thus share the bifurcations associated to the AdExp model [36,3], with the notable addition of a new bifurcation for rebound spiking which will be developed in the next section.

Rebound spiking mechanisms in the different models
Rebound spiking is a common property in neurons, with is potentially significant in epilepsy [49] and for information processing, be it in the striatum [50], the thalamocortical loop [51], or in auditory processing [52] and grid cells response generation [53,54]. This mechanism, though already available in several models such as AdExp [35], strongly restricts the responses of the neuron such that only a fraction of the typical dynamics of rebound-spiking neurons can be recovered. The reason is that, in the AdExp model, rebound bursting is always associated to a sag and significant adaptation -see conditions in [3] and Supplementary Information D -and therefore cannot reproduce either non-sag subthreshold responses or some spiking behaviors associated to excitatory inputs, cf. Figure 7.
The mAdExp model provides two new ways of extending the variety of rebound behaviors that can be modeled: a) by introducing a new mechanism for rebound spike generation without inhibitory sag and b) through the energy dynamics, leading to less significant sags and lower excitability compared to the adaptation mechanism -see also Figure S9 in the Supplementary Information.    Table S4 in Supplementary Information E.1 for detailed parameters.
Rebound spiking in mAdExp can occur through a new bifurcation for Eu − E 0 < ∆ T and V th sufficiently low (see Supplementary Information D for details) which leads to the positive divergence of the V -nullcline before V th and thus to the existence of a stable fixed point such that V th > V F P > V * , with Figure 7 shows how the mAdExp model can successfully reproduce complex behaviors found in the Allen Cell Types Database 4 such as rebound bursting with little to no sag 5 (A.2) or cells displaying both rebound spiking and rapid depolarization block 6 (B.2). Due to the mAdExp properties, the possibility of rebound dynamics is thus extended compared to the AdExp model and can be obtained with or without sag, as well as with or without spike adaptation.

Choices underlying the models
The eLIF and mAdExp models where chosen as integrate-and-fire models because of the analytic simplicity of such equations and their computational efficiency compared to conductance-based models. Indeed, the straightforward detection of spike times in such models makes them especially suited for simulations of large scale neuronal networks using standard spike-based simulators and their discontinuous dynamics makes bursting possible with only two equations instead of three for continuous models like the Hindmarsh-Rose neuron [55] or general conductance-based models.
Though our models are almost completely phenomenological, their parameters can be directly related to biological phenomena, often even in a quantitative manner, enabling precise predictions from their theoretical analysis. The objective of obtaining single-neuron models where the variables can be interpreted and mapped in a straightforward manner also prevented us from working with previous models such as the Epileptors [56,57] or Model 2 from [26] -see discussion below.
To obtain a model capable of reproducing all the behaviors that we deemed necessary, the mAdExp model was derived from the AdExp neuron [35] and not from other well known implementations such as the QIF, first proposed by Izhikevich in his seminal paper [34]. This choice was made because, despite some obvious drawbacks regarding the more complex analytics and slightly slower integration of its exponential term, the AdExp model best reproduces the I-V curve of neurons and the dynamics of spike initiation [58], and is exempt of some of the mathematical shortcomings of the QIF model [59].
The variable was designed to qualitatively reproduce biological mechanisms and behaviors associated to neuronal metabolism. However, the complexity of these mechanisms led us to choose a strongly reductionist approach to reproduce some of the features that came out of previous studies using more detailed conductancebased or multicompartmental models [9,25,26,60]. Therefore, though it can be qualitatively mapped to some specific mechanisms, cannot be quantitatively related to any biological measurement.

Novelty and biological relevance
The eLIF and mAdExp neurons are the first integrate-and-fire models to provide an unambiguous description of energy dynamics, enabling to investigate its consequences in single-cell or in recurrent-network configurations. Indeed, contrary to previous models where slow variables were usually introduced to model adaptation from calcium-gated potassium and bursting 7 , the eLIF and mAdExp models provide the variable as a way to explicitly model energy-related and spike-initiation constraints. Though other implementations of models including slow variables might be able to reproduce some of the behaviors examined here, to the best of our knowledge, the eLIF and mAdExp are the only phenomenological models that permit the investigation of feedback loops between energy levels, neuronal excitability, and spike emission. The fact that the model's variable are interpretable and directly linked to biological mechanisms enabled us to extend the AdExp model in a straightforward manner; this will also let others expand the models if they need to capture additional mechanisms or external interactions, e.g. with glial cells.
The only other examples of models with an explicit variable representing energy levels we found were developed in [61] and in Model 2 from [26]. While the former is quite simple and not connected to any specific biological mechanism, the latter explicitly presents the second variable as a proxy for the ATP concentration in the neuron. In Model 2, the interpretation of A as a proxy for ATP level notably stems from the fact that it was designed as a simplification of the more detailed Model 1, a conductance-based model that included K-ATP channels, where the A variable was defined as the ATP concentration. The model provides interesting dynamical properties and enabled the authors to develop a new way of modeling the neuro-glio-vascular system. However, if one's purpose is to investigate dynamics where both calcium-gated potassium adaptation and energetic constraints are involved, then one would not be able to use Model 2 as it does not provide a clear distinction between these two mechanisms. Indeed, the newly introduced A variable only influences the value of the threshold and is therefore quite close to a GLIF model [62], which makes it impossible to separate effects that would biologically stem from "standard" adaptation mechanisms, associated to calcium-gated potassium currents, and effects that would be specific to ATP-related dynamics.
The implementation chosen in the mAdExp model solves this issue by establishing a clear separation between the retroactions associated to adaptation and those related to changes in energy levels. In this model, the effects of and w on the membrane potential can be opposite and occur (in general) on different timescales. The variable also regulates the spike initiation mechanism of the neuron, meaning that, contrary to previous models, energy depletion may render a neuron totally unresponsive regardless of the input strength. In addition, our models consider all sources of energy consumption including spikes and subthreshold ion currents.
Thus, over long timescales, the parameter qualitatively accounts for energy availability as the ATP/ADP ratio to which pumps and channels are sensitive [63,39], with this sensitivity summarized in the I KAT P parameter of the mAdExp model. Contrary to slow current variables that can vary arbitrarily into the positive and negative realm, represents an energy stock that must remain positive for the neuron to survive: if the neurons encounter conditions where reaches zero with V ≥ V d , the models provide an explicit condition for neuronal death.
Over shorter timescales, sharp decreases in following spike emission can lead to depolarization block. This phenomenon is mostly associated with sodium channel inactivation in neurons, and is caused by a sodium accumulation that is too quick to be compensated by the Na/K pump. Though it is not directly related to energetic constraints, we consider that having this mechanism associated to the "energy" variable makes sense because the timescale of sodium channel inactivation depends on the resting sodium levels, which in turn depend on the energetic health of the neuron: a neuron with a very active pump would be able to sustain more spikes than one with a defective pump. This mechanism is related to the δ and c parameters in the models.
Finally, as represents the ATP/ADP ratio, the α parameter quantifies the metabolic and mitochondrial state of the neuron and can be used to investigate the transition from health to disease as exemplified in the Results. A decrease of the α parameter can be related to metabolic insults associated to either mitochondrial defects [64,65], a decrease of oxygen or glucose availability, or the buildup of various molecules such as reactive oxygen species (ROS) [66,67] that prevent proper metabolic homeostasis.

Consequences of the V/ relationship
One of the major features of the model is the interaction between the energy level and the resting potential of the neuron. This interaction can lead to a transition from "healthy" or "optimally responsive" neurons to "diseased", non-responsive neurons. Interestingly the neuron may go through a hyper-excitable state during this transition, meaning that disease progression can be marked by a broad range of neuronal dynamics and properties. Because changes in the energy level affect the neuronal excitability, the synchronizability and information processing properties of the neurons change significantly as their available energy decreases. This property of the model matches observations in various neurodegenerative diseases. Synchronizability notably changes in Parkinson's disease (PD), for instance, where oscillations in the beta range (13-30 Hz) become predominant and are thought to be involved in some motor symptoms. Though known variations in the connectivity strongly influence this dynamical change, modification of intrinsic neuronal properties due to metabolic insult are also likely to contribute to the transition towards more synchronized activity [32,68]. Even more obvious, epileptic seizure are characterized by excessive or hypersynchronous neuronal activity and their onset and termination are likely to be related to the metabolic state of the neurons [69,19,20]. Finally, the transition through an hyperactive phase before entering the non-responsive depolarized state has also be proposed for diseases such as ALS [25].
From an information transfer perspective, the positive retroaction between depolarization and energy depletion can lead to increased false positives due to hyperexcitable neurons in diseased conditions. Furthermore, because of the necessity of a minimum "metabolic level" for spike emission, this also means that energy-impaired neurons cannot sustain long-term responses, and would tend to display phasic responses. These combined effects could further drive bursty activity such as what is observed in PD, where the reliability of thalamic relay breaks down and the cells start emitting bursts of activity which could lead to tremor [70,71].
The mAdExp model can reproduce the main relevant dynamical properties in these phenomena and therefore enables detailed and potentially large-scale computational studies. Such simulations could lead to more realistic dynamical models and thus to new experimentally testable predictions.

Limitations
Due to their simplicity, the eLIF and mAdExp models still suffer from many of the limitations of the original LIF and AdExp models.
For example, the eLIF cannot reproduce bursting behavior and can only exhibit simple accelerating or decelerating spiking patterns. Though the dynamical richness of mAdExp is greater than the LIF and AdExp models, its adaptation mechanism also possesses the same drawbacks as the original model: the presence of a single adaptation timescale τ w . As for the AdExp and QIF, though, the mAdExp can be extended in a straightforward manner to account for multiple adaptation timescales by adding additional w-like variables.
Since multiple biological phenomena are associated to or can affect the variable (Na f inactivation, ATP/ADP ratio, oxygen concentration, pH, ROS. . . ), precise experimental verifications and relations to biochemical pathways can be quite complex or even impossible to predict, at least if several phenomena are occurring on similar timescales. The depolarization block, for instance, only stems from the combination of metabolic parameters δ and c in our model. It is indeed strongly related to sodium or potassium accumulation (though not 'metabolic' per se, these directly depend on the efficiency of the NKP), to general energetic considerations [72,73], and may have been selected due to energetic constraints [74]. Yet, other slow mechanisms that are not accounted for in our models also contribute to this behavior in biological neurons; notably chloride-related changes in cell volume [75]. However, the purpose of our models is to explore how changes in energy availability may increase or reduce the occurrence of specific behaviors such as the depolarization block. Thus, for the sake of simplicity, we did not include such passive mechanisms as they do not directly influence energy availability.
Eventually, complex interactions between sodium or calcium levels and ATP production [76,77] is only coarsely implemented in the model. In particular, because the adaptation variable w represents calcium-gated potassium, and not directly the calcium levels, interactions between and w would not capture precise biological mechanisms. Overall, calcium dynamics can have very different impacts on ATP production, depending on concentrations and timescales, which cannot be completely accounted for by the simple relationship present in the model.

Conclusion
The two models introduced in the present study provide a novel reductionist approach to include generic energetic constraints and energy-mediated dynamics to the models of single neurons. The low-dimensional nature of these two dynamical systems makes them suitable for analytical investigation of energy-based bifurcations in neuronal behaviors, as well as for large scale simulations.
The mAdExp model, in particular, is able to replicate a large range of biologically-relevant behaviors as well as their evolution under metabolic stress. Complex behaviors that are crucial for some brain regions and disorders, such as rebound spiking or depolarization block, now can be successfully reproduced. Since energetics plays a critical role in many disorders, this model is especially well suited to explore possible origins of the differences observed between normal and diseased activities in neuronal populations.
Finally, these new models are not limited to the comparison between specific healthy or diseased states, as they provide a tunable parameter to represent neuronal health. Thus, the continuous transition between states can now be investigated, as well as dynamical feedback between activity and resource consumption in resource-limited conditions such as in neuronal cultures or seizures.

A. Benchmarks
The runtime of the models was measured using NEST 2.20 [1] and compared with existing implementations. The neurons were parametrized to spike at 25 Hz during 60 s and compared to a baseline run of 60 s without any neuron model. Table S1 compares the runtime of all models mentioned in this papers, as well as conductance-based neurons.
As can be seen from Table S1, the runtime of the models are similar to or faster than those of the AdExp and conductance-based models, while accounting for energy dynamics and displaying a larger variety of behaviors.

B. Fixed points and bifurcations of the eLIF model B.1. Nullclines
The two nullclines of the model are given by:  Table S1. Runtime of various models in NEST. A "baseline" run with no neuron (None), compared to runs with one neuron of each of the mentioned models. For the new energy-based models (eLIF and mAdExp), two runs were performed: one using a naive implementation and another using slightly optimized implementation (numbers in parentheses). Conductance-based models are also included: a standard Hodgkin-Huxley (HH) model which can display regular spiking an depolarization block, and one with calcium and calcium-gated potassium (HH+Ca) to reproduce bursting dynamics

B.2. Saddle-node bifurcation via I e
For a state where 3 FPs are present (see Figure 2), the coalescence of the higher stable FP, S + , and the unstable FP, U , occurs at a point B = (V B , B ), when the V -nullcline touches the 3rd order polynomial, i.e. when the local slope of the tangent to the curve is equal to Using also the second equation for V B , one gets the two critical values for I e = ±I * Which can be further simplified to give Equation 3.

B.3. General solution for the fixed points
The FPs of the eLIF model are the intersection of the two nullclines given by Equation S4. Writing out the equation for the FPs results in the 3rd order polynomial. From [2], we can get the general solution for the roots of this 3rd order polynomial in the case where E u > E 0 . Let us write it under the form ax 3 +bx 2 +cx+d, given x = / 0 Coefficients here are given by: Note that, though δ was used for coherence with [2], it is not related to the δ parameter which appears in Equation 1 and is associated with the spiking cost in the neuronal model.
which leads to At the bifurcation points If I e = I * e± , one recomputes δ as − 3 y N 2a to get its correct sign. This gives Then Single real solution In the case where I e ∈ [I * e− , I * e+ ] or E u ≤ E 0 , the single real root and is obtained through Cardano's formula:

C. Fixed points and bifurcations of the mAdExp model Nullclines
The nullclines of the mAdExp model can be expressed in multiple ways, among which:

Approximation of the fixed points
In this section, we consider parameter sets where the effect of I KAT P is negligible. As long as the fixed points have a value of V F P which is lower than V th − ∆ T , their value can be well approximated by replacing g L by (g L + a) in the solutions of the eLIF model (see previous section), then considering: Numerically, on can then converge iteratively towards an improved solution, starting from this initial guess FP 0 , then correcting the external current that will be used to compute FP i+1 by I e,i+1 = I e − w F P,i +

D. Behaviors
This section provides some additional information regarding the behaviors that can be obtained through the eLIF and mAdExp models. Figure S8 shows how different parameters can give rise to both type I and type II I − f curves. for V th > V F P , the neuron has a continuous type I response curve whereas for V th > V F P the curve, though still continuous, becomes closer to a type II curve, with a sharp increase starting immediately at the bifurcation current I * e . See Table S2 in Supplementary Information E.1 for detailed parameters.  Table  S2 in Supplementary Information E.1 for detailed parameters.

Rebound spiking/bursting
The following paragraphs show an example of "rebound activity" with the eLIF model ( Figure S9), as well as details about the conditions leading to rebound activity for the AdExp and mAdExp models.
AdExp For the AdExp model, rebound spiking occurs [3] either: • for type I excitability (a/g L < τ m /τ w ) • in all situations for type II excitability (a/g L > τ m /τ w ) a) if τ m /τ w < 1 or b) if τ m /τ w ≥ 1 Cases I.b and II.b correspond to a neuron exhibiting dampened oscillations, so the presence of the sag is obvious. For cases I.a and II.a, the faster timescale associated to the membrane potential conditions the presence of a sag. Because type II excitability with τm 4τw 1 − τw τm 2 > a/g L is impossible, as τm 4τw 1 − τw τm 2 < τm τw < a g L for τ m /τ w ≥ 0, this covers all cases. Thus, rebound spiking in the AdExp model is always associated to a sag.
This can also be shown mathematically for I.a and II.a by looking at the eigenvector associated to the lowest eigenvalue: For τ m /τ w < 1 (I.a and II.a), the denominator d of the x component of e − gives its sign, and since one can see that, as expected from the ratio of timescales, there is always an overshoot and a sag for I.a and II.a.
mAdExp The new rebound bursting behavior is associated to a positive divergence of the V -nullcline (cf. Equation S10). Since the divergence occurs for the positive sign is obtained for To get the mAdExp model to display rebound spiking and no sag one must combine the previous condition with the constraints of the AdExp and eLIF models: • the condition for no overshoot is a type I neuron with either τm τw > 1 or a g L > τm 4τw 1 − τw τm 2 , or any neuronal type with a ≤ 0 (note that for small values of a, the sag, though technically present, can be neglected for all practical purposes), • the condition for no overshoot due to energy dynamics is E u ≥ E 0 (necessary for the ).

E. Parameters
Detailed parameter sets used in the different figures can be found in the following tables.    300  200  250  300 400  300  100  200  100  250  pA  Table S3. Parameters used for the different behaviors of the mAdExp model on Figure 6.

E.2. Network simulation
These are the parameters used for the network simulations of Figure 5.  Table S5. Static neuronal parameters used with the all eLIF neurons in Figure 5.  Table S6. Specific parameters used for each of the simulations at a different neuronal health in Figure 5.
Each health level, corresponding to a value of α is associated to the corresponding values for V th and V reset in the same column.
The networks were all Erdös-Renyi graphs containing 800 excitatory neurons and 200 inhibitory neurons, each having an average degree of 100. There were generated using the NNGT library [4]. Even for the healthy situation with α = 1, the network was placed in the inhibition-dominated regime, with an inhibitory synaptic strength 120 pA and a synaptic timescale of 2 ms against excitatory synapses with strength 30 pA and a timescale of 0.2 ms. The inhibitory strength and all timescales were kept unchanged for all conditions, only excitatory strength was scaled according to the values s e in Table S6.