## Abstract

We study the effect that the heterogeneity present among the elements of an ensemble of coupled excitable neurons has on the collective response of the system to an external signal. We consider two different interaction scenarios, one in which the neurons are diffusively coupled and another in which the neurons interact via pulse-like signals. We find that the type of interaction between the neurons has a crucial role in determining the response of the system to the external modulation. We develop a mean-field theory based on an order parameter expansion that quantitatively reproduces the numerical results in the case of diffusive coupling.

## 1. Introduction

Synchronized behaviour arising among the constituents of an ensemble is common in nature. Examples include the synchronized flashing of fireflies, the blossoming of flowers, cardiac cells giving rise to the pacemaker role of the sinoatrial node of the heart and the electrical pulses of neurons. This global behaviour can originate from a common response to an external stimulus or might appear in autonomous, non-forced, systems. The theoretical basis for the understanding of synchronization in non-forced systems was put forward by Winfree (1967), who showed that the interaction—i.e. coupling—between the constituents is an essential ingredient for the existence of a synchronized output. The seminal work on coupled oscillators by Kuramoto (1975) offered a model case whose solution confirms the basic hypothesis of Winfree: while interaction helps towards the achievement of a common behaviour, a perfect order can be achieved only in the absence of diversity—heterogeneity—among the components of an ensemble. In the Kuramoto model, diversity manifests when the oscillators have different natural frequencies, those they display when uncoupled from each other. While this result holds for systems that can be described by coupled oscillators, recent results indicate that in other cases diversity among the constituents might actually have a positive role in the setting of a resonant behaviour with an external signal. This was first demonstrated in Tessone *et al.* (2006) and has since been extended to many other systems (Acebron *et al.* 2007; Perc *et al.* 2008; Chen *et al.* 2009; Gosak 2009; Postnova *et al.* 2009; Tessone & Toral 2009; Ullner *et al.* 2009; Wu *et al.* 2009; Zanette 2009). In the case of non-forced excitable systems, a unifying treatment of the role of noise and diversity has been developed in the study of Tessone *et al.* (2007).

Many biological systems, including neurons, display excitability as a response to an external perturbation. Excitability is characterized by the existence of a threshold, a largely independent response to a suprathreshold input and a refractory time. It is well established that the dynamical features of neurons can be described by excitable models: when a neuron is perturbed by a single impulse, the neuron can generate a single spike, and when perturbed by a continuum signal, a train of spikes with a characteristic firing frequency can appear instead. Although the creation and propagation of electrical signals have been thoroughly studied by physiologists for at least a century, the most important landmark in this field is due to Hodgkin & Huxley (1952), who developed the first quantitative model to describe the evolution of the membrane potential in the squid giant axon. Because of the central importance of cellular electrical activity in biological systems and because this model forms the basis for the study of excitable systems, it remains to this date an important model for analysis. Subsequently, a simplified version of this model, known as the FitzHugh–Nagumo (FHN) model (FitzHugh 1961) that captures many of the qualitative features of the Hodgkin–Huxley model, was developed. The FHN model has two variables, one fast and one slow. The fast variable represents the evolution of the membrane potential and it is known as the excitation variable. The slow variable accounts for the K^{+} ionic current and is known as the recovery variable. One virtue of this model is that it can be studied using phase-plane methods, because it is a two-variable system. For instance, the fast variable has a cubic nullcline, while the slow variable has a linear nullcline, and the study of both nullclines and their intersections allows one to determine the different dynamical regimes of the model. The FHN model can then be extended by coupling the individual elements, and provides an interesting model, for example, of coupled neuron populations.

In the initial studies of the dynamics of the FHN and similar models, the individual units (e.g. neurons) were treated as identical. However, it is evident that real populations of neurons display a large degree of variability, both in morphology and dynamical activity; that is, there is a diversity in the population of these biological units. Then, it is natural to ask what role diversity plays in the global dynamical behaviour of these systems and a lot of activity has been developed along these lines in recent years. In general, it has been found, as noted above, that diversity can in fact be an important parameter in controlling the dynamics. In particular, it has been shown that both excitable and bistable systems can improve their response to an external stimulus if there is an adequate degree of diversity in the constituent units (Tessone *et al.* 2006).

