Nonlinear multiplicative dendritic integration in neuron and network models

Neurons receive inputs from thousands of synapses distributed across dendritic trees of complex morphology. It is known that dendritic integration of excitatory and inhibitory synapses can be highly non-linear in reality and can heavily depend on the exact location and spatial arrangement of inhibitory and excitatory synapses on the dendrite. Despite this known fact, most neuron models used in artificial neural networks today still only describe the voltage potential of a single somatic compartment and assume a simple linear summation of all individual synaptic inputs. We here suggest a new biophysical motivated derivation of a single compartment model that integrates the non-linear effects of shunting inhibition, where an inhibitory input on the route of an excitatory input to the soma cancels or “shunts” the excitatory potential. In particular, our integration of non-linear dendritic processing into the neuron model follows a simple multiplicative rule, suggested recently by experiments, and allows for strict mathematical treatment of network effects. Using our new formulation, we further devised a spiking network model where inhibitory neurons act as global shunting gates, and show that the network exhibits persistent activity in a low firing regime.


INTRODUCTION
A hallmark of neural computations is the interaction of excitatory and inhibitory drives (Wilson and Cowan, 1972;Brunel, 2000;Herz et al., 2006;Vogels and Abbott, 2009). Dendritic integration of excitatory and inhibitory synaptic potentials is a complicated biophysical process and involves many non-linearities (Mainen and Sejnowski, 1996;Borg-Graham et al., 1998;London and Häusser, 2005). In particular, the position of the presynaptic inputs to a particular neuron on the dendritic tree is very important for signal integration and potential spike initiation (Huang et al., 2007). For instance, if an excitatory input occurred in the distal region of the dendritic tree (far away from the soma) it is possible that a simultaneous inhibitory input to the proximal region (near to the soma) could effectively "shunt" the excitatory pulse on its way to the soma. Accordingly, this effect is called shunting inhibition [Barrett and Crill, 1974;Blomfield, 1974;Rall, 2011, see also review in Koch (2004)]. On the other hand, if the position of the two synapses on the dendritic tree would be exchanged, a simultaneous input would have a markedly different outcome, with the excitatory pulse reaching the soma and potentially causing the generation of a spike.
While the theoretical basis for the effect of shunting inhibition has been laid out several decades ago by analyzing passive electric membrane properties of dendrites (Barrett and Crill, 1974;Blomfield, 1974;Rall, 2011;Koch, 2004), neural network models today still often rely on single compartment models (point models), such as (quadratic) integrate-and-fire model (Gerstner and Kistler, 2002;Izhikevich, 2006), and thus ignore this potentially important non-linear integration of synaptic inputs (McLaughlin et al., 2000;Wang, 2002;Vogels and Abbott, 2009;Rasch et al., 2011).
Acknowledging this inaccuracy of common network models already decades ago, it was shown that incorporating non-linear dendritic effects, such as shunting inhibition, into the mean-field equations of a network of neurons is indeed possible (Abbott, 1991a,b). These non-linear effects are potential important as they e.g., explain persistent activity of low firing rates which cannot be well explained by common models as they typically fire in the saturation regime of the input-output relation with unreasonably high firing rates (>100 Hz; Gold and Shadlen, 2001). However, by starting from the full cable equation, the equations in Abbott (1991b) for incorporating the non-linear dendritic effects are not of a simple mathematical form.
Indeed, it is known that the effect of shunting inhibition can be (at least during steady state) mathematically conceptualized simply as a "dirty multiplication" of excitatory and inhibitory input conductances (Koch et al., 1982;Koch, 2004) in contrast to the common form of linear summation of any two inputs in simple neuron models. Accordingly, a recent study, re-investigating shunting inhibition experimentally (Hao et al., 2009), found that indeed the integration of simultaneous excitatory post-synaptic potentials (EPSP) and inhibitory post-synaptic potentials (IPSP) could be well described in a multiplicative form: where κ is a factor determining the strength of the shunting effect that depends on the spatial arrangement of the synaptic inputs on the dendritic tree. Given the good experimental agreement of this description, a follow-up study tried to incorporate this abstract formulation of shunting inhibition into a simple single-compartment neuron model (Zhou et al., 2013), which the authors called the "DIF"-model (dendritic integrate-and-fire). In fact, incorporating non-linear effects of shunting inhibition into simple neuron models would make neural network simulations based on point models more realistic while keeping the formulation and simulation simple and efficient: there would be no need to simulate complex multi-compartment neuron models instead where effects of shunting inhibition are very well understood (Hao et al., 2009).
However, in the study (Zhou et al., 2013), the derivation was based on phenomenological arguments to describe the effects of shunting as observed in experiments. Here we suggest a different derivation of a single compartment model starting from a biophysical more realistic conductance-based three-compartment model. The advantage of our derivation is that the shunting strength κ can be approximated analytically in biophysical terms. Moreover, our formulation naturally includes the correct distance dependence of the shunting strength κ and the experimental results of Hao et al. (2009) can be captured well.
The usefulness of integrating non-linear dendritic processing into a simple neuron model is that network effects induced by shunting inhibition can now be analyzed in much simpler mathematical form. For that, we further extend our model to include multiple excitatory and inhibitory sites with specific layout on the dendritic branches and devised a network of mutually connected excitatory and inhibitory groups of neurons. The structure of our network model is similar to those analyzed earlier (Wilson and Cowan, 1972;Brunel, 2000), and derivations of its topology were suggested for many brain areas, such as for memory circuits (Hopfield, 1984), pre-frontal cortex computations (Machens et al., 2005), or decision making (Wang, 2002;Wong and Wang, 2006). In contrast to previous models, however, in our network inhibitory neurons act as global shunting gates, and thus introduce a new multiplicative non-linearity. In fact, this arrangement of inhibitory inputs-specifically targeting the peri-somatic regions to act as global shunting gates-has been found to be a common motif for certain basket cell types in the brain (Markram et al., 2004).
We show that our simplified mathematical form of non-linear dendritic inputs is capable of generating persistent activity with relatively low firing rate in simulations of networks of spiking neuron (around 15 Hz), as has previously been found for firing rate models using a more detailed incorporation of dendritic effects (Abbott, 1991a). Moreover, firing levels can be "adjusted" in a gradual manner by e.g., changing the excitatory drive.
Non-linear dendritic processing might be an important underlying mechanism for commonly observed spiking dynamics in the brain. Our work suggests that non-linear dendritic effects can be simply incorporated in networks of spiking neurons.
Parts of the results were previously presented on a conference (Zhang et al., 2011).

