Predicting collective dynamics from individual neuron properties.

This is a short introduction and rationale to my PhD subject on organized behaviors in neuronal networks, focusing on collective bursting in neuronal cultures.

Contents:

• Introduction to the aEIF model to describe the neuronal dynamics
• Network topology in cultures
• Synchronization in homogeneous networks
• The mean-field model

## 1. Neuronal dynamics and the aEIF model

### Principle

As described in the previous sessions on neuronal models, we use simplified equations to describe the behavior of neurons. How far should we simplify them? For physicists, the best way is to keep the model 2D, so that you can always represent the evolution of the neuron’s state on a piece of paper, but make it discontinuous so that it can also have a chaotic behavior.

For that reason, we will here use the adaptive Exponential Integrate-and-Fire model (a.k.a. aEIF or AdEx) [Brette2005]. In addition to being 2D, this model also recovers a large number of biologically relevant features compared to the IAF, the Izhikevitch, or GIF models.

The model describes the behavior of a neuron using only two variables: the membrane potential $\tilde{V}$, and a slow current $\tilde{w}$, which conveys the capacity of the neuron to adapt, i.e. to modulate its excitability depending on the input it receives.

### Equations

The two variables are coupled through the following equations:

$$\begin{equation} \left\lbrace \begin{array}{l} \tilde{C}_m \frac{d\tilde{V}}{d\tilde{t}} = - \tilde{g}_L(\tilde{V}-\tilde{E}_L) + \tilde{g}_L \tilde{\Delta}T \exp\left(\frac{\tilde{V}-\tilde{V}{th}}{\tilde{\Delta}_T}\right) - \tilde{w} + \tilde{I}_e + \tilde{I}_s\ \tilde{\tau}_w \frac{d\tilde{w}}{d\tilde{t}} = \tilde{a} (\tilde{V} - \tilde{E}_L) - \tilde{w} \end{array}\right. \end{equation}$$

as in the integrate-and-fire model, the spike is terminated by a reset condition, which reads:

$$\begin{equation} \text{if } \tilde{V} > \tilde{V}_{peak} \quad \left\lbrace \begin{array}{r c l} \tilde{V} &\leftarrow & \tilde{V}_r\ \tilde{w} &\leftarrow & \tilde{w} + \tilde{b} \end{array}\right. \label{eq:AdEx_adim} \end{equation}$$

For simplicity, I will use dimensionless variables in the rest of the session. The equations will thus become:

$$\begin{equation} \left\lbrace \begin{array}{l} \dot{V} = - (V-E_L) + e^{V} - w + I_e + I_s\ \tau_w \dot{w} = a (V - E_L) - w \end{array}\right. \end{equation}$$

See Appendices – Variable change for details.

### Behavior

This model is able to reproduce almost all behaviors observed in electrophysiological recordings, from fast-spiking to intrisically bursting or chattering neurons [Naud2008].

## 2. Neuronal networks in cultures

Now that we have the neuronal model, let us have a look at the structure that will underly the system.\
The neurons are coupled by synapses, so this structure can be represented simply by a directed network, described by an adjacency matrix A, where: $$A_{ij} = s_{ij} \delta(i \rightarrow j)$$ i.e. the matrix entry is non-zero only of neuron $i$ extends an axon towards $j$, in which case it’s value represents the synaptic strength of the connection.

In cultures, the connectivity is much simpler than in the brain, since the structure is planar (at least in the cultures we study). Though the network cannot be considered as fully homogeneous, as what was considered in previous sessions, some studies [Cohen2010] suggest that its heterogeneity is rather limited.

Indeed, it seems that the in-degree connectivity can described by a Gaussian distribution $\mathcal{N}(\overline{k}, \sigma_k)$, with an average degree of the order of 100.

For a given average connectivity $\overline{k}$, the Gaussian in-degree (GID) network has the advantage of allowing us to test the influence of heterogeneity by varying $\sigma_k$.

## 3. Synchronization in homogeneous networks

Let us try to understand simply how synchronization can occur in networks containing oscillators.

### Simulations: what does it look like?

One of the most well-known models for this is the Kuramoto model [Acebron2005].

Let a complex oscillator follow any given trajectory in phase-space and let its phase along the trajectory be $\theta(t)$.

In this model, all that we consider is the phase of the oscillator, and the oscillator will interact with the rest of the world only because of the phase delay that might exist between it and some other oscillator. This is denoted, for any oscillator $i$ among $N$, by the equation:

$$\dot{\theta}i(t) = \omega_i + \frac{K}{N}\sum{j=1}^N sin(\theta_j(t) - \theta_i(t))$$

where $K$ is the non-normalized coupling constant.

Many studies on the subject exist, so we will just look at how synchronization can occur on a simpler system, where one oscillator sees a large periodic stimulus (e.g. a population of synchronized oscillators). This exercise is detailed in [Strogatz1994] section 4.5.

The periodic stimulus has a period $\Omega$ and a phase $\Theta$, hence:

$$\dot{\theta}(t) = \omega + K sin(\Theta - \theta(t))$$

looking at the phase difference $\varphi(t) = \Theta - \theta(t)$, and noting $\Delta\omega = \Omega - \omega$, we get:

