A minimal model for synaptic integration in simple neurons

Information processing in neurons is a complex phenomenon that is affected by, among other things, the kind of stimulus encoded, type of nerve cells involved, ion channels and dendritic morphology. Individual neurons are typically specialized enough to be treated as computational units, capable of performing a variety of tasks [1]. The most prominent of these is the integration of information from presynaptic inputs. Each synapse, when stimulated by the appropriate neurotransmitter, gives rise to excitatory postsynaptic potentials (EPSP) which propagate down the dendrites and are “summed” in the cell body of the postsynaptic neuron [2, 3, 4]. Subsequently, an action potential is generated when the membrane potential crosses a threshold. For many important classes of neurons, additional processes modulate synaptic integration. These include nonlinearities in dendritic responses to postsynaptic currents [5], variability in dendritic structure, inputs from inhibitory synapses and the active properties of dendrites [6]. The rules of summation usually involve nonlinearities themselves in that the net activation in the cell body of the neuron may be slightly more or less than the arithmetic sum of the processed dendritic inputs [7, 8]. Understanding how these factors affect synaptic integration requires detailed biophysical simulations incorporating several spatial and temporal scales. Nevertheless, synaptic integration in its bare minimum role as a summation of dendritic impulses can still demonstrate the functional advantages of postsynaptic neurons. In this article, we pursue such an approach and describe a minimal model of synaptic integration. The model is minimal in the sense that synaptic currents apply simple perturbations to the phase space of the postsynaptic neuron which simplifies the analysis of the resulting dynamics. Moreover, a small number of biophysically relevant parameters are used. We make use of the well known FitzHugh-Nagumo (FN) equations to simulate single neurons. The FN model is prototypical of many excitable systems and is able to capture different modes of neuronal spiking behaviour. In his seminal work, FitzHugh [9] showed that the phase space of the FN model could be partitioned into distinct regions, each of which corresponded to a certain physiological “state” of the neuron such as active, refractory, etc. He also showed that the FN equations could be reduced from the famous Hodgkin-Huxley model [10] which is a realistic description of neuron electrophysiology. Thus, the FN model trades detail for simplicity, yet retains the essential features of neuronal excitability.
As a first approximation to modelling excitatory postsynaptic currents, we make use of rectangular pulses. These are obtained from presynaptic membrane potentials by simple application of a Heaviside function. The Heaviside function is suitably parameterized to control the height and width of the pulse current. Each synapse thus contributes a series of rectangular pulses to the total current stimulating the postsynaptic neuron. In the next section, we briefly introduce the dynamics of the FN equations and investigate the excitation threshold for different numbers of synapses. Next, we formulate the model and explore its validity and limitations. Then, we illustrate a possible application by exploring two noise driven processes: stochastic and coherence resonance. Finally, we highlight the scope for extending the model to incorporate other aspects of synaptic integration without compromising the simplicity of implementation. Biologically relevant observations are highlighted throughout. In starting with a minimum set of parameters and a simple model for synaptic inputs, we take a bottom-up approach that allows us to isolate the functional properties of synaptic integration made possible by summation of excitatory inputs alone (as opposed to other mechanisms like dendritic processing for instance).

The FN equations consist of two state variables whose dynamics operates on significantly different timescales. The fast “voltage” variable, $v$, mimics a neuron’s membrane potential and is complemented by a slow “recovery” variable, $w$. The dynamical equations are given by:

where the parameters $a$ and $b$ have values $0.5$ and $0.15$ respectively and $\epsilon$ represents the ratio of the two timescales and is typically set to $0.005$. A time dependent current injection into the neuron is modelled by the function $I$. In contrast to the all-or-none nature of neuronal firing patterns, the FN model does not exhibit a well defined threshold for excitability [11, 12]. That is to say, action potentials with a range of amplitudes are possible depending on the choice of the current term $I(t)$. However, there is a sense in which one can approximate a firing threshold in the FN model and this requires the analysis of a family of trajectories known as canards [13, 14, 15, 16, 17]. We assume that the initial condition lies on the intersection of the nullclines and that the current term $I(t)$ is null at $t=0$ (Fig. 1(a)). Physiologically, this corresponds to a neuron at its resting membrane potential. If at some arbitrary time, a current is applied in the form of a rectangular step stimulus (of infinite duration), the cubic nullcline $(\dot{v}=0)$ shifts upwards and the formerly stable initial condition is now below the new fixed point. Depending on how much the nullcline shifts upward, the resulting trajectory can be either a small-amplitude excitation (Fig. 1(b)), an intermediate amplitude excitation (Fig. 1(c)), or a large amplitude excitation that resembles an action potential (Fig. 1(d)).