RESULTS
To get an intuitive idea of the non-linearity involved in integrating two synaptic inputs, let us first consider a conductance-based neuron model which receives a pair of excitatory and inhibitory inputs at the soma. The dynamics of the neuron can be written as Koch (2004) where v is the membrane potential of the neuron and C the membrane capacitance. g is the leaky conductance and E L the resting potential. g E and g I denote, respectively, the excitatory and inhibitory conductances, and E E and E I the corresponding reversal potentials. One finds from Equation (2) that the membrane potential in the steady state can be written as with the factor γ := g + g E + g I . Thus the excitatory conductance g E , as well as the inhibitory conductance g I is scaled by a factor (γ) involving both, the inhibitory and excitatory conductances. Therefore the integration of inhibitory and excitatory currents on the somatic potential is non-linear rather than an independent linear summation. Of course, shunting interaction in reality is more complicated than this simple case because its effect depends on the spatial configuration of excitatory and inhibitory inputs on the dendrite of a neuron. In the following, we will derive a simple single compartment neuron model which incorporates the effects of shunting inhibition via a multiplicative rule. After the derivation, we will show how different geometric arrangements can be modeled in a network and finally show how shunting inhibition can lead to persistent activity in a neural network.

On-path configuration
Let us consider a simple integrate-and-fire neuron model, which consists of a soma and two dendritic compartments. The neuron receives an excitatory and an inhibitory input at the locations E and I on the dendrite, respectively (see Figure 1). Let us assume that an inhibitory input is on the route from an excitatory input on its way to the soma, i.e., the position I is between E and the soma ( Figure 1A). This situation was called the "on-path configuration" (Koch et al., 1983). The subthreshold dynamics of this neuron can be written (compare to the equivalent circuit in Figure 1B): where v S , v I , and v E denote the local membrane potentials at the soma and the dendrite locations I and E, respectively. C S and C D are the membrane capacitances of the soma and the dendrite locations E and I, respectively. E L is the resting potential of the neuron. E I and E E are the reversal potentials of the inhibitory and excitatory currents, respectively. g S and g D are, respectively, the leaky conductances at the soma and the dendrite locations. g IS is the transfer conductance from the dendritic location I to the soma, and g SI the transfer conductance from the soma to the dendritic location I. g IE and g EI are defined accordingly. The excitatory and inhibitory synaptic conductances, g E and g I , respectively, are describing the opening of corresponding ion channels. If we assume simple exponential activation curves, the dynamics of the synaptic inputs can be written as Koch (2004) Thus, synaptic conductances are driven by presynaptic spike trains which are expressed as a sum of delta-functions, m δ(t − t m ), with t m denoting the moment of the m-th spike. τ E and τ I are the time constants for the excitatory and inhibitory synaptic conductances, respectively. w E and w I are the corresponding synaptic connection strengths.
The dynamic of the neuron model of Equations (4-6) is difficult to analyze because it involves the dynamics of three membrane voltages. Thus a further simplification is desirable especially when considering the dynamics of a large network of interacting neurons. In the following, we will further simplify the above model by a separation of time-scales approach.
From Equations (4-6), we note that the magnitudes of the time constants of the potentials at the soma and the dendrite locations I and E can be roughly estimated to be τ S ≈ C S /(g S + g IS ), τ DI ≈ C D /(g D + g SI + g EI ) and τ DE ≈ C D /(g D + g IE ). Since the capacitance increases linearly with the surface area of the membrane (see e.g., Koch, 2004), the membrane capacitance at the soma is much larger than that at the dendritic locations E and I, i.e., C S C D . Furthermore, we assume that the leak conductances g S and g D are of the same order and larger than the transfer conductances, i.e., g S g IS , g D g SI and g D g IE . Therefore, we have τ S τ DI and τ S τ DE . In fact, using the parameters from Table 1 it is τ DI ≈ 1.7 ms, τ DE ≈ 2.3 ms, and τ S ≈ 13.4 ms and therefore indeed τ S τ DE ≈ τ DI . This implies that the dynamics of v I and v E occur much faster than that of v S . Thus, we can effectively treat v S as a slow time variable, and v I , v E as fast time variables. The dynamics of the somatic potential can then be solved approximately by assuming that v I and v E reach their steady values instantly. This is achieved by setting the left-hand sides of Equations (5) and (6) to zero and solve for the membrane potentials. We obtain: Further, we assume that the input locations, E and I, and the soma are well separated on the dendritic branch (as this is the experimental condition during which (Hao et al., 2009) obtained their results). Substituting Equations (9) and (10) into Equation (4), Scaling factor α 5  we then get a simplified model for the dynamics of the somatic potential (see section 4.1 for a detailed derivation), where τ S = C S /(g S + g IS ), and To understand the behavior of the new neuron model (Equation 11), it is instructive to look at the steady state of the somatic potential in response to constant synaptic inputs. The steady statev S is obtained by setting the left-hand side of Equation (11) to be zero: Thus, when no inhibitory input is applied, i.e., g I = 0 and is the voltage change at the soma due to the excitatory input g E . Analogously, if no excitatory input is applied, i.e., g E = 0 and f d (g E ) = 0, then the voltage change at the soma is given by v S − E L = f p (g I ). Thus, in case of single inputs the multiplicative effects of the dendrite is reduced to an additive form as in the common formulation of an integrate-and-fire model. However, if both excitatory and inhibitory inputs are applied simultaneously, their joint effect on the somatic potential is given by that is a summation of the excitatory and inhibitory contribution when each of them was applied separately, and an additional product of their independent contributions, i.e., κ on f d (g E )f p (g I ).
The multiplicative term comes from the non-linear shunting process, and the coefficient κ on represents the shunting strength. The new neuron model (Equation 11) describes the subthreshold voltage dynamics at the soma if two synaptic sites, excitatory and inhibitory synapses, are arranged in on-path configuration. The summation of excitatory and inhibitory conductances agrees with the form found in a recent experiment [Hao et al. (2009); see also Equation (1)], that is a linear sum of individual excitatory or inhibitory effects, f d (g E ) + f p (g I ), and a multiplicative term involving both, excitatory and inhibitory currents f d (g E )f p (g I ), as well as a shunting strength factor κ on . The functions f d (g E ), for the distal excitatory synaptic site, and f p (g I ), for the proximal inhibitory synaptic site, correspond to the induced voltage change at the soma for a given conductance input at the respective synaptic sites. If input conductance g E is small, is approximately linear (see Figure 2A). For larger excitatory inputs, the function saturates to a positive value. Analogously, if the inhibitory conductance g I is small, g I g S , then the function f p (g I ) is approximately linear (Figure 2B), but similarly saturates to a negative value for larger inhibitory inputs.
Note that these functions relating the somatic effect of the synaptic conductances include the transfer conductances from one site to the other. They are thus dependent on the distance between the excitatory and inhibitory site, as well as the distance to the soma. If only passive cable properties are considered, transfer conductances in both directions, e.g., g EI and g IE , would be equal (Koch, 2004). In practice, however, depolarization of the membrane potential caused by excitatory currents is amplified by the existence of voltage-dependent ion-channels (Cook and Johnston, 1997;Lai and Jan, 2006). Thus transfer conductances of excitatory input locations to other parts of the dendrite  is increased in comparison to the transfer conductances from inhibitory sites. We thus set g EI = αg IE , with α > 1 to reflect this property. Moreover, since in on-path configuration excitatory currents will flow pass the inhibitory location, causing similarly an amplification through active channels, we further set g IS = αg SI . Note that the setting of α does not affect our qualitative results but will influence the size of the somatic voltage change in response to synaptic inputs (see below and Figure 2). In contrast to the approach of Zhou et al. (2013), in our neuron model (Equation 11) the shunting strength κ on is explicitly given in biophysical terms. From Equation (14) it can be seen that the shunting strength will be particular prominent if the resting potential and the inhibitory reversal potentials are similar E L ≈ E I , in agreement with experimental and early theoretical findings (Koch, 2004;Hao et al., 2009). Note further that κ on decreases with growing transfer conductance g IS . Thus if location I is set farther away from the soma (while still retaining the on-path configuration), the transfer conductance to the soma g IS will naturally decrease, and therefore the shunting will increase. Thus, κ on tends to have a larger value if inhibitory inputs are at a distal site of a dendrite in comparison to inputs at a proximal site, agreeing with experimental observations (Hao et al., 2009). Moreover, in contrast to earlier two-port analysis (Koch et al., 1983;Hao et al., 2009;Zhou et al., 2013), our approximation of κ on does not depend on the transfer conductance between E and I, therefore the shunting strength is approximately constant when the excitatory synapse location is increasingly distal but the inhibitory location is fixed. This again agrees well with experimental findings (Hao et al., 2009).
Taken together, we found that in on-path configuration a single somatic point model can be derived starting from a three-compartment model when assuming that dendritic processing is fast compared to the somatic integration and that locations E and I are well separated. In this case the arithmetic rule (Equation 1) as suggested by Hao et al. (2009) is well captured [compare to Equation (15)]. We tested the accurateness of the simplifications by comparing the dependence of κ on on the input conductances. Note that in the simplified model (as in the suggested rule of Hao et al., 2009) κ on is independent of g E and g I . If we instead attempt to compute the shunting strength directly from the full three-compartment model by assuming that the arithmetic rule (Equation 1) was correct, we find that κ on-full is indeed almost not dependent on the input conductances. In detail, we calculated the steady state somatic voltage of the full three-compartment model, and solve for κ according to the arithmetic rule (Equation 1). That is, we first subtract both individual excitatory and inhibitory contributions from the steady state voltage [by setting g I = 0 or g E = 0, respectively) and then divide the result by both individual contributions (compare to Equation (1)] to get the shunting strength (which might in this case still depend on the input conductances g E and g I ).
As plotted in Figure 3A (solid line), the computed κ on-full from the full model (using the parameters of Table 1) changed very little for a large range of g I (or g E , not shown). This indicates that the suggested arithmetic rule of Hao et al. (2009) (Equation 1) is applicable and that our single compartment model is a good approximation to the on-path configuration.
In summary, we derived a new formulation of the subthreshold dynamics of an integrate-and-fire model with integrated dendritic processing [Equation (11); note that our form of the single compartment model is different from the formulation of Zhou et al. (2013)]. Our model incorporates effects of shunting inhibition of two synaptic inputs with on-path configuration.