In this paper, we continue our study of the effect of diversity, investigating in some depth its role in the FHN model of excitable systems. We consider in detail a system of many neurons coupled either chemically or electrically. We show that in both cases a right amount of diversity can indeed enhance the response to a periodic external stimulus and we discuss in detail the difference between the two types of coupling. The outline of the paper is as follows. In §2, we present the FHN model and the coupling schemes considered. In §3, we present the main results we have obtained from numerical simulations. In §4, we develop an approximate theoretical treatment based on an order parameter expansion, which allows us to obtain a quantitative description of the behaviour of the model and compare its predictions with the numerical results. Finally, in §5, we summarize the conclusions.

## 2. FitzHugh–Nagumo model

### (a) Dynamical equations

Let us consider a system of excitable neurons described by the FHN model. The dynamical equations describing the activity of a single neuron are
2.1
and
2.2
where *x*(*t*) and *y*(*t*) represent, respectively, the fast membrane potential and the slow potassium gating variable of a neuron. We assume that the time scales of these variables are well separated by the small parameter, *ϵ*=0.01. Other parameters are fixed to *b*=0.5, *c*=4.6 and *d*=0.1 (Glatt *et al.* 2008), while the value of *a* will fluctuate from one neuron to another, thus reflecting the intrinsic diversity in the neuronal ensemble.