The height of the instantaneous step current therefore determines the amplitude of the resulting excitation. Moreover, the transition from small amplitude excitations to large amplitude excitations occurs in an exponentially small ($O(e^{-\epsilon})$) region of the parameter range $I$. Within this range, trajectories known as canards track the middle branch of the cubic nullcline for a relatively long ($O(\frac{1}{\epsilon})$) time (Fig. 1(c)). The middle branch of the cubic nullcline is called the repelling slow manifold since nearby trajectories diverge sharply from it. The right and left branches in contrast are the attracting slow manifolds since nearby trajectories converge onto them. Informally, canards are trajectories that follow both the attracting and repelling branches of the slow manifold for a considerable amount of time [15]. Fig. 1(c) shows a maximal canard, i.e. a trajectory that tracks the unstable branch of the cubic nullcline for the longest time. This trajectory passes through the local maximum of the cubic nullcline. The maximal canard is also called the quasithreshold separatrix since neighbouring trajectories diverge sharply to its left or right [9]. The quasithreshold separatrix can loosely be thought of as the neuron’s threshold if we treat all trajectories to its left as subthreshold responses and trajectories to its right as action potentials.

If the parameter $I$ is made sufficiently large so that the fixed point lies on the middle branch of the cubic nullcline, a stable limit cycle results (via a supercritical Hopf bifurcation). Here too, for an exponentially small range of $I$ values beyond the Hopf bifurcation point, small amplitude limit cycles transition into large amplitude limit cycles in much the same way that small amplitude trajectories transitioned into large amplitude ones in Fig. 1. For the rest of this paper, we assume that the fixed point lies on the left branch of the cubic nullcline. In other words, we will only be concerned with the FN model in its excitable regime (as opposed to the limit cycle regime beyond the Hopf bifurcation point).
Suppose that an instantaneous step current of infinite duration (as applied in Fig. 1) is applied to the system. We can obtain a stimulus-response curve (Fig. 2) by plotting the voltage difference from the resting value ($v_{\rm{max}}-v_{0}$) that a given step current $I$ is able to generate. We notice a sharp transition just beyond $I=0.02$ from small to large amplitude excitations. This is sometimes referred to as a canard transition/explosion. The sharpness of the transition allows us to approximate neuronal firings in the FN model as all or none, since intermediate amplitude excitations are sandwiched into a small threshold manifold that is rarely ever encountered by trajectories.