Out-of-path configuration
Analogous to the on-path configuration, we can derive a simplified model for the out-of-path configuration, that is, when the excitatory synapse lies on the route of the inhibitory current to the soma. Formally in (Equation 4-6), excitatory and inhibitory synapses exchange position, but the equations otherwise remain the same (that is all labels E and I have to be interchanged in Equation 4-6). The resulting equations are similar, except that now f d (g E ) and f p (g I ) are replaced by f d (g I ) and f p (g E ) (with appropriately adapted labeling): To arrive at these equations we used analogous assumptions as in the on-path case. In particular, we assumed that both dendritic locations are well separated. Note, that in out-of-path configuration, the shunting strength is negligible, because the absolute difference of the somatic reversal potential (E L ≈ −80 mV) and the excitatory synaptic reversal potential (E E ≈ 0 mV) is very large, thus κ out ≈ 0 mV −1 . This result agrees well with previous theoretical analysis (Koch et al., 1983), where it was found that the shunting strength in out-of-path configuration is negligible compared to that in on-path configuration.
In the derivation of the model for the out-of-path configuration we assumed that all synaptic locations are well separated and found that the arithmetic rule (Equation 1) is approximately satisfied. However, note that our approximation is less successful than in on-path configuration. In particular, if compared to the shunting strength derived from the full model in out-of-path configuration (analogous to the calculation of κ on-full as described above), we find that κ out-full is strongly dependent on the input conductances, specifically on g I . This dependence is plotted in Figure 3A (dashed line).
Taken together, since both, the arithmetic rule (Equation 1) and our single compartment model (Equation 16), predict a constant shunting strength κ out in respect to the inputs, they are both less valid in out-of-path configuration than in on-path configuration, where κ on-full can indeed be regarded as constant.