Let us first concentrate on the dynamics of a single neuron as described by equations (2.1) and (2.2). This dynamics has a strong dependence on the parameter *a*. Three different operating regimes are identifiable: for , the system has a stable focus in the right branch of the cubic nullcline leading the system to an excitable regime; for a limit cycle around an unstable focus appears (oscillatory regime); and for a stable focus appears again, now at the left side of the cubic nullcline (excitable regime). Figure 1*a* shows the nullclines *f*(*x*)=*x*(1−*x*)(*x*−*b*)+*d* and *g*(*x*,*a*)=(*x*+*a*)/*c* of the system in the three operating regimes described above, for *a*=−0.1, 0 and 0.06, respectively. In the excitable regime, spikes (also known as pulses) can appear as a result of an external perturbation of large enough amplitude. A convenient definition is that a spike appears when the membrane potential of the neuron exceeds a certain threshold value, e.g. *x*≥0.5. In the oscillatory regime, spikes appear spontaneously with an intrinsic firing frequency *ν* which, as shown in figure 1*b*, does not depend much on the value of *a*.

To illustrate the dynamical behaviour of *x*(*t*) and *y*(*t*), we show in figure 2 the phase portrait for three different values of *a*=−0.1, 0 and 0.06, corresponding to the three nullclines represented in figure 1*a*.

### (b) Coupling scenarios: electrical versus chemical interaction

We consider now that we have an ensemble of *N* coupled neurons. Each one is described by dynamical variables *x*_{i}(*t*), *y*_{i}(*t*), , obeying the FHN equations. The neurons are not isolated from each other, but interact mutually. The full set of equations is now
2.3
and
2.4
where is the coupling term accounting for all the interactions of neuron *i* with other neurons. To take into account the natural diversity of the units, we assume that the parameter *a*_{i}, which controls the degree of excitability of the neuron, varies from one neuron to another. In particular, we assign to *a*_{i} random values drawn from a Gaussian distribution with mean 〈*a*_{i}〉=*a* and correlations 〈(*a*_{i}−*a*)(*a*_{j}−*a*)〉=*δ*_{ij}*σ*^{2}. We use *σ* as a measure of the heterogeneity of the system and, in the following, we use the value of *σ* as an indicator of the degree of diversity. If *σ*=0, all neurons have exactly the same set of parameters, while large values of *σ* indicate a large dispersion in the dynamical properties of individual neurons.

The most common way of communication between neurons is via chemical synapses, where the transmission is carried by an agent called a neurotransmitter. In these synapses, neurons are separated by a synaptic cleft and the neurotransmitter has to diffuse to reach the post-synaptic receptors. In the chemical coupling case, the synaptic current is modelled as
2.5
In this configuration, we consider that neuron *i* is connected to *N*_{c} neurons randomly chosen from the set of *N*−1 available neurons. Once a connection is established between neuron *i* and *j*, we assume that the reciprocal connection is also created. Then, the connection fraction of each neuron is defined as *f*=*N*_{c}/(*N*−1). In equation (2.5), *K* determines the coupling strength and *g*_{ij} represents the maximum conductance in the synapse between the neurons *i* and *j*. For simplicity, we limit our study to the homogeneous coupling configuration, where *g*_{ij}=1 if neurons *i* and *j* are connected and *g*_{ij}=0, otherwise. The character of the synapse is determined by the synaptic reversal potential of the receptor neuron, *E*^{s}_{i}. An excitatory (respectively inhibitory) synapse is characterized by a value of *E*^{s}_{i} greater (respectively smaller) than the membrane resting potential. We consider *E*^{s}_{i}=0.7 for the excitatory synapses and *E*^{s}_{i}=−2.0 for the inhibitory ones. We also define the fraction of excitatory neurons (those that project excitatory synapses) in the system as *f*_{e}=*N*_{e}/*N*; *N*_{e} being the number of excitatory neurons.

Finally, *r*_{j}(*t*) is a time-dependent function representing the fraction of bound receptors and it is given by
2.6
where *α*=2.5 and *β*=3.5 are the rise and decay time constants, respectively. Here, *t*_{on}=0.1 represents the time the synaptic connection remains active and *t* is the time from the spike generation in the presynaptic neuron.

There is another type of synapse where the membranes of the neurons are in close contact and thus the transmission of the signal is achieved directly (electrical synapses). In this case of electrically mediated interactions, also known as diffusive coupling, the total synaptic current is proportional to the sum of the membrane potential difference between a neuron and its neighbours, and it is given by 2.7

The last ingredient of our model is the presence of an external forcing that acts upon all neurons simultaneously. Although our results are very general, for the sake of concreteness, we use a periodic forcing of amplitude *A* and period *T*. More precisely, the dynamical equations under the presence of this forcing are modified as
2.8
and
2.9
which is the basis of our subsequent analysis.

## 3. Results

We are interested in analysing the response of the global system to the external forcing. We will show that its effect can be enhanced under the presence of the right amount of diversity in the set of parameters *a*_{i}, i.e. a conveniently defined response will reach its maximum amplitude at an intermediate value of the r.m.s. value *σ*.

In order to quantify the global response of the system with respect to the diversity, we use the spectral amplification factor defined as 3.1 where is the global, average collective variable of the system and 〈⋅〉 denotes a time average. We will analyse the cases of electrical and chemical coupling separately.

### (a) Electrical coupling

In this subsection, we concentrate on the situation in which the neurons are electrically (diffusively) coupled in a random network, where a neuron *i* is connected randomly with *N*_{c}=*f*(*N*−1) other neurons. The mean value of the Gaussian distribution of the parameters *a*_{i} is fixed to *a*=0.06 and the coupling strength to *K*=0.6. Figure 3 shows the spectral amplification factor *η* as a function of the diversity *σ* for fixed values of the amplitude *A*=0.05 and two different values of the period *T* of the external forcing, for an increasing connection fraction *f*. It can be seen from figure 3 that intermediate values of *σ* give a maximum response in the spectral amplification factor. Moreover, the maximum value shifts slightly to smaller values of *σ* as the fraction of connected neurons *f* increases. We have also observed that a period *T* of the external forcing close to the inverse of the intrinsic firing frequency of the neurons (*ν*≈0.9, according to figure 1*b*) yields the largest response.

In order to further illustrate the response of the system to the external sinusoidal modulation, figure 4 shows the raster plot of the ensemble (lower panels) and the time traces of 10 randomly chosen neurons (upper panels) for different values of the diversity parameter. It can be seen in both the top and bottom panels of this figure that an intermediate level of diversity gives a more regular behaviour than either smaller or larger values of *σ*. This fact is more evident in the time traces, where the amplitude of the oscillation varies randomly for large *σ* values.

### (b) Chemical coupling

We consider in this subsection the situation in which the units are chemically coupled. Figure 5 shows the spectral amplification factor *η* as a function of the diversity *σ* for fixed values of the amplitude *A*=0.05 and two different periods of the external modulation when the fraction of randomly connected neurons *f* increases. The coupling strength is fixed to *K*=1.5. The fraction of excitatory neurons in the system is set to *f*_{e}=0.8. It can be seen from figure 5 that the spectral amplification factor *η* increases as *f* increases, reaching the maximum response at *f*≃0.05. Interestingly, beyond this value *η* does not change significantly, indicating that the response of the system does not improve when the percentage of connected neurons is larger than 5 per cent. Or, put in another way, with a 5 per cent connectivity, the system already behaves as being fully connected as far as the response to the external forcing is concerned. It is also worth noting that the maximum response is always at the same value of *σ*, independent of *f*. The effect of changing the ratio of excitatory/inhibitory synapses in our system is shown in figure 6 in the globally coupled case *f*=1. The spectral amplification factor increases as the fraction of excitatory connections *f*_{e} increases, while the position of the maximum shifts slightly to larger values of *σ*≃0.05.

Comparing the results from both electrical and chemical coupling schemes, it can be seen that the electrical coupling gives a larger value of *η* but requires, at the same time, a larger diversity. The electrical coupling also exhibits a larger range of diversity values for which the system has an optimal response when compared with the chemical coupling. In contrast, the optimal response in the chemical coupling scheme occurs for small values of the diversity and does not significantly change in amplitude and width when the percentage of connected neurons *f* is increased above 5 per cent.

## 4. Order parameter expansion

It is possible to perform an approximate analysis of the effect of the diversity in the case of diffusive (electrical) coupling. The analysis allows one to gain an insight into the amplification mechanism by showing how the effective nullclines of the global variable *X*(*t*) are modified when varying *σ*. The theoretical analysis is based upon a modification of the so-called order parameter expansion developed by Monte & D’Ovidio (2002) and Monte *et al.* (2005) along the lines of Komin & Toral (2010). The approximation begins by expanding the dynamical variables around their average values and as , and the diversity parameter around its mean value, . The validity of this expansion relies on the existence of a coherent behaviour by which the individual units *x*_{j} deviate in a small amount from the global behaviour characterized by the global average variable *X*(*t*). It also assumes that the deviations are small. We expand equations (2.8) for (d*x*_{i}/d*t*) and (2.9) for (d*y*_{i}/d*t*) up to second order in , and ; the resulting equations are
4.1
and
4.2
where
4.3
and we used the notation *f*_{x} to indicate the derivative of *f* with respect to *x*, and so forth.

Note that equation (4.2) is exact as it is linear in all the variables. If we average equations (4.1) and (4.2) using we obtain 4.4 and 4.5 where we have used and defined the second moment . We also defined , , and the shape factors , and . The evolution equations for the second moments are found from the first-order expansion , so that and , were the dot stands for time derivative. 4.6 4.7 4.8 4.9 4.10 4.11 and 4.12

The system of equations (4.4) and (4.5) together with equations (4.8)–(4.12) forms a closed set of differential equations for the average collective variables *X*(*t*) and *Y*(*t*):
4.13
4.14
4.15
4.16
4.17
4.18
and
4.19
Numerical integration of this system allows us to obtain *X*(*t*), from which we can compute the spectral amplification factor *η*. The value of *η* obtained from the expansion is later compared with that obtained from the numerical integration of equations (2.8) and (2.4) (figure 9).

We can obtain another set of closed equations for *X*(*t*) and *Y*(*t*) if we perform an adiabatic elimination of the fluctuations, i.e. , yielding
4.20
with *H*(*x*)=−3*x*^{2}+2(1+*b*)*x*−*b*−*K*. Substituting *Ω*^{x} in equations (4.13) and (4.14), we find a closed form for the equations describing the evolution of the mean-field variables, *X*(*t*) and *Y*(*t*):
4.21
and
4.22
These equations provide a closed form for the nullclines of the global variables *X* and *Y* for the non-forcing case, *A*=0. They also reflect how diversity influences these variables. Figure 7 shows these nullclines *Y*_{1}(*X*,*σ*) and *Y*_{2}(*X*,*a*) of equations (4.21) and (4.22), respectively, for *a*=0.06 and different values of the diversity *σ*. It can be seen in the figure that diversity changes the shape of the cubic nullcline *Y*_{1}(*X*,*σ*) leading to a loss of stability of the fixed point of the system that, for a certain range of *σ*, becomes a limit cycle. To schematize the behaviour of the global variables *X* and *Y* when the diversity changes, we show in figure 8 the phase portrait for different values of *σ*=0, 0.5, 0.8, 1.0, 1.2 and 1.4 (corresponding to the values represented in figure 7). It can be seen that there is a range of *σ* for which the system exhibits a collective oscillatory behaviour even in the absence of the weak external modulation.

With the collective variable *X*(*t*) obtained from the adiabatic elimination, we can estimate the spectral amplification factor *η*. Figure 9 shows the results obtained from the numerical integration of equations (4.13)–(4.19), together with the numerical simulation of the full system (equations (2.8) and (2.4)) and the adiabatic approximation obtained from equations (4.21) and (4.22). It can be seen that our order parameter expansion is in good agreement with the numerical integration of the full system, even in the case in which the second moments are adiabatically eliminated.

## 5. Conclusions

We have studied the effect of diversity in an ensemble of coupled neurons described by the FHN model. We have observed that an intermediate value of diversity can enhance the response of the system to an external periodic forcing. We have studied both electrical and chemical coupling schemes finding that the electrical coupling induces a larger response of the system to an external weak modulation, as well as existing for a wider range. In contrast, the chemical coupling scheme exhibits a smaller optimal amplitude and narrower range of response, but for smaller values of diversity. We have also found that the response of the system in the electrical coupling scheme strongly depends on the fraction of connected neurons in the system, whereas it does not improve much above a small fraction of connected neurons in the chemical coupling scheme. We have also developed an order parameter expansion whose results are in good agreement with those obtained numerically for the electrically (diffusively) coupled FHN system. By an adiabatic elimination of the fluctuations, we have found a closed form of the effective nullclines of the global collective variables of the system obtaining a simple expression of how the diversity influences the collective variables of the system.

The microscopic mechanism leading the system to a resonant behaviour with the external signal is as follows. In the homogeneous situation, where all the units are identical, the weak external modulation cannot induce any spike in the system. When the diversity increases, a fraction of the neurons enters into the oscillatory regime and, owing to the interactions, pull the other neurons with them. This leads the system to an oscillatory collective behaviour that follows the external signal. For larger values of the diversity, the fraction of neurons in the oscillatory regime decreases and the rest of the neurons offer some resistance to being pulled by the oscillatory ones; thus, the system cannot respond to the external signal anymore. These results suggest that the diversity present in biological systems may have an important role in enhancing the response of the system to the detection of weak signals. A physiological example of this kind of functionality can be found in the cochlear nuclei of the kangaroo rat, where the diversity in the different dynamic characteristics that neurons exhibit plays a positive role in the response of the whole system to different external signals (Moushegian & Rupert 1970). In this sense, we hypothesize here that evolution could have selected the right degree of diversity for the system to work ‘optimally’; i.e. one could argue that being able to optimize the response to a varying environment is indeed beneficial from a Darwinian point of view.

## Acknowledgements

We acknowledge financial support from the following organizations: National Science Foundation (grant DMR- 0702890); G. Harold and Leila Y. Mathers Foundation; European Commission Project GABA (FP6-NEST Contract 043309); EU NoE Biosim (LSHB-CT-2004-005137); and MEC (Spain) and Feder (project FIS2007-60327).

## Footnotes

One contribution of 13 to a Theme Issue ‘Complex dynamics of life at different scales: from genomic to global environmental issues’.

- This journal is © 2010 The Royal Society