Depending on the neurotransmitter and number density of its receptor, an action potential in the presynaptic neuron will lead to a small current injected into the postsynaptic dendritic membrane. This in turn generates a small amplitude excitatory postsynaptic potential (EPSP). The amplitude of a single EPSP in most cases is insufficient to elicit postsynaptic firing. This is largely due to attenuation in passive dendrites. Typical EPSP values from individual synapses range anywhere from $\frac{1}{15}$ to $\frac{1}{400}$ of the firing threshold [18]. Figure 2 shows that there is a monotonic relationship between current amplitude and the EPSP height for subthreshold responses (before the canard transition). Thus, if a sequence of rectangular pulses of identical heights overlap in step, at some point the current will exceed the threshold manifold. We scale the height of each current pulse such that a minimum number, $N_{\rm{min}}$, of current pulses can generate an action potential. A threshold value for firing is chosen (to the right of the canard transition in Fig. 2), denoted $I_{\rm{T}}$. The height of each current pulse is therefore $\frac{I_{\rm{T}}}{N{\rm{min}}}$. We denote the width of the pulse by $\tau$. A higher value for $\tau$ indicates a longer effective time available for a sequence of pulses to overlap and cross the firing threshold. As in real neurons, the current pulses must overlap quickly enough to prevent the subthreshold membrane potential from falling. We can visualize this by letting $N_{\rm{min}}$ pulses each of height $\frac{I_{\rm{T}}}{N{\rm{min}}}$, overlap with an identical time gap denoted by $\mu$. Such a “staircase current” is an idealization but the width of the staircase ($\mu$) can be thought of as the average time between successive overlaps in the real case. Figure 3 shows the corresponding threshold curves for different settings of $N_{\rm{min}}$. We notice that there is effectively a maximum time gap exceeding which, the pulse current sequence fails to evoke a postsynaptic action potential. The maximum time gap is seen to decrease as the number of synaptic inputs ($N_{\rm{min}}$) increases. This is because a large number of small height pulses apply tiny upward shifts or nudges to the cubic nullcline. Consequently, the trajectory responds by drifting towards the new fixed point (since the net current is still subthreshold). For higher $N_{\rm{min}}$, the nullcline undergoes many small magnitude, upward nudges. The smaller the nudge, the better the trajectory is able to track the shifted nullcline. This means that the component of the trajectory’s drift in the $w$-direction is higher for high $N_{\rm{min}}$. As a result, by the time the final pulse arrives, the voltage of a neuron with larger $N_{\rm{min}}$ would have climbed less than one with a lower $N_{\rm{min}}$ (due to increased drift in the $w$-direction). This leads to the canard transition for higher $N_{\rm{min}}$ occurring earlier than that for lower $N_{\rm{min}}$.

Low time gap requirements are especially relevant to neurons that operate as coincidence detectors [1, 2]. The temporal resolution of coincidence detection is largely determined by parameters intrinsic to the neuron such as its membrane time constant. Figure 3 shows that even with these parameters remaining constant (no change in $a$, $b$ and $\epsilon$), modulating synaptic currents towards lower values can improve this temporal resolution.

Neurons are also stimulated by inhibitory, hyperpolarizing currents called anodal currents. In some cases, it is possible for an anodal current to cause an action potential. The mechanism by which this happens is referred to as anodal break excitation [19]. In the FN phase space, anodal break excitation is possible because the cubic nullcline shifts downwards when a hyperpolarizing current is applied (negative values for $I(t)$). The trajectory responds by settling onto the new fixed point. When the anodal pulse terminates, the nullcline shifts upwards to its former position which can leave the state point sufficiently below the original fixed point to evoke firing. Although our focus here is on the analysis of excitatory stimuli, threshold curves analogous to Fig. 3 can be obtained for anodal currents. This is shown in Fig. 4.

The observations are qualitatively similar to the excitatory case. There are certain features of the threshold curves in Fig. 3 that we have not explained yet. For instance, we notice that in the suprathreshold regime (low values of $\mu$), the threshold curve goes through an inverted S-shape/inflection point (just before the canard transition). This can be explained as follows. The maximum voltage ($v_{\rm{max}}$) attained by any trajectory corresponds to the point where the trajectory crosses the cubic nullcline from the left (eg. Fig 1d), since to the right of this point, the flow on the v-axis reverses direction. For low enough values of $\mu$, the trajectory hits the nullcline corresponding to the sum of all pulses. However, for values of $\mu$ near the canard transition, the first pulse ends before the trajectory manages to hit the final nullcline. Now, $v_{\rm{max}}$ is determined by values that range from the trajectory’s intersection with the penultimate nullcline, to the trajectory’s intersection with the final nullcline. This is what causes the small inflection point. In the subthreshold regime (beyond the canard transition), we notice points where the monotonic decrease of $v_{\rm{max}}$ is non-smooth. This can be explained similarly: $v_{\rm{max}}$ is determined by the penultimate and lower nullclines. In other words, the trajectory reverses direction after crossing the nullcline of the $n\textsuperscript{th}$ pulse and traverses far enough that all pulses $>n$ fail to get it past its previously attained maximum voltage.