$$\dot{\varphi}(t) = \Delta\omega - K sin(\varphi(t))$$

Thus we can attain a permanent phase lock solution for $\dot{\varphi} = 0$ if $\Delta\omega \leq K$.

### Relaxation oscillators

Though the Kuramoto model allows us to understand qualitatively many features of synchronization, it is quite different from what happens in neuronal assemblies, especially since these continuous phase models are:

• characterized by rather slow synchronizations, whereas the discrete coupling and almost discontinuous dynamics of neurons leads to much faster synchronization phenomena,
• not leading to zero phase-lag for finite coupling which is very different from the behavior observed in neuronal populations.

As can be seen on Figure 1, when looking at the adaptation variable $w$, bursting neurons display a behaviour which is very similar to that of a relaxation oscillator.

Contrary to the phase oscillator, the relaxation oscillator can have a discontinuous trajectory. This feature is what allows not only the fast transition from the asynchronous to the synchronous phase, but also the existence of a zero phase-lag state for finite couplings.

## 4. The analytic mean-field model

Our model is based on the hypothesis that many neuronal cultures contain intrisically oscillating neurons (the majority of hte pyramidal neurons according to [Penn2016]), and that the collective bursting behavior is observed even in the absence of inhibitory neurons.

Therefore, we use the AdEx model to describe periodic spiking adaptive neurons and tried to understand how adaptation alone could explain the existence of the collective dynamics.

To that end, we set the persistent current $I_e$ to a positive value (mimicking biological persistent sodium currents), so that the neurons spike periodically even in the absence of external input.

The relaxation oscillations of $w$ can be understood as follows: the adaptation variable increases rapidly during the burst towards a peak value $w^∗$, then undergoes a quasi-exponential decrease until it reaches its minimum value $w_{min}$.\
Since the increase of $w$ is simply due to the spike-driven adaptation linked to the $b$ increment after a spike, all that’s left to understand is what determines the minimal and maximal values for $w$.

Figure 2 shows the origin of the phenomena: as the persistent current $I_e$ progressively drives the depolarisation of the membrane, $w$ decreases while $V$ increases. Once $w$ goes below the $V$-nullcline, the potential can diverge and the first spike occurs: the burst begins.

During the burst, the neuron receives an average input from the other neurons, which drives a subsequent series of burst. This phenomenon comes to an end when $w$ reaches a threshold value given by the position at which the spike trajectory (red line + arrow) will intersect the $V$-nullcline (red cross), thus preventing the spike divergence to occur and leading to the termination of the burst.

Using a detailed mean-field model, we could translate this condition to a dimensionless self-coherent equation:

$$\begin{equation} w^* = w_{min} + b \left[ \overline{t_s}(w^*) - d \right] + \overline{k} Q_s. \end{equation}$$

where $\overline{t_s}(w^*)$ is the average interspike during a burst, $d$ is the propagation delay of an action potential between the source and target neuron and $Q_s$ is the charge delivered by a spike on the post-synaptic neuron.

As $\overline{t_s}(w)$ is an increasing function of $w$, we understand that the spike-driven adaptation, described by $b$, will have a significant importance on the properties of the collective dynamics.

**[Acebron2005]** J. A. Acebron, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena

**[Brette2005]** R. Brette and W. Gerstner, Adaptive Exponential Integrate-and-Fire Model as an Effective Description of Neuronal Activity

**[Cohen2010]** O. Cohen, A. Keselman, E. Moses, M. R. Martínez, J. Soriano, and T. Tlusty, Quorum Percolation in Living Neural Networks

**[Ladenbauer2013]** J. Ladenbauer, M. Augustin, L. Shiau, and K. Obermayer, Impact of adaptation currents on synchronization of coupled exponential integrate-and-fire neurons

**[Naud2008]** R. Naud, N. Marcille, C. Clopath, and W. Gerstner, Firing patterns in the adaptive exponential integrate-and-fire model

**[Penn2016]** Y. Penn, M. Segal, and E. Moses, Network synchronization in hippocampal neurons

**[Strogatz1994]** S. Strogatz, Nonlinear Dynamics and Chaos](http://www.stevenstrogatz.com/books/nonlinear-dynamics-and-chaos-with-applications-to-physics-biology-chemistry-and-engineering)

## Appendices

### Variable changes for the neuronal model

$$\begin{eqnarray*} V = \frac{\tilde{V}-\tilde{V}_{th}}{\tilde{\Delta}_T}, & \quad E_L = \frac{\tilde{E}_L - \tilde{V}_{th}}{\tilde{\Delta}_T} \quad & \text{(general relation for voltages)}\ w = \frac{w}{\tilde{g}_L \tilde{\Delta}_T}, & I = \frac{\tilde{I}}{\tilde{g}_L \tilde{\Delta}_T} & \text{(general relation for currents)}\ t = \frac{\tilde{t}}{\tilde{\tau}_m}, & \tau_w = \frac{\tilde{\tau}_w}{\tilde{\tau}_m} & \text{(general relation for times)}\ g_L = 1 & a = \frac{\tilde{a}}{\tilde{g}_L} & \text{(general relation for conductances)} \end{eqnarray*}$$