Shunting strength depends on the distance between the synaptic sites and the soma
Since we derived our single compartment model (Equation 11) starting from a three-compartment model, distance dependence of the shunting strength can be qualitatively assessed by assuming distance dependence transfer conductances, namely g EI , g EI , g IS , and g SI . Biophysically, resistance will grow linearly with increasing length of a passive cable (Koch, 2004), thus we assume that each transfer conductance between two points on the dendrite is reversely linearly dependent on their distance x.
For estimating the distance dependence of the shunting strength κ, we follow the experiments of Hao et al. (2009) (their Figure 3), and first fix the inhibitory location I at three positions, I 1 , I 2 , and I 3 , and set g SI to 2.2, 3, and 5 nS, respectively. The other transfer conductances are given by their distance dependence, g trf = g max /(μx + 1). Here, g max is the conductance when two points are co-localized together (x = 0), and is much larger than the leaky conductance, we set g max = 300 nS. The spatial scaling factor μ depends on the electrotonic properties of dendrite (Koch, 2004) and determines the spatial scale, we set arbitrarily μ = 3 as it does not change the qualitative picture. Other transfer conductances are set according to the spatial arrangement of the input sites I and E, e.g., for the out-of-path configuration it is 1/g SE + 1/g IE = 1/g SI . The rest of the parameters are set according to Table 1. Figure 3B shows the variation of the shunting strength of the full three-compartment model for the on-path configuration (Equations 4-6) and the out-of-path configuration as a function of the distance of the excitatory site from the soma. Note that when varying this distance, the model switches from out-ofpath to on-path configuration at the site of the inhibitory input (marked with dashed-dotted lines in Figure 3B). We found that the shunting strength increases sharply for out-of-path configurations the nearer the excitatory and inhibitory sites are, whereas the shunting strength remains approximately constant for the on-path configuration. This distance dependence of our threecompartment model reproduces the experimental findings as well as complex model simulations using dendrites having 200 compartments (Hao et al., 2009). Our single compartment model assumed that inhibitory and excitatory sites are well separated and thus approximates the full model in the asymptotic regime (solid lines in Figure 3B).

NEURON HAVING MULTIPLE DENDRITES
Up to now, we analyzed the shunting contributions to a neuron receiving only two inputs, one excitatory and one inhibitory synapse. However, in real situations as well as in neural network models, a neuron will typically receive hundreds to thousands of input synapses. To investigate the effect of multiple input synapses, we here analyze three possible configurations of synapses on dendritic branches exemplifying potential excitatory and inhibitory interaction patterns (see Figure 4).