Synaptic integration of inputs exhibits two modes known as spatial and temporal summation [2, 3, 4]. Spatial summation involves the simultaneous overlap of two or more dendritic inputs and a subsequent firing. In temporal summation, high frequency stimuli from the same dendritic input can summate to evoke firing. This requires that the time lag between stimuli be much smaller than the characteristic decay time of the postsynaptic potential. We remark here that the model we use can only exhibit spatial summation. This stems directly from the dynamics of the FN system, which belongs to a class of neuronal models termed resonators [15, 20]. Real neurons can temporally integrate subthreshold excitatory currents if these currents arrive with sufficiently low time gaps. In our model though, even if two or more subthreshold pulses follow each other with zero time lag, then that simply corresponds to a single subthreshold pulse of finite duration which anyway cannot cause firing.

We shall treat voltages obtained from FN trajectories as action potentials whenever a threshold value, $v_{\rm{T}}$, is crossed from below. We choose $v_{\rm{T}}=0.9$ which is close to the action potential peak. Because we include noise in the presynaptic FN equations, it is possible that during the course of an action potential $v$ may cross $0.9$ several times by virtue of stochastic fluctuations. To avoid such overcounting, we specify that a low voltage ($v=0.2$) be crossed from below in between two crossings of $v_{\rm{T}}$.

Equation (1) is modified to include an additive noise term as follows:

where $\xi(t)$ is white noise corresponding to a normalized Wiener process with autocorrelation $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$, and $\sigma$ (which we shall call the noise intensity) is the standard deviation of the Wiener increments. Equation (2) remains unchanged. How the presynaptic neurons themselves are stimulated depends on the type of phenomenon studied. For instance, we can let the presynaptic neurons be driven by a sinusoidal current stimulus as follows:

Denote the firing times of the $i^{\rm{th}}$ neuron by the set $\{t_{\rm{j}}\}$. The sequence of current pulses that the $i^{\rm{th}}$ synapse generates is therefore

where $\tau$ is the pulse width and $H$ denotes the Heaviside function. The net stimulus from all synapses is then

where $I_{T}$ is the firing threshold described in section 2, $N$ is the total number of synapses and $N_{\rm{min}}\in\mathbb{N}$ is the minimum number of presynaptic firings requires to elicit a postsynaptic action potential. In other words, if $N_{\rm{min}}$ pulses overlap in close succession, then the total postsynaptic current equals $I_{\rm{T}}$ which evokes firing. Postsynaptic neurons are kept noiseless for the simulations in this paper (i.e., $\sigma=0$). We choose this because as stated in the introduction, we are concerned with synaptic integration in its minimum role as a summation of excitatory stimuli. This will however require us to state assumptions on the source of noise that justify our choice of keeping postsynaptic neurons noiseless. For instance, we may choose to study synaptic integration of inputs from sensory neurons whose stimuli are affected by environmental noise. This would then justify keeping the corresponding postsynaptic neuron noiseless as it has access only to its presynaptic neurons and not the external environment.

In a sequence of input pulses overlapping to generate an action potential, there is always a time lag between one pulse and the next. These time lags cause a net drift of the trajectory in the leftward region of the quasithreshold separatrix of the nullcline set by the final pulse. The final pulse pushes it to the right of the separatrix and causes firing. While setting $I_{\rm{T}}$, we must avoid values close to the canard transition in Fig. 2. It is possible that for such values of $I_{\rm{T}}$, even if $N_{\rm{min}}$ pulses overlap in very close succession, the net drift can still bring the trajectory to the left of the quasithreshold separatrix of the final nullcline and prevent firing.

The stochastic differential equations corresponding to Equations(3,2) respectively are:

The numerical integration of Equation(7) requires a choice of interpretation of the stochastic integral used to integrate the noise term $\frac{\sigma}{\epsilon}dW$ (called the Itô and Stratonovich interpretations). For additive white noise however, both interpretations evaluate to the same stochastic integral so this does not affect the choice of which algorithm we use. We choose the simple first order Euler-Maruyama scheme [21, 22] which gives:

