## Abstract

Homeostatic intrinsic plasticity (HIP) is a ubiquitous cellular mechanism regulating neuronal activity, cardinal for the proper functioning of nervous systems. In invertebrates, HIP is critical for orchestrating stereotyped activity patterns. The functional impact of HIP remains more obscure in vertebrate networks, where higher order cognitive processes rely on complex neural dynamics. The hypothesis has emerged that HIP might control the complexity of activity dynamics in recurrent networks, with important computational consequences. However, conflicting results about the causal relationships between cellular HIP, network dynamics, and computational performance have arisen from machine-learning studies. Here, we assess how cellular HIP effects translate into collective dynamics and computational properties in biological recurrent networks. We develop a realistic multiscale model including a generic HIP rule regulating the neuronal threshold with actual molecular signaling pathways kinetics, Dale's principle, sparse connectivity, synaptic balance, and Hebbian synaptic plasticity (SP). Dynamic mean-field analysis and simulations unravel that HIP sets a working point at which inputs are transduced by large derivative ranges of the transfer function. This cellular mechanism ensures increased network dynamics complexity, robust balance with SP at the edge of chaos, and improved input separability. Although critically dependent upon balanced excitatory and inhibitory drives, these effects display striking robustness to changes in network architecture, learning rates, and input features. Thus, the mechanism we unveil might represent a ubiquitous cellular basis for complex dynamics in neural networks. Understanding this robustness is an important challenge to unraveling principles underlying self-organization around criticality in biological recurrent neural networks.

## Introduction

Homeostatic intrinsic plasticity (HIP) stabilizes nervous systems functions against environmental or pathological perturbations, through a negative feedback of neuronal firing on intrinsic excitability (Schulz, 2006). For instance, HIP maintains stereotyped network activity patterns in invertebrates (LeMasson et al., 1993; Marder and Goaillard, 2006). In vertebrates, HIP was proposed to regulate the average level of neuronal activity (Maffei and Turrigiano, 2008; Remme and Wadman, 2012). However, perceptual, motor, or cognitive processes rely on neural dynamics; its complexity largely exceeds simple features such as target levels or stereotyped patterns of activity (Vogels et al., 2005). HIP was thus suggested to regulate the complexity of activity dynamics in vertebrate recurrent networks (e.g., neocortical columns, the hippocampus; Lazar et al., 2007, 2009; Siri et al., 2007, 2008; Wardermann and Steil, 2007; Markovíc and Gros, 2012) with possibly important functional consequences. However, the situation remains largely obscure for several reasons.

First, the effect of HIP on network dynamics is unclear as it was found to either increase or decrease the complexity of the dynamics, depending on conditions (Lazar et al., 2007; Wardermann and Steil, 2007; Markovíc and Gros, 2012). Second, HIP was suggested as a mechanism preventing the pathological saturation of firing frequency and synchronization that is typical of Hebbian synaptic plasticity (SP; Turrigiano, 1999; Abraham, 2008; Watt and Desai, 2010). Yet, the proposition (Siri et al., 2007, 2008) that HIP may oppose these pathologies, by counteracting the underlying decrease in the complexity of the dynamics, has remained unexplored. Third, recurrent networks display improved recognition performances at the critical boundary between irregular and regular dynamics (the “edge of chaos”; Bertschinger and Natschläger, 2004; Siri et al., 2007, 2008). However, in networks with HIP, critical dynamics does not imply improved performance (Lazar et al., 2007), which generally arises at subcritical (regular) regimes (Steil, 2007; Wardermann and Steil, 2007; Lazar et al., 2009), urging for a better understanding of these relations (Lazar et al., 2007, 2009; Steil, 2007; Schrauwen et al., 2008).

Here, we focus on the unsolved issue of the transmission of HIP effects along scales, i.e., how can one explain HIP effects at the global network level, based on its original action, which actually operates at the level of individual neurons? To evaluate the biological relevance of this transmission, we model cardinal architectural and plasticity features characterizing vertebrate recurrent networks. In particular, we derive a generic biophysical HIP rule modeling the ubiquitous activity-dependent kinase/phosphatase regulation of sub-threshold conductances (Delord et al., 2007), which sets the current threshold for firing (Naudé et al., 2012), in contrast to previous models based on rules that are phenomenological or derived to optimize information transfer. Using a dynamical mean-field approach and numerical simulations, we unravel a new cellular mechanism of HIP action, the transduction of inputs by large derivative ranges of the transfer function, which ensures increased dynamics complexity, robust interaction with SP at the edge of chaos, and improved pattern separability at the network scale.

## Materials and Methods

##### Neuronal dynamics.

We study the effects of HIP on the activity dynamics of random recurrent neural networks. Our neural network model includes *N* firing-rate neurons with discrete dynamics. To account for the slower timescale (approximately minutes to days) of synaptic and intrinsic plasticity, compared with neuronal activity (∼10 ms; Siri et al., 2007, 2008), we update the synaptic weights and firing thresholds at a timescale slower than that of neuronal activities. We call a “learning epoch” the period of τ successive updates of neuron activities separating two updates of the synaptic weights and thresholds. In real neurons, the induction of plasticity (through signaling pathways) can be quite fast (approximately seconds to minutes), but its expression (as effective receptor or voltage-gated channels changes) operates at lower timescales (approximately minutes to days). Here, a learning epoch (the time step for learning) corresponds to 5.10^{3} or 10^{4} activity time steps of 10 ms, i.e., ∼ 1 min, encompassing the timescale of induction as well as a minimal estimate of the timescale of expression of plastic processes. In simulations, learning typically operates within about *T* = 1000 epochs, i.e., ∼15 h. We could have considered higher numbers of time steps for the learning epoch, to encompass larger timescales for the expression of plastic processes. However, this would have dramatically slowed down simulations, which were already quite demanding. Actually, the present results are independent of the exact timescale of the learning epoch, as long as it is sufficiently slower, compared with the timescale for neuronal activation, so that neuronal dynamics can converge toward its attractor within each learning epoch (which was the case here). In the following, *t* denotes the update step of neuronal activation, and *T* the update step of synaptic weights and neuronal thresholds (i.e., the learning epoch). Let *x*_{i}^{(T)}(*t*) be the firing frequency of neuron *i* at time, during learning epoch *T*, ranging between 0 (quiescence) and 1 (saturation). The dynamics of the model is given by the following:
where *f*(*u*) = (1 + tanh (*Gu*))/2 is a sigmoidal transfer function, *G*, the gain of the transfer function, and θ_{i}^{(T)} the neuronal threshold of neuron *i* during the epoch *T*. The quantity *u*_{i}^{(T)}(*t*) is called the local field:
i.e., the weighted sum of synaptic inputs due to recurrent connectivity in the network (where *w*_{ij}^{(T)} is the synaptic weight from the presynaptic neuron *j* to the postsynaptic neuron *i*), plus ξ* _{i}*, a constant external input (the vector ξ = {ξ

*} is the input pattern applied to the network). The timescale separation between neuron activity and plasticity updates allows the analysis of network dynamics based on numerical estimation of the largest Lyapunov exponent*