Single synapses on individual dendritic branches
In the first configuration, individual excitatory or inhibitory inputs are scattered on different dendritic branches (as illustrates in Figure 4A). In this case, the synapses are physically separated on parallel branches and thus shunting interaction between excitatory and inhibitory currents can be ignored. As shown in the section 4.2.1, the simplified model for the somatic potential v S is given by where g T i denotes the corresponding synaptic input conductance of the ith synapse and T is a reminder of the type of the ith synapse, either inhibitory or excitatory. The time constant of the somatic voltage is given by τ S = C S /(g S + Ng TS ), where g TS denotes the transfer conductance from an input site to the soma. We assume that the transfer conductances from synaptic sites to the soma are equal for all N dendritic branches. The function relating the synaptic conductances to the somatic voltage change is given by where g T i is the synaptic conductance on the ith dendritic branch and E T i the corresponding reversal potential. Thus, according to Equation (19), if inhibitory and excitatory synapses are located on individual branches, their individual somatic contributions can simply be added.

On-path configuration on each dendritic branch
In the second configuration (Figure 4B), each dendritic branch has a pair of excitatory and inhibitory synapses in on-path configuration. In this case, the simplified model can be written as (see section 4.2.2) where τ S = C S /(g S + Ng IS ), and Note that in this configuration, individual shunting components of each branch can be simply added together.

Global shunting
In the third configuration, only excitatory synapses are distributed on dendritic branches and a single inhibitory synapse is located in the peri-somatic region (see Figure 4C). In this configuration, the single inhibitory input might shunt all incoming excitatory currents, and we therefore call it global shunting. For simplicity, we assume that the inhibitory synapse is located very close to the soma, so that the inhibitory compartment can be identified with the somatic compartment. This targeting of the perisomatic region can be commonly observed for certain types of basket-cells (Markram et al., 2004). In this case, the following single compartment model can be derived (see section 4.2.3 for details) were τ S ≈ C S /g S and Thus, the general multiplicative form of the shunting effect remains the same as in the on-path configuration of our initial model (Equation 11), but with the difference that now multiple excitatory inputs are first added together before being multiplied by the inhibitory component. Note that if more than one inhibitory input is considered (all similarly targeting the perisomatic region), all contributions can be simply added up. That is g I in Equation (25) can then be replaced by a sum of all inhibitory inputs, j g I j .

NETWORK DYNAMICS WITH GLOBAL SHUNTING INHIBITION
In the above we have developed simple models describing the sub-threshold dynamics of point-neurons while integrating nonlinear dendritic processing. These simplified single-compartment models are valuable to analyze the effects of shunting inhibition on the dynamics of large-scale networks.
In the following, we will show by using our simplified neuron model that shunting inhibition might lead to persistent activity in spiking neuron networks as has been suggested previously for a different type of dendritic non-linearity (Abbott, 1991b). Our network structure of mutually interconnected groups of inhibitory and excitatory neurons (as illustrated in Figure 5) follows a common excitation-inhibition-type network layout (Wilson and Cowan, 1972;Brunel, 2000).

Persistent activity by shunting inhibition
In detail, we built a network model of mutually connected N E excitatory and N I inhibitory neurons, with N E = 4N I . Each neuron is sparsely connected to others with a small probability p = 0.1. For simplicity, the connection strengths between neurons are constants, with w E and w I denoting the excitatory and inhibitory strengths, respectively. We assume that each neuron in the circuit has pN E branches, one for each excitatory input. All inhibitory inputs target the per-somatic region. In our model,

FIGURE 5 | A network model for generating persistent activity.
Neurons in excitatory (Exc.) and inhibitory (Inh.) groups are mutually connected in random and sparse fashion. Inhibitory effect are governed by a (peri-somatic) global shunting inhibition.
both, excitatory and inhibitory neurons, are governed by nonlinear dendritic processing and thus include shunting inhibition. Transfer conductances are assumed identical for each branch (for detailed parameters settings, see Figure 6).
According to Equation (25), the dynamics of the ith neuron in the network is Here, the superscript T can either be T = E or I and denoted the type of the ith neuron. J T i represents the total synaptic input received by the neuron, which consists of excitatory, inhibitory and shunting components. τ T S , T = E or I, are the somatic time constants of excitatory or inhibitory neurons, receptively, depending on the corresponding capacitance; we set C E S = 740 pF and C I S = 370 pF. To model background noise, we include a Gaussian noise term of zero mean and variance σ 2 . g ext refers to a small external input.
We further assume that a neuron fires a spike if its membrane potential reaches a threshold of E thres = −50 mV. After a spike the membrane potential instantaneously is reset to E reset = −70 mV. The dynamics of the excitatory and inhibitory synaptic conductance are given by Equations (7) and (8), respectively. We adjust the time constant of excitatory conductance to τ E = 100 ms to account for the prominent role of slow NMDA conductances in the generation of persistent activity (Wang, 2002;Wong and Wang, 2006). Inhibitory conductance time scale is set to τ I = 10 ms. When parameters are chosen properly, we found that the network can retain sustained activity at a low firing rate even after an external stimulus is removed. An example of the network activity during persistent activity is shown in Figure 6A. Note that after a brief input (in the range 50-250 ms) the population firing rate remains high and does not return to the spontaneous activity level (which is induced by the noise term in the region 0-50 ms). Remarkably, the average population firing rate of the network does not reach very high firing rate but instead continuous to fire with a relatively modest rate of around 15 Hz for our parameter setting. An increase in the excitatory synaptic weight w E causes the firing level to gradually increase (from 13 to 18 Hz in our example, see Figure 6B).

