Energy-based models to describe complex neuronal activities and metabolic constraints

In this work, we introduce new phenomenological neuronal models ( e LIF and mAdExp) that account for energy supply and demand in the cell as well as their interactions with spiking dynamics. Through energetic considerations, these 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 e LIF and mAdExp enable large-scale simulations that are necessary 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
ion concentrations which, from NKP-driven coupling between metabolism and activity, can then lead to major 23 changes in the neuronal dynamics. 24 Outside of neuroscience, the influence of Na/K pump and energy consumption on activity and disorders   (Büeler 2009; Haddad and Nakamura 2015). 37 It is therefore in the context of neuronal diseases that one can find the few studies that really focused on these model (Brunel 2008). In order to provide an intuitive and analytically tractable implementation that would 85 illustrate the consequences of energy dynamics and the constraints it places on spike-emission, we developed a two-dimensional dynamical system describing the evolution of a) the membrane potential V of a point neuron 87 and b) the available amount of energy ϵ that the neuron can access. To make the equations more readable 88 and the parameters easy to interprete, the model is displayed using three equations; However, it can be easily 89 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 ca-91 pacitance C m , and a leak conductance g L , the combination of the last two defining the membrane timescale to I e . When either of these inputs brings the neuron above its threshold potential V th , provided that there is 94 enough energy (ϵ > ϵ c ) a spike is emitted and the voltage is instantaneously reset to V r .

95
The available energy ϵ varies with a typical timescale τ e and is regulated by a production term (which tries 96 to maintain it close to the typical energy value ϵ 0 ) and two consumption mechanisms. The first consumption 97 mechanism is associated to the fluctuations of the membrane potential. It is responsible for the nonlinear shape 98 of the ϵ-nullcline (see Figure 2, red line). The parameters defining the shape of the nullcline are: the flex 99 potential E f (that corresponds to the inflection point, or flex, of the curve) and the energy-depletion potential 100 E d , that is a potential at which ϵ-nullcline crosses the x-axis -E d thus corresponds to the lowest voltage-clamp 101 potential that will lead to complete energy depletion and therefore neuronal death. The second source of energy 102 consumption is the energetic cost δ of the spike generation mechanisms. The ability of the neuron to maintain 103 its energy levels close to ϵ 0 depends on its "energetic health" described by the α parameter: a healthy neuron 104 would have a value of α equal to one, while diseased neuron would see their α parameter decrease towards zero.

105
Contrary to most previous models, the leak potential is not constant, as it depends on the energy level of the 106 neuron. The steady-state value E L of the membrane potential thus varies linearly, starting from E u when the 107 energy is zero and decreasing as ϵ increases, crossing the potential E 0 for ϵ = ϵ 0 (see Figure 2) for details).

108
The behavior of the standard LIF is recovered when E u = E 0 and δ = 0.
Compared to the eLIF implementation, the presence of the spike initiation mechanism through the exponential 118 function removes the necessity of a hard threshold for spike prevention due to energy limitation: the (ϵ − ϵ c )/ϵ 0 119 factor suppresses the exponential divergence as soon as the amount of available energy goes below ϵ c .

120
The dynamics of the ϵ variable remains mostly unchanged except for the addition of a new consumption term 121 associated with the adaptation current w: biologically γ −1 corresponds to the energetic cost of bringing back 122 the potassium ions which exited the cell per pA unit of the adaptation current. where spikes are elicited upon entrance, is marked by the light grey filling; the energy-limiting region (ϵ < ϵ c ) is marked by the grey filling and overlaps with the super-threshold region in the dark grey area, where energy limitations prevent spiking though the neuron is above threshold.
Compared to the original AdExp model, the w dynamics includes an additional term, ϵc ϵc+2ϵ I KAT P , to account 124 for ATP-sensitive potassium channels that trigger potassium outflow when the ATP/ADP ratio becomes small. for the exponential decay of the adaptation current with timescale τ w is considered to be energy-independent.

128
Thus, only E L and K-ATP induce energy-dependent changes in the adaptation current.

144
The new eLIF and mAdExp models enable us to obtain a variety of new dynamics such as rebound spiking,  Table 2 in Appendix A.6 for detailed parameters.
number or stability of the fixed points, and which make the neuron change its behavior, for instance from resting  The eLIF model, like the integrate-and-fire (LIF) neuron, has only two dynamical states: quiescent or active 155 (spiking). Due to the energetic constraints, the model has two possible quiescent states which are the "normal" 156 resting state, with a membrane potential located below threshold, and a super-threshold state where depleted 157 energy levels prevent spike emission. The finite energy resources also imply that, contrary to the LIF neuron, 158 the active state can be transient, as the neuron transits from its resting state to a quiescent, super-threshold 159 state through an active period.