_{i}*L*

^{(T)}. The activity

*x*

_{i}

^{(T)}(

*t*) of each neuron

*i*during each learning epoch

*T*is thus computed for each time

*t*using constant values for the synaptic weights

*w*

_{ij}

^{(T)}and firing threshold θ

_{i}

^{(T)}. Neuronal activities, averaged over τ time steps, are used to update the synaptic weights and the neuronal thresholds for the next learning epoch, according to the plasticity rules (see below). The activity at the end of a learning epoch is then used to set the initial conditions of the next learning epoch

*x*

_{i}

^{(T+1)}(1) =

*x*

_{i}

^{(T)}(τ).

##### Network architecture.

The network connectivity is set at random using a biologically realistic setting with sparse connectivity and segregation of excitatory and inhibitory neurons (following Dale's principle; Siri et al., 2007). Each neuron is inhibitory (with probability *p _{I}*) or excitatory (

*p*= 1 −

_{E}*p*) and projects to

_{i}*N*=

_{c}*p*randomly chosen postsynaptic neurons with uniform probability, where

_{c}N*p*is the connection probability, independent of the polarity of the synapse. We chose

_{c}*p*and

_{I}*p*to reflect typical cortical architectural constraints. The initial synaptic weights

_{c}*w*

_{ij}

^{(1)}of connections are drawn from a gamma distribution, which ensures that excitatory (respectively, inhibitory) neurons only project to synapses with positive (respectively, negative) weights. If neuron

*j*is inhibitory,

*w*

_{ij}

^{(1)}follows Γ( − μ

*/*

_{w}*n*; σ

_{I}*/*

_{w}*n*) where Γ(

_{I}*m*,

*s*) is a gamma distribution with mean

*m*and SD

*s*, with

*n*=

_{I}*p*being the total number of synapses from an inhibitory neuron. Similarly, if neuron

_{I}p_{c}N*j*is excitatory, then

*w*

_{ij}

^{(1)}follows Γ( − μ

*; σ*

_{w}/n_{E}*) with*

_{w}/n_{E}*n*=

_{E}*p*. Normalization terms (

_{E}p_{c}N*n*and

_{E}*n*) guarantee that during the first epoch, the expected total excitation received by a postsynaptic neuron equals the expected total inhibition (balanced network). This initial global balance between excitation and inhibition can eventually be broken by SP. Note that with this setting, the initial weight matrix

_{I}*W*(1) is asymmetric, i.e.,

*w*(

_{ij}*1*) differs from

*w*(1) (with probability 1). Autapses were prohibited throughout learning epochs (

_{ji}*w*= 0 for all

_{ii}*T*).

##### SP.

We used a Hebbian SP rule to update synaptic weights at the end of each epoch (Siri et al., 2007, 2008), depending on the average firing frequencies of presynaptic and postsynaptic neurons, as follows:
where α is the Hebbian plasticity rate, and λϵ[0;1] is a passive forgetting coefficient for weights' evolution. In this equation *m*_{i}^{(T)} denotes the reduced average activity of neuron *i* during epoch *T* as follows:
where *d* = 0.1 is the threshold from which the neuron *i* is considered as active during epoch *T*, i.e., when its average activity is >10% of its maximal value. *H* () is the Heavyside function (*H*(*m*) = 0 if *m* < 0, 1 otherwise), and prevents heterosynaptic plasticity. The term *m*_{i}^{(T)} *m*_{j}^{(T)} *H*(*m*_{j}^{(T)}) represents the correlative principle of Hebb's rule: weights increase when both the presynaptic and postsynaptic neurons are active above *d*. The Heaviside function requires that the presynaptic neuron is active to induce synaptic weights' changes, preventing heterosynaptic depression and potentiation when both presynaptic and postsynaptic neurons are inactive. Finally, *s _{ij}* represents the sign of the synapse:

*s*= +1 for presynaptic excitatory neurons,

_{ij}*s*= −1 for inhibitory ones. The excitatory or inhibitory nature of the presynaptic neuron also affects the Hebbian plasticity rate α. Indeed, to ensure the global balance of excitation and inhibition in our setting, where excitatory neurons are three times more numerous than inhibitory ones, we used α

_{ij}*= α*

_{I}p_{I}*, where α*

_{E}p_{E}*(respectively, α*

_{E}*) stands for excitatory (inhibitory) SP rates. This global excitation/inhibition balance has dramatic consequence on the effects of intrinsic plasticity on network dynamics (see Results).*

_{I}##### HIP.