Theoretical analysis of the persistent activity induced by global shunting inhibition
An intuitive picture for the underlying mechanism of the persistent activity induced by non-linear dendritic processing is as follows. When the firing rates of excitatory neurons are small, the self-excitation of excitatory neurons dominates and the firing rates of all neurons increase; however, with the excitation increasing further, the effect of shunting inhibition starts to grow dramatically in a non-linear manner due to its multiplicative dependence on the excitatory current. In case of persistent activity these two opposite effects are approximately balanced, so that neurons fire persistently at relatively low values without the need of external inputs. To further elucidate these ideas, we carry out mean-field approximation by considering the average synaptic inputs to the excitatory and inhibitory populations, and approximate the network dynamics. The mean-field approach has previously been used to analyze the dynamical processes of neural networks with integrate-and-fire neurons (Amit and Brunel, 1997a,b;Brunel and Wang, 2001;Renart et al., 2003;Wong and Wang, 2006). Basically, as the network is composed of identical neurons, the synaptic input to each neuron in the network can be treated as a Gaussian random process. Thus we can use a single variable to represent population firing rates; let r E and r I denote the firing rates of excitatory and inhibitory neural populations, respectively. The population firing rates depend on synaptic currents, which in turn depends on firing rates. Thus, the population firing rate of a steady state can be determined by assuming self-consistency. In the analysis, we only consider the mean synaptic inputs to a neuron and neglect the input variance as it is of no importance in our setting (Wong and Wang, 2006).
From Equations (7) and (8) we find the dynamics of the excitatory and inhibitory synaptic conductances in the rate limit as If one then approximates the non-linear spiking mechanism by a threshold linear function (see e.g., Mongillo et al., 2008), If one now uses linear functions to approximate f Gd and f Gp in the range of their average inputs according to Equation (30) one can straightforwardly solve for the population firing rate (see section 4.3.1). In our simulations, the mean-field approximation matched the simulation quite well (see Figure 6B, dashed line). Finally, the stability of the fixed-point of the population firing rate can be analyzed (see section 4.3.2). It turns out that both eigenvalues are negative in our parameter settings, indicating stability. Moreover, the shunting strength κ G has an stabilizing effect: with increasing κ G , eigenvalues will be decreased further. In summary, we confirm that our form of global shunting inhibition can cause spiking network models to exhibit persistent firing activity in a low rate regime.

DISCUSSION
In the present study, we have presented a new derivation of a "dendritic-integrate-and-fire" single compartment neuron model. We show that the model captures well the non-linear integration of excitatory and inhibitory inputs at the soma with an arithmetic rule, where the shunting effect is expressed as a product between the contributions of excitatory and inhibitory inputs, as was suggested by recent experimental finding (Hao et al., 2009). In contrast to attempts relying on phenomenological modeling (Zhou et al., 2013), we start the derivation from a biophysical description of a three-compartment model, which allows for a biophysical interpretation of the equations. We find that the dependence of the shunting strength on the distance of the synaptic locations is well captured with the three-compartment model. In particular, in on-path configuration, the shunting strength decreases with the distance of the inhibitory site from the soma and stays constant with the distance of the excitatory site to the soma. These properties are found in experiments (Hao et al., 2009), but are not well captured by a simple two-port analysis, where the derived shunting strength is still dependent on a term related to the transfer conductance between excitatory site and the soma (Zhou et al., 2013). The reason for this result of a two-port analysis is that currents from both synaptic inputs are assumed to reach the soma. However, the real situation of the on-path configuration as captured by our three-compartment model is that currents from the distant excitatory compartment have to pass through the inhibitory compartment on its way. Thus there is no direct connection from excitatory site to soma as assumed by the simple two-port analysis as derived by Zhou et al. (2013).
In our derivation of the multiplicative rule of dendritic integration, we made two main assumptions: (1) the transfer conductances between the dendritic compartments in the antidromic direction are negligible small and (2) dendritic computations are fast in respect to the soma. The first assumption is only valid if the inhibitory and excitatory compartments are well separated on the dendritic tree. Thus the multiplicative rule refers only to the limit of well separated sites. Note that in the simulation (Figure 3B), the shunting strength is not constant for nearby input sites and the multiplicative rule is not applicable. Moreover, because of the depolarization for excitatory inputs, active channels increase the antidromic transfer to the inhibitory site in out-of-path configuration (which we modeled with the parameter α). Thus in comparison to on-path configuration, assumption (1) is less valid in out-of-path configuration (as illustrated in Figure 3A). Thus we conclude that the multiplicative rule can be mainly applied in the on-path configuration for well separated input sites. Second, we further assumed that time scales of dendritic and somatic processing are separated, which is a good approximation because the large area of the soma in comparison to dendritic compartments. Whether the leak conductances per area of the somatic and dendritic compartments are on the same order or somewhat higher in the dendrite [as suggested experimentally (Golding et al., 2005;Omori et al., 2006Omori et al., , 2009] is of no consequence to our analysis because the much larger membrane area of the soma will still ensure that the total leak is larger than in dendritic compartments. For our parameter setting, the time constant for the somatic membrane voltage is about 10 times slower than that of the dendritic compartments. Thus one can indeed assume that dendritic compartments instantly relax to steady state. Note that our assumption does not mean that the dynamics of synaptic conductances has to be fast. The dynamics of the synaptic input conductances were in fact not approximated and are still governed by their original form Equations (7) and (8). Thus our model can accommodate both, fast and slow ion-channels, such as AMPA and NMDA, respectively.
In a further analysis, we found that the arithmetic rule does much better apply to the on-path configuration than to the outof-path configuration. In fact, in the latter case the shunting strength itself is dependent on the size of the inhibitory input and thus the arithmetic rule is only partly valid. This is admitted, but not further analyzed, also in the experimental study (Hao et al., 2009) where the authors state that the rule is valid only up to a certain range of input conductances.
Since our general aim was to derive a simple neuron model based on biophysical properties to be used for investigating the effects of shunting inhibition on network dynamics, we included active properties only as a multiplicative factor for certain transfer conductances to arrive at simple models. To model additional effects of active channels on shunting inhibition the voltage dependent dynamics of conductances have to be incorporated has been attempted in a recent study (Jadi et al., 2012) suggesting that information processing capabilities of dendritic integration might be even richer than thought previously for passive dendrites (Vu and Krasne, 1992).

PERSISTENT ACTIVITY
Our main aim of this study was to derive a simple model for analysis of network effects of non-linear dendritic processing. As an example of this approach, we investigated an simple network consisting of mutual connected groups of inhibitory and excitatory neurons and show that global shunting inhibition can naturally induce persistent activity.
Neuronal persistent firing has been widely observed in neural systems and is believed to play important roles in cognitive functions. For instance, persistent activity might hold information of a memory trace thereby storing information about recent inputs (Wang, 2001;Machens et al., 2005). It has been observed that during persistent activity, neurons fire irregularly at low rates, typically less than 100 Hz (Gold and Shadlen, 2001). From a dynamical systems point of view, to maintain self-sustained activity in a network, it needs, on one hand, strong excitatory interactions between neurons to retain the excitation, and, on the other hand, inhibitory interactions to avoid the divergence of neuronal responses. In a firing-rate based network model, one often assumes a sigmoidal input-output (IO) function to avoid the explosion of the network activity (e.g., Hopfield, 1984). In a network model of spiking neurons, the saturation of synaptic current mediated by NMDA receptor is often assumed to account for bounded neuronal responses (Wang, 2002). Although these models are capable of sustained activity, it seems energetically unlikely that the nervous system maintains activity in a regime where the saturation of neural properties is essential. Thus it is probably that other non-linear mechanisms ensure that run-away excitation does not occur.
In an earlier study, it was shown that a non-linear form of dendritic processing is indeed enough to generate persistent firing (Abbott, 1991b). Similar to our findings, the author accounted for non-linear dendritic processing in a network model and analyzed the resulting dynamics with a mean-field approach, concluding that persistent activity can occur in the low firing regime. The difference to our analysis is that the author starts from an even more biophysical detailed view (from the full cable equation) resulting in more detailed equations. The derivation thus does not yield the intuitive picture of a multiplicative rule weighted with a shunting strength parameter κ, as suggested experimentally. We here focused on the biophysically meaning of this rule and thus provided a different derivation starting from a threecompartment model and investigated its validity for different input configurations. The advantage of this rule [in contrast to equations suggested by Abbott (1991b)] is that two different non-linear effects of dendritic processing are clearly separated: First, the non-linear somatic response functions (see Figure 2), which comprise non-linear effects related to the leaky signal transduction of dendrites. Second, the multiplicative term of excitatory and inhibitory currents describing the shunting inhibition. By explicit weighting of the multiplicative term by a shunting strength κ both non-linearities can be individually analyzed.
We show that the multiplicative non-linearity of the shunting inhibition is needed to generate sustained activity without saturation. Note that the persistent activity found in our network does not rely on the non-linear form of the functions relating the input conductances to the somatic voltage change Figure 2, because further analysis shows that this saturation occurs at a firing rate much larger than 100 Hz and hence is not relevant (not shown).
In our formulation we used a factor α to effectively incorporate active channels. For global shunting, α mainly regulates the slope of the inhibitory somatic response equation (Equation 26) through the transfer conductance g ES , while the total excitatory response (summed over all N dendritic branches) is approximately independent of g ES (see Equations 27 and 25). Although synaptic weights may need to be adjusted for a given persistent firing rate when α is varied, the qualitative picture of the network dynamics does not change.
Besides generating a sustained population activity of low firing rate, activity levels could also gradually adjusted by changing excitatory synaptic weights. Note that this graded change of activity level is not directly related to the "graded persistent activity" as found e.g., in the entorhinal cortex (Egorov et al., 2002). In the latter case, persistent activity level is changed by the amount of input given during a brief period. The underlying mechanism remains largely unknown (Brody et al., 2003). In contrast, in our case firing rate levels are adjusted by changing the recurrent synaptic weights.
This property of gradual adjustment of firing rates by varying the excitatory weight is induced by the strong multiplicative nonlinearity during shunting. Note that the inhibitory effect mediated through shunting inhibition is multiplicative in both the excitatory and inhibitory rates. This square dependence induces a strong and robust change in the activity level when the synaptic weight is varied. On the other hand, if persistent activity is induced by a saturation of the firing rate IO-function, a change of the synaptic weight will have very little effect because the IO-function will be flat in the saturation limit.
In summary, our results suggest that shunting inhibition might be an important underlying mechanism in the dynamics of neural networks.

THE SIMPLIFIED NEURON MODEL
We consider first there is only excitatory synaptic input. By setting g I = 0 in Equation (5) Substituting Equation (6) into (35), we obtain Re-organizing the above equation, we get Assuming g IE g EI (g D ) 2 , we get Substituting Equation (38) into (4), we obtain where f d (g E ) is given by Equation (12). Similarly, we can calculate the case when there is only an inhibitory synaptic input by setting g E = 0, and the result is where f p (g I ) is given by Equation (13). Finally, by combing the above results, we obtain the dynamics of the somatic potential when both excitatory and inhibitory synaptic inputs are applied (see Equation 11) where κ is given by Equation (14).
The key algebraic insight in the derivation of the multiplicative rule for the integration of the synaptic input conductances, Equation (41), is the following. We exemplify the key step in a simplified system. Note first that dendritic processing in the inhibitory compartment is linearly transmitted into the somatic compartment, given by Equation (4). We thus can focus on the steady state equation of the inhibitory compartment, as given by Equation (9) and set for now g SI = 0. If we assume that the transfer conductance back to the excitatory compartment is negligible [that is g IE g EI (g D ) 2 as above], the excitatory voltage is not dependent on the inhibitory voltage. Thus it is where we set h(g E ) ≡ g EI (v −E − E L ), γ ≡ g D + g SI + g EI , and β ≡ E I − E L , for the time being. Crucially, inhibition and excitation terms in Equation (43) are mixed due to the g I in the denominator. The key insight is that the term (Equation 43) can be expressed as a multiplicative rule separating the effects of g I and g E : which has the desired multiplicative form with f 1 (g I ) ≡ βg I / (γ + g I ) and f 2 (g E ) ≡ h(g E )/γ. The derivation of Equation (41) was done analogously albeit starting from the full model (including the somatic compartment).

Single synapses on individual dendritic branches
We first consider that single synaptic inputs are distributed on parallel dendritic branches, i.e., each branch receives only one type of synaptic input (see Figure 4A). In this case, the dynamics of the neural compartments can be written as Here, we use T = E or I to denote the type of an input synapse. v T i represent the potential at the location of the T-type input on the ith dendrite. g ST is the transfer conductance from the soma to an input site, which we assume equal for different dendritic branches. g TS is defined analogously.
Given that the dynamics of the voltage of the synaptic compartments, v T i , happens on much faster time scale than the somatic voltage v S , we can assume that they instantaneously relax to the steady state. By thus setting the left-hand side of Equation (51) Substituting Equation (52) into (50) and assuming that g TS g ST (g S ) 2 leads to the model as stated in the main text (Equation 19).

On-path configuration on each dendritic branch
We now consider the second configuration in which each dendrite branch comprises of a pair of excitatory and inhibitory inputs in on-path arrangement, that is the inhibitory synapse is located proximal and the excitatory synapse distal (in respect to the soma; see Figure 4B). The voltage dynamics of all compartments is then given by where v E i and v I i are the membrane potentials at the locations E and I, and g E i and g I i are the excitatory and inhibitory synaptic conductances of the ith dendritic branch, respectively.
Consider again that the voltages at the dendritic locations, v E i and v I i , can be regarded as fast variables, Equations (54) and (55) can be solved for the stationary state of v E i and v I i , respectively. Substituting the results into Equation (53), and again assuming that the product of the transfer conductances are much smaller than the leaks, g SI g IS (g S ) 2 and g IE g EI (g D ) 2 , we find the final form of the model as given in the main text (Equation 21).

Global shunting
In the third configuration, we consider that an inhibitory input targets the soma directly and excitatory inputs are distributed on parallel dendritic branches (see Figure 4C). The dynamics of the compartments can be written as Note that we do not need to explicitly model an inhibitory compartment here because the inhibitory input directly affects the somatic potential. Thus the somatic time constant, given by τ S = C S /(g S + g IS ) when explicitly modeling the inhibitory compartment, is now approximated by τ S ≈ C S /g S since the transfer conductance from inhibitory compartment to the soma g IS is neglected. Following analogous assumptions as above we get the Equation (25) in the main text.

Fixed point of population firing rates
We first approximate the functions f Gd and f Gp in the range of their average inputs, that isg E and N I pg I [compare to Equation (30)]. This yields f Gd (g E ) ≈ ag E + b and f Gp (N I pg I ) ≈ Frontiers in Computational Neuroscience www.frontiersin.org May 2013 | Volume 7 | Article 56 | 12 cN I pg I + d, where a > 0, b ≥ 0, c < 0, d ≤ 0 are constants (see Figure 6 for the numerical values). Now we find for the inputs: J T = N E p(ag E + b) + cN I pg I + d (58) + κ G N E p(ag E + b)(cN I pg I + d) Since the conductances and population rates have to match during a steady state (self-consistency), we can in the following calculate the fixed point of the population rates,r E andr I . Using Equation (58) together with the steady states of Equations (31) and (32), we find first for the inputs with the shortcuts A 3 = κ G aw E τ E pN E cw I τ I pN I (62) Combining Equation (59) with Equations (33) and (34) finally yields the fixed point of the population firing rates where B 1 = μ E A 4 − μ E β, B 2 = μ E A 1 + μ I A 2 − 1, and B 3 = μ I A 3 .

Stability
Let us analyze the stability of the stationary states, Equations 64 and (65). We thus have to consider the following dynamical equations τ I dg I dt = −g I + w I τ I r I =: G 2 (g E ,g I ) The fixed points of the above equations arē The stability of the population activity is determined by the eigenvalues of the Jacobian matrix J at the fixed points. The Jacobian matrix is given by whose eigenvalues are calculated to be Thus, when λ 2 < 0 the solution is stable. Note that since only c < 0 and d ≤ 0 the increasing shunting strength κ G will decrease λ 2 and thus tend to stabilizing the system.