where $\Delta W_{n}$ is sampled from a Gaussian distribution with unit variance and $\delta t=0.001$. Because the Euler-Maruyama method is stiff, interspike interval histograms of neuronal firings were compared with those obtained using a second order stochastic Runge-Kutta method [23] for typical noise intensities used in this paper. The same was repeated for step size an order of magnitude lower. No significant differences were found. Interspike interval histograms (normalized) are a good approximation to the probability distribution functions of waiting times between firings. For obtaining postsynaptic trajectories, we use a standard deterministic Runge-Kutta method ($4^{\rm{th}}$ order) with step size $0.001$ (since they are kept noiseless). Sample trajectories are shown in Fig. 5 for a postsynaptic neuron stimulated by $5$ synapses.

There are two caveats with defining the postsynaptic input function as in Equation(6). Suppose we choose a relatively high value for $I_{\rm{T}}$, say $0.08$ and that $N_{\rm{min}}=2$. Now if a single pulse alone stimulates the postsynaptic neuron, then $I_{\rm{post}}(t)$ will be $\frac{I_{\rm{T}}}{2}=0.04$ which is suprathreshold (see Fig. 2). Thus, a single pulse will cause firing when it is not supposed to. The second problem arises when stimulating the neuron with high number of inputs, $N$. Suppose $N=100$, $N_{\rm{min}}=5$ and that $I_{\rm{T}}=0.03$. If $50$ inputs undergo spatial summation, then $I_{\rm{post}}(t)=\frac{50*0.03}{5}=0.3$ which is well beyond the Hopf bifurcation point for $I$ (leading to limit cycles and further on, into an unstable node). Thus, large values of $N$ and small values of $N_{\rm{min}}$ can lead to bifurcations in the dynamics. This is not necessarily undesirable from a modelling perspective. For example, a Hopf bifurcation caused by large $N$ and small $N_{\rm{min}}$ might help model interesting features such as bursting dynamics [18]. We will however not deal with such a scenario here and assume that the neuron is always in its excitable regime (where the fixed point lies on the left branch of the cubic nullcline). It is possible to artificially correct for high values of $I_{\rm{post}}(t)$ by composing it with another Heaviside function. However, we avoid this and simply point towards the validity of the model by outlining the constraints on the parameters $I_{\rm{T}}$, $N$ and $N_{\rm{min}}$: First, we would like to prevent any number of presynaptic inputs less than $N_{\rm{min}}$ from exciting the postsynaptic neuron. The worst case corresponding to this is when $N_{\rm{min}}-1$ inputs fire (nearly) simultaneously. Then, $I_{\rm{post}}(t)=\frac{(N_{\rm{min}}-1)I_{\rm{T}}}{N_{\rm{min}}}$ which we would like to limit to some subthreshold value $I_{\rm{sub}}$(chosen appropriately from Fig. 2). Thus, we have

If bifurcations are to be strictly avoided, then we limit the maximum value of $I_{\rm{post}}(t)$ to the value $I_{\rm{Hopf}}$ where the Hopf-bifurcation occurs ($I_{\rm{Hopf}}\approx 0.11$):

This constraint can be relaxed if such high values for postsynaptic currents are rare, or if limit cycle spikes are desired. We have said that $N_{\rm{min}}$ is a positive integer. There is no particular reason why $N_{\rm{min}}$ should be an integer except that doing so would provide a natural interpretation to the current term for postsynaptic neurons. The value of $N_{\rm{min}}$ along with $I_{\rm{T}}$ and the summation term in Eq. (6) determines how much the cubic nullcline shifts upward and so a non-integer value for $N_{\rm{min}}$ should also work. However, we will use integer values for $N_{\rm{min}}$ throughout. In the biological context, constraints (11-12) might correspond to simple neurons such as cerebellar granule cells which on average receive only 4 excitatory inputs from short dendrites [24].

In the following section, we explore the effects of excitatory stimuli only. However, by the simple addition of a term to equation (5), we can incorporate the effects of inhibitory synapses. For a neuron with $N_{\rm{1}}$ excitatory inputs and $N_{\rm{2}}$ inhibitory inputs, the postsynaptic current will be as follows:

where $S_{\rm{i}}^{+}$ and $S_{\rm{i}}^{-}$ are the sequence of Heaviside pulses from excitatory and inhibitory synapses respectively, $n_{\rm{1}}$ and $n_{\rm{2}}$ are the minimum number of presynaptic firings required to elicit a postsynaptic action potential, when excitatory or inhibitory inputs are considered alone. The parameters $I_{\rm{1}}$, $I_{\rm{2}}$ are the corresponding firing thresholds. How to set the inhibitory firing threshold, $I_{\rm{2}}$, depends on the kind of phenomena we would like to investigate. For example, if we want to explore firings generated by inhibitory synapses via anodal break excitation, then $I_{\rm{2}}$ can be chosen appropriately from Fig. 4. On the other hand, if we are concerned only with the role of inhibiting/ hyperpolarizing currents, then $I_{\rm{2}}/n_{\rm{2}}$ can be replaced with a single parameter that describes the height of each inhibiting pulse.

To illustrate some use cases of the model, we briefly investigate two well known noise driven phenomena: stochastic and coherence resonance. We compare presynaptic and postsynaptic responses in these effects and identify postsynaptic features that may be biologically significant.

We have investigated a simple model for synaptic integration using rectangular pulses. The model is minimal since the phase space of the FN system is perturbed by a finite sequence of instantaneous changes (via rectangular pulses) as opposed to continuous changes in parameter values. We have illustrated the use of the model in probing noise assisted/noise driven phenomena where the noise source is assumed to be external to the system and relatively larger than intrinsic sources of noise such as channel noise and synaptic noise [26]. A natural extension would be to let the postsynaptic responses be affected by noise as well. This is particularly relevant to internal neurological systems such as cortical networks. Channel noise for the postsynaptic case can be simulated with nonzero $\sigma$ in Eq. 3. Synaptic noise can be implemented by using an additional noise term to perturb the pulse width, $\tau$. This lets the effective time gap for spatial summation fluctuate, mimicking the stochastic opening and closing of ion channels.

The neurocomputational properties of a model depend on the types of bifurcations it undergoes [15, 20]. The FN model is capable of exhibiting interesting features such as tonic spiking, phasic spiking, class-1 excitability, rebound spikes, accommodation, etc [41]. A further point of inquiry would be to see how postsynaptic responses vary for each of these properties as a function of parameters $N$, $N_{\rm{min}}$, $\tau$. For neurons stimulated by suprathreshold signals, different $m:n$ phase locking patterns are possible (where $m$ spikes occur every $n$ cycles) [42]. It would be interesting to see how postsynaptic neurons respond to these deterministic phase locking patterns of presynaptic neurons.

For the stochastic resonance simulations, we have used signals that are in phase with each other. This is under the assumption that the presynaptic (sensory) neurons are spatially close so that they are effectively stimulated by identical signals. We included phase differences so that $I_{\rm{pre}}(t)=A\sin(2\pi ft+\phi)$ where $\phi$ is randomly sampled and independently set for the presynaptic neurons. If $\langle\phi\rangle$ is not small, the ISIH peaks of the postsynaptic response was found to decrease. Thus, anatomical features become relevant for the ability to exploit background noise if the neurons are separated enough.

We have used uncorrelated, white noise for our simulations. It has been shown both with experiment and analysis that coloured noise such as $1/f$ and $1/f^{2}$ noise can exhibit stochastic resonance and moreover enhance it [43, 30]. Implementing this to the postsynaptic case is an important extension of the results shown here.

It is known that single spikes produce unreliable synaptic transmission, with the probability of transmission ranging from $<0.1$ to $>0.9$ [18]. This suggests the functional importance of bursting patterns which can enhance transmission. A neural code that involves presynaptic neurons firing coincident bursts is thus likely to be robust. Bursts can be simulated in the model by adding a constant term to $I(t)$. This places the fixed point of the FN system closer to the Hopf bifurcation point. For stimulus frequencies that are smaller than the frequency of the limit cycle, we can have several spikes per stimulus cycle, mimicking the bursting dynamics. Note that $I_{\rm{T}}$ will have to be obtained from the threshold plot computed by including the constant term added to $I(t)$. Probabilistic synaptic transmission is incorporated by introducing a random variable $X\in(0,1]$: Spikes are generated when the neurons cross the firing threshold $I_{\rm{T}}$ and when $X(t)<p$ (uniform sampling) where $p$ is the probability of transmission.