In the framework of firing rate coding, intrinsic plasticity affects neuronal excitability by modifying the transfer function, i.e., the relationship between the input and the resulting firing frequency; here, we consider the homeostatic plasticity of the maximal conductance of a subthreshold voltage-gated ionic current (Desai et al., 1999; Gibson et al., 2006; van Welie et al., 2006), i.e., an ionic current activating below the threshold for action potential (Naudé et al., 2012). This includes, e.g., the leak, persistent and slowly inactivating sodium, low-threshold calcium, or muscarinic, and slowly inactivating potassium currents, which affect the current threshold of the neuronal transfer function (Bekkers and Delaney, 2001; Brickley et al., 2001; Vervaeke et al., 2006; O'Leary et al., 2010). We consider an intrinsic learning rule preserving essential features of intrinsic plasticity: (1) HIP modifies the threshold θ_{i}^{(T)} between learning epochs, shifting the transfer function and modifying the input range to which the neuron responds, and (2) θ_{i}^{(T)} evolves as a function of the neuron mean firing frequency <*x*_{i}^{(T)}>, through the activation of intracellular biochemical signaling pathways.

In the model, the threshold of the function transfer is modified by HIP according to the following:
where θ̄ is the maximal threshold value and *F*_{i}^{(T)} the phosphorylated fraction of the plastic subthreshold conductance, assumed to constitute the functional pool of channels affecting excitability. Initially, all neuron thresholds are set to 0, i.e., *F*_{I}^{(0)} = − *A*_{θ}/*B*_{θ}. Parameters *A*_{θ} and *B*_{θ} are set so as to ensure a homeostatic regulation of the threshold by neuronal activity: an increase (respectively, a decrease) of the firing frequency leads to a decrease (respectively, an increase) of neuronal excitability. The phosphorylated fraction of the subthreshold plastic conductance *F*_{i}^{(T)}, which determines the threshold value, is regulated by activity-dependent kinases and phosphatases (Ganguly et al., 2000; Cudmore and Turrigiano, 2004; Zhang et al., 2004). We model these intracellular pathways using the activity-dependent kinase and phosphatase (aKP) framework (Delord et al., 2007), adapted to discrete time dynamics as follows:
where *K*_{i}^{(T)} and *P*_{i}^{(T)} represent mean activation levels of protein kinases and phosphatases K and P at learning epoch *T* (Delord et al., 2007). This biochemical pathway is activated by intracellular calcium concentration, as observed experimentally (Ganguly et al., 2000; Cudmore and Turrigiano, 2004; Zhang et al., 2004) as follows:
and
where *p* is the Hill constant, *k*_{max} and *p*_{max} are kinetic coefficients, *K _{p}* and

*K*are the phosphatase and kinase calcium half-activations, respectively.

_{K}*Ca*

_{i}

^{(T)}represents the mean concentration of the intracellular calcium during learning epoch

*T*. At the equilibrium state, intracellular calcium concentration linearly depends on the mean firing rate (Helmchen et al., 1996; Wang, 1998), which we simulate with the following: with basal concentration

*Ca*

_{0}(Wang, 1998; Delord et al., 2007). A non-null basal calcium concentration ensures non-null rate of plasticity for very low neuronal activities (Delord et al., 2007). This HIP model is sufficiently realistic to take into account essential dynamical properties of intrinsic plasticity, i.e., the gradation of threshold modifications and the activity-dependent time constant governing its dynamics (Delord et al., 2007).

##### Input patterns.

To assess the ability of the network to learn and discriminate inputs, we used *n*_{ξ} arbitrary input patterns ξ* _{k}*(

*k*= {1, …,

*n*

_{ξ}}), defined as

*i*is the neuron index,

*N*the total number of neurons, and

*f*

_{k}_{,1}and

*f*

_{k}_{,2}two trigonometric functions (see below). We used either

*n*

_{ξ}= 2, 3, or 4 input patterns to assess how separability behaves with increasing numbers of inputs. In a given condition, the network was trained by repeating the sequence formed by the

*n*

_{ξ}inputs at each successive learning epoch (e.g., ξ

_{1}, ξ

_{2}, ξ

_{1}, ξ

_{2}, … for

*n*

_{ξ}= 2). At the end of every learning epoch, synaptic weights and thresholds were modified according to SP and HIP rules, respectively, using the neuronal activities evoked by the specific input presented. Besides, at each learning epoch, the network was presented off-line with the

*n*

_{ξ}− 1 other input patterns. These responses were not used to update weights and thresholds, but only to compute the input pattern separability of the network. Specifically, the distance between the response of the network to two given inputs ξ

*and ξ*

_{k}*was computed, at each learning epoch*

_{I}*T*, as

*i*at epoch

*T*during presentation of the input pattern ξ. This distance was normalized by the distance between the corresponding inputs

*D–F*,

*S*was averaged over simulations run for 10 random realizations of the synaptic weight matrix and computed at the end of the simulation, once learning was achieved. The functions {

*f*

_{k,1},

*f*

_{k,2}} were chosen as {sin, cos}, {cos, sin}, {cos, cos}, and {sin, sin} for the four inputs ξ

_{1}, ξ

_{2}, ξ

_{3}, ξ

_{4}, respectively.

##### Numerical estimation of the Lyapunov exponent from numerical simulations.

Lyapunov exponents measure the local tendency of nearby trajectories living on an attractor of the dynamics: (1) to get closer, on average, in some directions of the phase space, corresponding to a negative Lyapunov exponent and (2) or to move apart in other directions, corresponding to a positive Lyapunov exponent. There are as many Lyapunov exponents as the phase space dimensions (*N* in our model). The sum of all Lyapunov exponents gives the local average volume contraction in the tangent space of an attractor. This sum is always negative from the definition of an attractor. When all Lyapunov exponents are negative the attractor is a fixed point, or, in the discrete time case, a finite union of points constituting a stable periodic orbit. If at least one exponent is positive, then the attractor is chaotic. The largest Lyapunov exponent gives an efficient and rapid way of characterizing the system dynamics: fixed point or periodic orbit when negative, chaos when positive. When it is zero the attractor has some neutral direction neither expanding nor contracting. This is the case, e.g., when the attractor is a closed invariant curve (typically a distorted circle where dynamics acts as a rotation with an irrational winding number), or an invariant torus.

For a given neural network, both synaptic weights and neuronal thresholds evolve through learning epochs, according to Hebbian SP and HIP rules, respectively, from their initial values (random synaptic weight matrix, null thresholds; see above). The largest Lyapunov exponent, *L*_{1}^{(T)}, was estimated numerically from simulations of the neural network models with a QR-based method (Von Bremen et al., 1997) at each learning epoch *T*, i.e., using the synaptic weights *w*_{ij}^{(T)} and the neuronal thresholds θ_{i}^{(T)} values. Employing these values as fixed parameters for the learning epoch *T*, a random initial condition of the network activity was drawn and the network dynamics was first iterated for 2000 activity update steps to allow convergence to the (fixed-point, limit-cycle, quasiperiodic, or chaotic) attractor. The largest Lyapunov exponent was then computed over successive periods of 1000 activity update steps, until reaching an absolute convergence criterion of 10^{−3} between successive Lyapunov exponents. This procedure was repeated at every learning epoch *T*, to compute the learning-induced evolution of the largest Lyapunov exponent, *L*_{1}^{(T)}. To ensure robustness of this computation, *L*_{1}^{(T)} was averaged over 10 realizations of the random synaptic weight matrix (see Results). When relevant, the steady-state value of the largest Lyapunov exponent, *L*_{1}^{(∞)}, was estimated by considering the value of *L*_{1}^{(T)} at large learning epochs (*T* = 2000). Note that both *L*_{1}^{(T)} and *L*_{1}^{(∞)} differ from the theoretical largest Lyapunov exponent, *L*_{1}^{MF}, for two reasons. First, whereas *L*_{1}^{(T)} and *L*_{1}^{(∞)} are obtained from the standard biological model, *L*_{1}^{MF} is derived from a simplified network model characterized by synaptic homogeneity, i.e., devoid of Dale's principle and sparse connectivity. Second, *L*_{1}^{(T)} and *L*_{1}^{(∞)} are estimated numerically from simulations of finite-sized networks, while *L*_{1}^{MF} is computed through an analytical expression obtained from the dynamical mean-field theory valid at the thermodynamical limit (*N* → + ∞). Thus, a direct quantitative comparison of *L*_{1}^{(T)} or *L*_{1}^{(∞)} with *L*_{1}^{MF} is not relevant. However, qualitative forms of agreement (e.g., similar parameters dependence) indicate the ability of the theory to account for the mechanisms underlying HIP effects.

##### Standard network model.

In the present study, we define the standard model as a recurrent network endowed with balanced excitatory/inhibitory drive, Dale's principle, sparse connectivity, and no plasticity. In the different sections, either HIP, SP, or both forms of plasticity are added to study their effects on emerging properties. Unless stated, the parameters used are *N* = 200, *G* = 5, τ = 10^{4}, *p _{i}* = 0.25,

*p*= 0.15, μ

_{c}*= 50, σ*

_{W}*= 1,*

_{W}*A*

_{θ}= 1.1,

*B*

_{θ}= −3.635,

*Ca*

_{0}= 0.1,

*p*= 4,

*p*

_{max}=

*k*

_{max}= 10

^{−3},

*K*= 1/3,

_{p}*K*= 2/3, θ̄ = 20.

_{K}## Results

### Mean-field analysis of HIP effects

In this first section, we study the effects of the sole HIP, i.e., in the absence of Hebbian SP, to assess the mechanism by which it affects neuronal dynamics at the global scale of the network. To this aim, we use a dynamic mean-field theory formerly developed for this type of network (Cessac, 1995; Doyon et al., 1995). This method allows the characterization of the dynamics averaged over the synaptic weight distribution. Here, in absence of SP, we suppose that synaptic weights are independent and identically distributed, as usually considered in dynamic mean-field theories with random interactions (Sompolinsky et al., 1988; Arous and Guionnet, 1995; Faugeras et al., 2009). In addition, HIP is a local plasticity rule that does not affect the correlation between synaptic weights. Hence, it does not correlate initially uncorrelated weights and it is possible to extend the results obtained in (Moynot and Samuelides, 2002). We denote *E*( ) the expectation over synaptic weight matrix realizations. Our mean-field analysis assumes synaptic homogeneity, i.e., we considered a fully connected network composed of a single population of neurons, in the thermodynamical limit (*N* → + ∞). In this case, the local field (i.e., the input − the threshold) of neurons, *u _{i}*, are Gaussian (Cessac, 1995; Moynot and Samuelides, 2002), with the following expectation:
where the neuronal threshold expectation is obtained from Equation 5 as follows:
The expectation of the threshold is itself altered by modifications of the plastic variable

*F*, which evolves according to (Eq. 6, see Materials and Methods) the following: where the expectation of neuronal activity is given by (Moynot and Samuelides, 2002) the following: with and where

*V*

^{(}

^{T}^{)}is the variance of the local field. Let

*V*

_{i}

^{(T)}be the variance of

*u*

_{i}

^{(T)}at the epoch

*T*. The mean-field theory establishes that

*V*

_{i}

^{(T)}is given, in the thermodynamic limit, by the following: A specific case holds when μ

_{j}

^{(T)}does not depend on

*j*(the input is independent of

*j*). Then,

*V*

^{(T)}=

_{j}

^{(T)}(Eq. 9) which in turn determines its variance

*V*

^{(T)}(Eq. 14). Second, both μ

_{i}

^{(T)}and

*V*

^{(T)}affect the complexity of the network dynamics (Eq. 15). Third, these variables act through the derivative of the transfer function of neurons,

*f*′, which determines neurons' sensitivity to their inputs (Eq. 15). According to Equation 15, decreasing the absolute value of the expectation μ

_{i}

^{(T)}should reduce the network's dynamics complexity. However, it is not straightforward to predict from the mere equation of how the combined effect of μ

_{i}

^{(T)}and

*V*

^{(T)}of local fields should affect this complexity. To evaluate this effect, we computed the maximum Lyapunov exponent predicted by the dynamical mean-field theory. That is, we iterated Equations 9–14 until they reached a stable fixed point, corresponding to a steady-state distribution of the network activity (Cessac, 1995). Then, from Equation 15, we were able to compute the predicted maximal Lyapunov exponent,

*L*

_{1}

^{MF}, as the steady-state value of

*L*

_{1}

^{MF(T)}. Specifically, we computed

*L*

_{1}

^{MF}for a large range of the external input ξ (i.e., [ − θ̄, θ̄], i.e., [ − 20, 20]) in the presence or absence of HIP, to evaluate the robustness of HIP on the dynamics of the network.

We observed that, in the absence of HIP, *L*_{1}^{MF} reaches negative values when increasing the amplitude of the input ξ. A negative maximal exponent corresponds to a regular dynamics (stable fixed-point or periodic orbit; Fig. 1*A*, red). Actually, without the homeostatic compensation provided by HIP, large inputs (Fig. 1*B*, vertical red band) are transformed through the transfer function (plain black curve) into large outputs (*x _{i}* =

*f*(

*u*) ≈ 1; horizontal red stripe). However, the derivative of the transfer function (dashed black curve) at large inputs is small (

_{i}*f′*(

*u*) ≈ 0). Therefore, fluctuations of the input (vertical red band) have little impact on postsynaptic neurons because they are contracted into small output fluctuations (horizontal red band). Hence, fluctuations of activity are dampened in the network (low variance of the activity; Fig. 1

_{i}*C*, red) and regular dynamics are promoted. In the dynamic mean-field theory, Equation 15 confirms that this reduced derivative is indeed the factor decreasing

*L*

_{1}

^{MF}, bringing the network's activity into more regular regime, where dynamics complexity is lower.

In contrast, in the presence of HIP, *L*_{1}^{MF} was positive, i.e., dynamics was chaotic, and virtually constant, for a large range of input amplitude (of the order of [ − θ̄, θ̄], the range of possible threshold modifications; Fig. 1*A*, yellow). In this case, because of the homeostatic compensation provided by HIP, even large inputs are compensated by a shift of the threshold, resulting in a null expectation of the local field distribution (Fig. 1*B*, vertical yellow band). As a consequence, the output is balanced (*x _{i}* =

*f*(

*u*) ≈ 0.5; horizontal yellow band), and the derivative is maximal (

_{i}*f′*(

*u*) ≈

_{i}*G*/2). Therefore, fluctuations of the input are amplified in the output and thus maintained in the network. The activity presents a large variance (Fig. 1

*C*, yellow), which should favor complex dynamics. Again, the mean-field theory we have developed evidences that the maximal derivative observed in the presence of HIP is indeed the factor increasing

*L*

_{1}

^{MF}, bringing the network's dynamics into a chaotic regime.

Thus, the dynamic mean-field theory suggests that within the range of possible threshold plastic modifications [ − θ̄,θ̄] HIP (1) exerts a compensation to external inputs that is robust to their exact amplitude and (2) guarantees complex dynamics operating through a mechanism relying on a large derivative of the transfer function, and translating as large fluctuations of neuronal activity. In the following sections, we compare the predictions of the dynamic mean-field theory, which are exact in the thermodynamical limit *N* → + ∞, to simulations in finite sized models.

### Cellular effects of HIP

Here, we first checked that in a single neuron model endowed with HIP, homeostatic compensation of the local field indeed maintains the activity variability characterizing the mechanism we have unveiled. To test this proposal, we consider a single plastic neuron receiving inputs from a simulated standard recurrent model network (i.e., Dale's principle, sparse connectivity, balanced excitatory/inhibitory drive, no plasticity, see Materials and Methods), with chaotic dynamics. We observed that, starting from an initial state where the neuron was saturated due to a particular realization of its incoming synaptic weights (Fig. 2*A*, upper histogram, firing rate as a function of the time *t* within learning epoch *T* = 1), the presence of HIP was able to restore a large variability of activity (Fig. 2*A*, lower histogram, learning epoch *T* = 250). Actually, when the activity of the plastic neuron was initially large, HIP caused a progressive increase of the firing threshold, leading to a decrease of the time-averaged activity toward a median value (〈*x*(*t*)〉^{(T)} = 0.5; Fig. 2*B*, upper trace). Symmetrically, when the initial activity was low, HIP decreased the firing threshold, balancing activity (Fig. 2*B*, lower trace). In each case, HIP adapted the threshold and canceled the neuron's local field so that the neuron exploited the linear part of its transfer function (Fig. 2*C*), raising firing variability (Fig. 2*D*).

### Network effects of HIP

We next considered the effect of HIP on the standard recurrent network model, addressing the joint effect of HIP and SP in the next section. This standard model was designed as a biologically realistic recurrent network with sparse connectivity and segregation between excitatory and inhibitory neurons (Dale's principle). When no plasticity was present, we observed that, at the end of the simulations (*T* = 2000), the dynamics of the network was regular, e.g., a periodic orbit, as illustrated by the firing rate histogram of a individual neuron (Fig. 3*A*, upper, the bar indicates one period of the oscillation) and the raster of network activity (Fig. 3*B*, upper). In contrast, the activity in the presence of HIP led to chaotic dynamics (Fig. 3*A*, lower; *B*, lower).

Using this network and a constant external input applied to all neurons, we tested the theoretical predictions made by the dynamical mean-field approach developed above. We performed these simulations with *N* = 200 and computed numerically the maximal Lyapunov exponent (see Materials and Methods), averaging its value over 10 realizations of the synaptic weight matrix, for each external input amplitude tested. In the following, we denote *L*_{1}^{(T)} the maximal Lyapunov exponent computed in the simulations at the learning epoch *T*, and *L*_{1}^{(∞)} its steady-state value at large simulation times (*T* = 2000). Since the mean-field theory is exact for a fully connected model in the limit *N* → + ∞, one cannot expect a perfect match between the Lyapunov exponent derived from the theory (*L*_{1}^{MF}) and those estimated in the simulations (*L*_{1}^{(T)} or *L*_{1}^{(∞)}). However, qualitative agreement indicates the heuristic capacity of the dynamical mean-field theory.

We observed that the steady-state largest Lyapunov exponent *L*_{1}^{(∞)} was higher with HIP (Fig. 3*C*, yellow) than without (red). On average, *L*_{1}^{(T)} increased over time in the presence of HIP and stabilized at positive values, indicating a chaotic behavior (Fig. 3*D*, yellow). Moreover, HIP led the network to chaotic dynamics independently of the initial network dynamics (chaotic (*L*_{1}^{(T)} > 0) or not (*L*_{1}^{(T)} < 0); Figure 3*D*, black traces), which depended on the random realizations of synaptic weights (that possibly induced strong variations of individual neuronal local fields). The resulting chaotic behavior lives on a strange attractor (Eckmann and Ruelle, 1985). It is possible to visualize attractors by plotting the average network activity during consecutive simulation times (*t* and *t* + 1), using the Takens method (i.e., with time delay 1 here; Takens, 1981; Doyon et al., 1994). Such a representation in the plane {*x̄*(*t*),*x̄*(*t* + 1)} reveals the topological structure of the attractor underlying network dynamics. For instance, a single point reveals that dynamics is stuck at a fixed-point attractor, separate points indicate a periodic orbit attractor within which the system migrates, and a closed 1D dense variety (curve) indicates a quasiperiodic attractor, while a cloud of points reveals a chaotic attractor, as found here (Fig. 3*E*). Note that HIP increased the dynamics complexity of biological networks regardless of the value of the external input (i.e., even with ξ = 0; Fig. 3*C*), whereas the dynamical mean-field approach predicted identical *L*_{1}^{MF} values at ξ = 0 (Fig. 1*A*). Despite this discrepancy, attributable to the different hypotheses considered (sparse vs full connectivity, Dale's principle or not, and finite or infinite size), we found a good qualitative agreement between the theory and the simulations. Indeed, the similar dependence upon input amplitude and HIP presence demonstrated that the theory correctly accounts for the mechanisms underlying HIP effects.

For the range of parameters we used, we found that, independently of (1) the amplitude of the input current considered and (2) the random realization of the weight matrix, HIP was able to strongly decrease the saturation of neuronal activities (Fig. 3*F*, yellow), even when an important part of the network was initially saturated (red). The results obtained with network models (i.e., in the present and following sections) were robust to large variations of the network size (data not shown). Thus, HIP was able to maintain dynamical and chaotic fluctuations in the network and provided large variability of neuronal activity, confirming the mechanism unveiled by the dynamical mean-field approach.

### HIP and SP interaction in the recurrent network

We addressed the important issue of whether HIP could prevent runaway effects of SP, such as pathological levels of firing frequency and synchronization resulting from the positive feedback between neuronal activity and synaptic connectivity induced by SP (Abraham, 2008; Watt and Desai, 2010). In recurrent network models, these “runaway” effects manifest as a decrease of dynamics complexity toward regular dynamics (i.e., from positive to negative maximal Lyapunov exponents (Siri et al., 2007, 2008). It was proposed that HIP might oppose this dead-end evolution (Siri et al., 2007, 2008). However, this issue has remained unexplored hitherto. Unfortunately, the issue of HIP and SP interaction cannot be addressed in the realm of dynamic mean-field theory. Indeed, SP creates correlations between synaptic weights, whereas mean-field theory assumes statistical independence between synaptic weights. To our knowledge, no current theoretical tool yet exists to handle this case, so we relied on numerical simulations to investigate these issues.

We typically found that the dynamics of the network corresponded to a fixed point at the end of the simulations (*T* = 2000) when SP was present, as illustrated for the activity of an individual neuron (Fig. 4*A*, upper) and of the network (Fig. 4*B*, upper). In contrast, the activity in presence of HIP in addition to SP led to complex patterns of activity (Fig. 4*A*, lower firing rate histogram; *B*, lower raster). To quantify this apparent increase in complexity, we computed the largest Lyapunov exponent in standard recurrent networks endowed with SP or SP and HIP (see Materials and Methods). We followed the evolution of the maximal Lyapunov exponent (averaged over 10 realizations of the network) as a function of the learning epoch *T*. We first observed that in the presence of the sole SP, the steady-state value of the largest Lyapunov exponent, *L*_{1}^{(∞)}, was negative (Fig. 4*C*, purple). This occurred whatever the value of α/(1 − λ), which quantifies the magnitude of SP (Siri et al., 2007, 2008). Note that the largest Lyapunov exponent increased with α/(1 − λ). This arose because at lower values of this ratio, passive forgetting dominates the Hebbian term (see Eq. 3), so that synaptic weights strongly decrease and neurons tend to be disconnected. In such a situation, the network encountered the most static situation possible, i.e., the lowest largest Lyapunov exponent. Above the range illustrated (α/(1 − λ) < 200), the Hebbian rate became so large that the network was unstable, i.e., the synaptic matrix did not converge anymore to a steady state and the Lyapunov exponent exhibited oscillations. We therefore limited our exploration to this range. Even when initial conditions gave rise to an initially positive largest Lyapunov exponent *L*_{1}^{(1)}, the SP rule led to a decrease of *L*_{1}^{(T)} to negative steady-state values *L*_{1}^{(∞)} (Fig. 4*D*, purple). Accordingly, the network attractor observed at long learning times was a fixed point (Fig. 4*E*, purple). Moreover, our simulations indicated that this fixed point corresponds to a divergence of neuronal activities toward saturation (Fig. 4*F*, purple histogram), an illustration of the “runaway” effects induced by SP.

When HIP was added to SP, the largest Lyapunov exponent *L*_{1}^{(T)} converged toward steady-state values (*L*_{1}^{(∞)}) close to zero or even positive in a large domain of synaptic rates tested (Fig. 4*C*; *D*, green). Hence, the presence of HIP led the networks to chaotic dynamics (*L*_{1}^{∞} > 0) or to the so-called edge of chaos (*L*_{1}^{∞} ≈ 0; Fig. 4*E*, green dots). To evaluate the influence of SP on the effects of HIP in the network model, we ran simulations, in which we varied systematically the rate of SP, α/(1 − λ). For a large range of this parameter, network dynamics reached the boundary between irregular and regular dynamics (i.e., near-zero largest Lyapunov exponent; Fig. 4*C*, green) and the system maintained itself in this regime (Fig. 4*D*, green). Finally, in these networks endowed with both SP and HIP, neuronal activities were closer to 0.5, with large variability (Fig. 4*F*, green). Hence, even in the presence of SP, the effects of HIP rely on the mechanism unveiled in the dynamical mean-field theory. Moreover, our results support the hypothesis of HIP as an effective regulator of the runaway effects of SP, since it counteracts the divergence of neuronal activities and the subsequent decrease in dynamics complexity (Siri et al., 2007, 2008; Abraham, 2008; Watt and Desai, 2010).

### Biological robustness of HIP effects

We assessed the effect of several important biological constraints, related to network architecture and SP balance, on the ability of HIP to control the network dynamics through its homeostatic influence.

We first compared different network architectures to investigate the impact of Dale's principle (i.e., the separation of excitatory-only and inhibitory-only neurons) and of the sparseness of the standard recurrent network (see Materials and Methods). We built the surface map of the average steady-state largest Lyapunov exponent at the end of simulations, *L*_{1}^{(∞)}, as a function of the rates of SP (α/(1 − λ)) and HIP (θ̄), for the standard network model (synaptic balance, Dale's principle, and sparse connectivity). This surface was essentially composed of a flat region that spanned large ranges of the rates space, at *L*_{1}^{(∞)} values close to 0 (Fig. 5*A*). Periodic orbit (PO), critical (C; quasiperiodic, i.e., marginally stable), and chaotic dynamics all occupied significant domains in this region, while fixed-point (FP) dynamics was confined to its borders, at low plasticity rates (Fig. 5*B*). This illustrates the strong robustness of HIP effects on dynamics complexity to large variations of plastic rates. Discarding both Dale's principle and sparse connectivity in the network did not qualitatively affect the planar shape of the *L*_{1}^{(∞)} surface, indicating that the robustness of HIP effects extended to the presence or absence of these two important architectural determinants of recurrent networks (compare Fig. 5*A*,*B*).

We then investigated whether the balance between excitatory and inhibitory synaptic drive was important for the action of HIP. In the standard network model (see Materials and Methods), the parameter α, which determines the rate for SP, is set to different values for excitatory and inhibitory synapses to compensate for the larger proportion of excitatory neurons *p _{E}*α

*=*

_{E}*p*α

_{I}*). To destabilize this balance in these networks, we used identical values of the learning rate α for excitatory and inhibitory synapses, i.e., we set α*

_{I}*= α*

_{E}*= α. In this case, the excitatory drive is three times larger than that of the inhibition. Contrarily to balanced networks, unbalanced networks failed to stabilize at chaotic dynamical regimes or at the edge of chaos, systematically displaying fixed-point dynamics (FP; Fig. 5 compare*

_{I}*C*with

*A*,

*B*). As a general rule, neuronal activities were largely saturated so that the variability of neuronal activity remained low. SP caused a divergence of local fields toward positive values that could not be compensated by HIP (Fig. 5

*D*, bluish distribution), whereas this was possible in balanced networks (lime distribution).

Hence, these simulations specify some of the conditions necessary for a proper regulation of dynamical regime by HIP in recurrent networks. We pinpoint the negligible influence of static architectural constraints (i.e., Dale's principle and sparse connectivity), suggesting that HIP could serve a similar regulatory role in very different network types (e.g., cortical or not) in vertebrates. In contrast, HIP requires excitatory SP to be balanced by its inhibitory counterpart, which seems a ubiquitous property of vertebrate networks (Shadlen and Newsome, 1994; Tsodyks and Sejnowski, 1995; Gaiarsa et al., 2002). Hence, our modeling study stresses the necessary synergy between intrinsic and synaptic homeostatic processes to generate complex activity dynamics. Beyond that important requirement, our results indicate that the maintenance of complex dynamics by HIP is a very robust property that neither depends on (1) particular architectural constraints, (2) the presence of SP, and (3) specific rates of HIP and SP.

### Functional impact of HIP on pattern recognition

Previous studies have shown that computational performance on the recognition of static or time-varying inputs is improved when recurrent networks operate at the edge of chaos (Bertschinger and Natschläger, 2004; Siri et al., 2007, 2008). Machine-learning studies have assessed whether HIP may improve the recognition of time-varying stimuli by stabilizing recurrent network dynamics complexity at this border. However, in these studies, critical dynamics and improved performance have proved rather disjoint (see Introduction; Lazar et al., 2007, 2009; Steil, 2007; Wardermann and Steil, 2007). In the present model embedded with a more biologically realistic rule based on the actual enzymatic signaling pathways responsible for intrinsic plasticity in real neurons, we evaluated the possible functional impact of HIP in learning to discriminate static input patterns. The results obtained were robust to large variations in the scaling or spatial structure of input patterns. To readout the network behavior in a very general manner (i.e., without training any particular downstream output layer as in reservoir computing), we simply considered the ability to separate input patterns, based on its activity. To do so, we computed input separability, *S*, the network, and time-averaged distance between individual neuronal activity in response to pairs of arbitrary input patterns {ξ* _{k}*, ξ

*} (see Materials and Methods).*

_{I}In the absence of HIP, input separability increased during learning, as exemplified in Figure 6*A* (three realizations of the initial synaptic weights, purple traces), being maximal when *L*_{1}^{(T)} was close to zero (Fig. 4*D*, purple trace, i.e., around *T* ≈ 500). At this edge of chaos, the network lost its stability with respect to small perturbations, exhibiting strong nonlinear responses to the presence of inputs, and was therefore more sensitive to their presentation. As the learning epoch *T* increased, input separability decreased after this transient behavior and stabilized below its maximal value in the present balanced network (Fig. 6*A*, purple traces). However, inspection of the reconstructed attractors of the dynamics at large *T* values, in response to the two input patterns showed that the average activity level was comparable (Fig. 6*B*; i.e., two fixed-point attractors). In the presence of HIP (in addition to SP), input separability was higher than in the sole presence of SP, especially at long learning times (Fig. 6*A*, green traces). The reconstructed attractors of the dynamics occupied a larger region in the phase space and were more discriminable from each other when HIP was present (Fig. 6*C*; a quasiperiodic attractor (input ξ_{1}, light green) and a periodic attractor (input ξ_{2}, dark green)). Therefore, input separability was enhanced in the presence of HIP.

To better assess this computational effect of HIP, we evaluated separability in different learning conditions, i.e., for increasing values of θ̄ (the magnitude of HIP) and when teaching the network with *n*_{ξ} = 2, 3, or 4 input patterns. Consistent with our initial findings (Fig. 6*A–C*), we found that separability increased with the magnitude of HIP, compared with when SP only was present (Fig. 6*D*). Moreover, whereas separability was independent of the number of input patterns in the absence of HIP (Fig. 6*E*, purple), the enhanced separability due to HIP was remarkably increased when larger numbers of input patterns were presented to train the network (Fig. 6*E*, green).

Finally, to address the issue of how input separability relates to dynamics complexity in the biologically realistic networks considered in the present study, we systematically plotted separability as a function of the maximal Lyapunov exponent *L*_{1}^{(∞)} evaluated at the end of simulations (i.e., after learning), for all the conditions and network realizations simulated (Fig. 6*F*, increasing HIP magnitude: graded color tints from purple to green; *n*_{ξ} = 2 (circles), 3 (squares), or 4 (triangles)). We found that increasing values of θ̄ led to more complex dynamics (larger *L*_{1}^{(∞)} values) after learning, up to the edge of chaos (*L*_{1}^{(∞)} ≈ 0; no strongly supracritical dynamics was found, in contrast to simulations run in the absence of external inputs, where *L*_{1}^{(∞)} could rise up to ∼0.1; Fig. 5*B*). Linear regression of the data (Fig. 6*F*) indicated that, as an overall trend, input separability increased with the complexity of the dynamics (as a larger magnitude of HIP was considered). Thus, separability was enhanced as the complexity of dynamics approached criticality. In comparison, in realizations of the network devoid of any form of plasticity (Fig. 6*F*, black dots), we found that whereas the complexity of dynamics ranged in a large region encompassing criticality (*L*_{1}^{∞}ϵ [−0.2, 0.4]), the separability remained low (*s* ≈ 0.01), compared with when plasticity was present (*S* ϵ [0.08, 0.18]). Therefore, the dependence between the separability and the complexity of the dynamics was an acquired property of the network, which emerged from the interaction between synaptic and intrinsic and plasticity.

Together, these results indicate that criticality cannot be considered as a sufficient condition for improved computational performance (as found in some reservoir networks endowed with artificial plastic rules; Lazar et al., 2007). However, in networks endowed with realistic biophysical plastic processes, while SP alone increases separability at subcritical dynamics, separability significantly improves when dynamics approaches criticality, thanks to the combination of SP and HIP.

## Discussion

To our knowledge, we propose the first network model including a biologically realistic HIP rule, which is neither phenomenological (Lazar et al., 2006, 2007, 2009; Remme and Wadman, 2012) nor designed to maximize information transmission (Stemmler and Koch, 1999; Steil, 2007; Wardermann and Steil, 2007; Schrauwen et al., 2008; Markovíc and Gros, 2012). Here, HIP regulates the neuronal threshold through modifications of subthreshold voltage-gated conductances (Naudé et al., 2012), as frequently observed following HIP induction (Desai et al., 1999; Gibson et al., 2006; van Welie et al., 2006). Moreover, HIP modifications operate through a realistic description of signaling pathways' kinetics regulating intrinsic conductances in real neurons (Delord et al., 2007).

Functionally, it was initially proposed that HIP of the threshold is able to control the level of neuronal activity, avoiding neuronal saturation (Desai et al., 1999; Fig. 7*A*). Indeed, depending on its activity-dependence (e.g., calcium kinetics, kinase/phosphatase calcium dependence), a HIP rule can bring neurons to any arbitrary working point in the linear range of the transfer function. Besides, it was suggested that HIP actually maximizes information transmission (i.e., mutual information) between input and output distributions (Stemmler and Koch, 1999). Since, an artificial HIP rule was derived to implement this principle (Triesch, 2005) and improve reservoir learning (Steil, 2007; Schrauwen et al., 2008). In practice, such a rule adjusts the threshold so that low inputs are mapped evenly in the output distribution (Fig. 7*B*). In vertebrates, this principle makes sense at early sensory stages to transmit maximal information; however, in higher structures, neurons rather specifically respond to large input fluctuations, losing large amounts of information, to ensure neuronal selectivity (Ringach and Malone, 2007). HIP may either favor selection (Fig. 7*C*) or transmission in different neurons, depending on the HIP rule and transfer function (e.g., overall gain, symmetry).

The present study suggests that, beyond setting the overall level of activity or input processing, a prominent effect of HIP is to increase the complexity of activity dynamics, which is fundamental for higher order vertebrates' computations. This increase is consistent with that observed in reservoir networks, which share recurrent architecture with vertebrate networks (Lazar et al., 2007; Wardermann and Steil, 2007; Markovíc and Gros, 2012). In Lazar et al. (2007), HIP decreased complexity, but only when dynamics were spontaneous and SP present. Thus, our results strengthen the principle that increased dynamics complexity represents a general and robust (see below) consequence of HIP in biological recurrent networks. The mean-field approach we have devised indicates that increased complexity relies on a yet unexplored cellular effect of HIP: setting the neuronal working point so that input distributions are transduced by ranges of the transfer function with large derivative (Fig. 7*D*). Consequently, this principle may have important functional consequences, which we discuss now.

It has been suggested that HIP may compensate runaway effects of SP, i.e., destabilization to saturated activity levels arising from the positive feedback between activity and connectivity (Turrigiano, 1999; Watt and Desai, 2010). Previous studies from our group unraveled that SP-mediated pathological activity actually results from a general decrease of networks dynamics complexity due to SP (Siri et al., 2007, 2008). Here, we show that HIP increases dynamics complexity (through its effect on the transfer function derivative) and opposes the SP-mediated decrease, bringing network dynamics at the edge of chaos. This effect was only observed when SP preserved synaptic excitation/inhibition balance, an important determinant of network dynamics with uncorrelated synaptic weights (Shadlen and Newsome, 1994; Tsodyks and Sejnowski, 1995; van Vreeswijk and Sompolinsky, 1996). Our results thus suggest that HIP might act in concert with SP to control the dynamics complexity of recurrent networks, as was suggested for their activity level (Remme and Wadman, 2012). Such a synergistic regulation could arise from the coactivation of HIP and SP, through common signaling pathways (Misonou et al., 2004; van Welie et al., 2004; Fan et al., 2005; Karmarkar and Buonomano, 2006; Turrigiano, 2011). In two studies of the Triesch group, dynamics of the reservoir networks endowed with HIP and SP were clearly subcritical (Lazar et al., 2007, 2009). Such regular dynamics presumably emerges from the interaction between time-structured inputs and the spike timing-dependent SP rule considered in these studies. In contrast, here, inputs are static and the rules are not time dependent. The discrepancy might also possibly arise from the crude description of the inhibition (global k-winner-take-all dynamics) used in Lazar et al. (2007).

Moreover, with larger HIP magnitudes, performance increased as dynamics approaches the edge of chaos. In this region, the dynamics did not systematically exhibit exact criticality, but rather covered distinctive behaviors situated around this particular frontier, which included long periodic orbits and quasiperiodic or chaotic dynamics. The specific regime encountered depended on the input pattern, learning parameters, and initial conditions of the weights and thresholds. This result is consistent with previous studies of our group, where sensitivity to a single static input pattern was maximal at critical dynamics in recurrent networks endowed with SP (Siri et al., 2007, 2008). It also parallels the previous observation that in static reservoir recurrent networks (no plasticity, trained readout), computational performances on time-varying inputs are improved at criticality (Bertschinger and Natschläger, 2004). However, it was shown that HIP-endowed recurrent networks with critical dynamics do not systematically perform better than static networks (Lazar et al., 2007). Moreover, HIP noticeably improves reservoir networks performances that clearly exhibit subcritical dynamics (i.e., regular; Steil, 2007; Wardermann and Steil, 2007; Lazar et al., 2009). This discrepancy between the present study and these studies (Lazar et al., 2007, 2009; Steil, 2007; Wardermann and Steil, 2007) may stem from the different underlying operation, i.e., steady-state dynamics upon static inputs here (activity level, attractor nature) versus transient dynamics (i.e., “fading memory”) in response to time-varying inputs in reservoir networks. In particular, in reservoir models endowed with spike timing-dependent plasticity (Lazar et al., 2007, 2009), time-dependent SP and time-structured inputs presumably represent a common condition for both generating sub-critical dynamics (see above) and ensuring input recognition. Here, inputs are static and SP is not time dependent so that what is learned is the spatial structure of inputs. Thus, one might expect that best performances should be observed at very regular dynamics, e.g., fixed-point attractors. Performance is indeed enhanced at such regular regimes with SP. However, when HIP is added, performance culminates at criticality. This discrepancy might arise from the rules considered: actual kinetics of signaling pathways regulating excitability in real neurons (Delord et al., 2007; Naudé et al., 2012), here, versus phenomenological (Lazar et al., 2007, 2009) or optimized rules (Steil, 2007; Wardermann and Steil, 2007; Schrauwen et al., 2008). This divergence may require future efforts to be fully understood, as the mechanisms underlying improved performance have remained a constant question (Lazar et al., 2007, 2009; Steil, 2007; Schrauwen et al., 2008).

Here, we observed that input separability increases with the number of input patterns presented to train the network. However, for practical reasons, we restricted our study to limited pattern numbers (i.e., 4). In the future, it would be interesting to assess the effects of HIP on the memory capacity of networks for larger numbers of input patterns. Another appealing possibility would be to address possible improvements in the learning of complex activity patterns in FORCE networks that endow supervised network SP and a trained readout (Sussillo and Abbott, 2009), when incorporating HIP.

Our results indicate that the maintenance of complex dynamics by HIP in biological networks is a very robust property that neither depends on the size of the network, particular architectural constraints (Dale's principle and sparse connectivity), the presence of SP, specific learning rates, or the scaling and spatial patterning of inputs. In particular, the extreme robustness of HIP in bringing network dynamics at the edge of chaos in the presence of SP is stunning (Fig. 4*C*, 5*A*,*B*). Understanding this robustness is an important challenge for future studies that might unravel principles underlying self-organization around criticality in recurrent neural networks. Recently, it was suggested that HIP of both the threshold and gain of the transfer function might bring complex activity patterns in recurrent networks (Markovíc and Gros, 2012). Unfortunately, these behaviors were not extremely robust to changes in architectural or HIP parameters and were underlain by fast kinetics bringing rapid oscillations of excitability with no current experimental support. However, examining how the robustness we have observed with threshold HIP would extend to both threshold and gain HIP using slow kinetics in biological networks appears as an exciting possibility.

Together, the present results suggest that HIP effects may represent a ubiquitous determinant of biological networks across neuronal systems. Beyond the impressive robustness we have observed in recurrent architectures, the attractive perspective emerges that HIP processes share a general role in setting the working point of neural networks, independently of their specific function, e.g., the organization of rhythmic activity in invertebrates or the maintenance of complex dynamics for higher order computations and nontrivial behavior formation in vertebrates.

## Footnotes

This work was supported by the INRIA, ERC-NERVI number 227747, European Union Project FP7-269921 (BrainScales), and European Union Project FP7-318723 (MATHEMACS) to B.C.

The authors declare no competing financial interests.

- Correspondence should be addressed to Bruno Delord, Institut des Systèmes Intelligents et de Robotique, UMR CNRS 7222, Université Pierre et Marie Curie, 4, place Jussieu, 75005 Paris, France. bruno.delord{at}upmc.fr