160
In the language of dynamical systems, the quiescent states are associated to FPs inside the continuous region , whereas the active state is associated to the absence of a stable FP that can 162 be accessed continuously in the region of phase-space where the neuron lies -see Figure 3. 163 We will focus here on the situation that is most relevant for the study of neuronal disorders, i.e. the case where 164 E u > E 0 , meaning that decrease in energy levels leads to increase in membrane potential. This situation leads 165 to a neuronal behavior which is that of an integrator; another type of behavior, closer to that of a resonator, 166 with dampened oscillations is also possible for E u < E 0 and is discussed in section 3.4 and in the Appendix.

167
In this situation, due to the nonlinearity of the ϵ-nullcline, the biophysically acceptable domain for steady 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.
bifurcations can have two separate kinds of consequences, that can potentially happen simultaneously: a) a 180 change in the steady-state behavior of the neuron such as the switch from a unistable to a bistable state or 181 vice-versa, b) a transition from a quiescent to an active state.

182
Let us discuss these bifurcation in response to an external stimulation associated to an applied current I e .

183
The consequence of I e is to shift the V -nullcline horizontally (towards more negative potentials if I e < 0, or 184 towards more positive if I e > 0), which can lead to transition between the unistable and bistable states as one 185 stable FP either splits into one stable and one unstable FP or, on the contrary, merges with the unstable FP 186 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 188 and Appendix A.3.3 for the analytic derivation of the FPs.

189
As I e increases, the transition from three FPs to one FP can also lead the neuron to fire, either transiently    Table 2 in Appendix A.6 for detailed parameters. . 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 3 in Appendix A.6 for detailed parameters.  Table 4 in Appendix A.6 for detailed parameters.  Rebound spiking in mAdExp can occur through a new bifurcation for Eu − E 0 < ∆ T and V th sufficiently low (see Appendix A.5 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 Allen Institute for Brain Science (2015). Allen Cell Types Database. Available from: celltypes.brain-map.org 5 cell ID 566978098: celltypes.brain-map.org/mouse/experiment/electrophysiology/566978098 6 cell ID 570896413 celltypes.brain-map.org/mouse/experiment/electrophysiology/570896413

Biological roots of the energy variable 246
Because of the strong reductionist approach chosen in designing these models, the ϵ variable cannot directly, 247 and especially not quantitatively, be related to any biological measurement. However, since the models were 248 made to reproduce biological mechanisms and behaviors, some qualitative analysis is possible.

249
Indeed, with respect to the transition from health to disease, as well as to the K-ATP channels (associated to

304
The two models introduced in the present study provide a novel reductionist approach to include generic ener- of 60 s without any neuron model. Table 1 compares the runtime of all models mentioned in this papers, as 326 well as conductance-based neurons.

327
As can be seen from Table 1, the runtime of the models are similar to or faster than those of the AdExp and 328 conductance-based models, while accounting for energy dynamics and displaying a larger variety of behaviors. The two nullclines of the model are given by:

A.3.2. Saddle-node bifurcation via I e 333
For a state where 3 FPs are present (see Figure 2), the coalescence of the higher stable FP, S + , and the unstable Using also the second equation for V B , one gets the two critical values for I e = ±I * Coefficients here are given by: Note that, though δ was used for coherence with Nickalls 1993, it is not related to the δ parameter which

At the bifurcation points
If I e = I * e± , one recomputes δ as − 3 y N 2a to get its correct sign.

354
This gives Then ϵ 1 = ϵ 0 (1 + r), ϵ 2 = ϵ 0 (1 − 2r) 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 355 Cardano's formula: The nullclines of the mAdExp model can be expressed in multiple ways, among which: 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 2 in Appendix A.6 for detailed parameters.

363
In this section, we consider parameter sets where the effect of I KAT P is negligible. As long as the fixed points 364 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

365
(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 367 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 + 368

370
This section provides some additional information regarding the behaviors that can be obtained through the 371 eLIF and mAdExp models. 372 Figure 7 shows how different parameters can give rise to both type I and type II I − f curves.

373
Rebound spiking/bursting 374 The following paragraphs show an example of "rebound activity" with the eLIF model (Figure 8), as well as 375 details about the conditions leading to rebound activity for the AdExp and mAdExp models.  Table  2 in Appendix A.6 for detailed parameters. This can also be shown mathematically for I.a and II.a by looking at the eigenvector associated to the lowest 388 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 391 II.a.

392
mAdExp The new rebound bursting behavior is associated to a positive divergence of the V -nullcline (cf.

393
Equation 10). 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 396 with the constraints of the AdExp and eLIF models:

397
• 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 398 neuronal type with a ≤ 0 (note that for small values of a, the sag, though technically present, can be 399 neglected for all practical purposes),

400
• the condition for no overshoot due to energy dynamics is E u ≥ E 0 (necessary for the ).      (Nov. 1993). "A New Approach to Solving the Cubic: Cardan's Solution Revealed". en. In: