Simple model of complex dynamics of activity patterns in developing networks of neuronal cultures

TL;DR

This paper introduces a simple, biologically meaningful model of neuronal network activity that explains complex patterns like avalanches and bursts, supported by empirical data and capturing the balance of energy and activity in developing cultures.

Contribution

The paper presents a minimal two-variable model that accounts for complex neuronal activity patterns and self-organized criticality in developing neuronal networks.

Findings

Model explains emergence of network spikes and bursts.

Supports empirical observations of energy-dependent activity.

Networks balance on the edge of percolation transition.

Abstract

Living neuronal networks in dissociated neuronal cultures are widely known for their ability to generate highly robust spatiotemporal activity patterns in various experimental conditions. These include neuronal avalanches satisfying the power scaling law and thereby exemplifying self-organized criticality in living systems. A crucial question is how these patterns can be explained and modeled in a way that is biologically meaningful, mathematically tractable and yet broad enough to account for neuronal heterogeneity and complexity. Here we propose a simple model which may offer an answer to this question. Our derivations are based on just few phenomenological observations concerning input-output behavior of an isolated neuron. A distinctive feature of the model is that at the simplest level of description it comprises of only two variables, a network activity variable and an exogenous variable corresponding to energy needed to sustain the activity and modulate the efficacy of signal transmission. Strikingly, this simple model is already capable of explaining emergence of network spikes and bursts in developing neuronal cultures. The model behavior and predictions are supported by empirical observations and published experimental evidence on cultured neurons behavior exposed to oxygen and energy deprivation. At the larger, network scale, introduction of the energy-dependent regulatory mechanism enables the network to balance on the edge of the network percolation transition. Network activity in this state shows population bursts satisfying the scaling avalanche conditions. This network state is self-sustainable and represents a balance between global network-wide processes and spontaneous activity of individual elements.

Equations

$$B=\pi R^{2}n.$$

$$q_{t+1}=1-(1-p{)}^{z{q}_t}.$$

$${\overline{z}}_1=\frac{-1}{\log (1-p)}.$$

$$p=\frac{1}{{\overline{z}}_1}+O\left(\frac{1}{{{\overline{z}}_1}^2}\right).$$

$$(1-q_2^{\ast })<e(1-p{)}^z\Rightarrow q_2^{\ast }>1-e(1-p{)}^z.$$

$$\begin{array}{rl}{q}_{t+1}=& 1-(1-p\cdot \sigma (E_t,\underline{E}){)}^{z{q}_t}\\ E_{t+1}=& (1-\epsilon )E_t+\epsilon \overline{E}-r{q}_t\mathcal{H}\left(E_t-r{q}_t\right),\end{array}$$

$$\sigma (E_t,\underline{E})=\frac{1}{2}(\tanh (w{E}_t-\underline{E})+1),$$

$$\{(q,E)|q\in [0,1],E\in [0,\overline{E}]\}$$

$$(q_1^{\ast },E_1^{\ast })=(0,\overline{E}).$$

$$\begin{array}{rl}& (q_2^{\ast },E_2^{\ast }),0<q_2^{\ast }\le \frac{\overline{E}}{r\left(1+\frac{1}{\epsilon }\right)},q_2^{\ast }<1,\\ & q_2^{\ast }=1-{\left(1-p\sigma \left(\overline{E}-\frac{r{q}_2^{\ast }}{\epsilon },\underline{E}\right)\right)}^{q_2^{\ast }z},E_2^{\ast }=\overline{E}-\frac{r{q}_2^{\ast }}{\epsilon }.\end{array}$$

$$\begin{array}{rl}& (q_3^{\ast },E_3^{\ast }),1>q_3^{\ast }>\frac{\overline{E}}{r},\\ & q_3^{\ast }=1-(1-p\sigma (\overline{E},\underline{E}){)}^{q_3^{\ast }z},E_3^{\ast }=\overline{E}.\end{array}$$

$$E^{\ast }={\sigma }^{-1}(q^{\ast },\underline{E})=\frac{1}{w}{\tanh }^{-1}\left(\frac{2}{p}\left(1-(1-q^{\ast }{)}^{\frac{1}{q^{\ast }z}}\right)-1\right)+\frac{\underline{E}}{w}$$

$$E^{\ast }=\overline{E}-r\frac{q^{\ast }}{\epsilon },0<q^{\ast }\le \frac{\overline{E}}{r\left(1+\frac{1}{\epsilon }\right)},q^{\ast }<1$$

$$E^{\ast }=\overline{E},\frac{\overline{E}}{r}<q^{\ast }<1.$$

$${\overline{z}}_2=\frac{-1}{\log (1-p\sigma (\overline{E},\underline{E}))}\text{mbox}and{\overline{z}}_3=\frac{\log \left(1-\frac{\overline{E}}{r}\right)}{\frac{\overline{E}}{r}\log (1-p\sigma (\overline{E},\underline{E}))},\frac{\overline{E}}{r}<1.$$

$$c_{ij}=\left\{\begin{array}[]{ll}1,&\mbox{if there is a link from the}\ i\mbox{-th}\ \mbox{node to the}\ j\mbox{-th}\\ 0,&\mbox{otherwise}.\end{array}\right.$$

$$⟨N_{\mathrm{in}}⟩=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N{c}_{ji},⟨N_{\mathrm{out}}⟩=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N{c}_{ij}.$$

$$q_{i,t+1}$$

$$E_{i,t+1}$$

$$E_{i,t}^{\mathrm{fire}}=\frac{r_1{\sum }_{j=1}^N{c}_{ij}}{⟨N_{\mathrm{out}}⟩}+\frac{r_2{\sum }_{j=1}^N{c}_{ji}{q}_{j,t}}{⟨N_{\mathrm{in}}⟩},$$

$$a_{i,t}$$

$$p_{i,t}^{\mathrm{fire}}$$

$$p=0.01,\underline{E}=2,\overline{E}=4,w=1.5,\epsilon =0.0025,N=625,$$

$$f(q_t)=1-(1-p{)}^{z{q}_t}$$

$$\tilde{q}=f(\tilde{q})={\int }_0^{\tilde{q}}\frac{df}{dq}(q)dq.$$

$$E_{t+1}=(1-\epsilon )E_t+\epsilon \overline{E}-r{q}_t\ge \epsilon (\overline{E}-E_t)\ge 0$$

$$E_{t+1}\le (1-\epsilon )\overline{E}+\epsilon \overline{E}-r{q}_t\le \overline{E}.$$

$$J(0,\overline{E})=\left(\begin{array}[]{cc}-z\log(1-p\sigma(\overline{E},\underline{E}))&0\\ -r&1-\varepsilon\end{array}\right)$$

$$\begin{array}{rl}{q}^{\ast }=& 1-(1-p\sigma (E^{\ast },\underline{E}){)}^{q^{\ast }z}\\ E^{\ast }=& (1-\epsilon )E^{\ast }+\epsilon \overline{E}-r{q}^{\ast }\mathcal{H}(E^{\ast }-r{q}^{\ast }).\end{array}$$

$$\{\begin{array}{rl}{q}^{\ast }=& 1-(1-p\sigma (E^{\ast },\underline{E}){)}^{q^{\ast }z}\\ E^{\ast }=& \overline{E}-\frac{r{q}^{\ast }}{\epsilon },\text{mbox}if{E}^{\ast }\ge r{q}^{\ast },\end{array}\{\begin{array}{rl}{q}^{\ast }=& 1-(1-p\sigma (E^{\ast },\underline{E}){)}^{q^{\ast }z}\\ E^{\ast }=& \overline{E},\text{mbox}if{E}^{\ast }<r{q}^{\ast }.\end{array}$$

$$\{\begin{array}{rl}{E}^{\ast }=& g^{-1}\left(\frac{1}{p}\left(1-(1-q^{\ast }{)}^{\frac{1}{q^{\ast }z}}\right)\right),\\ E^{\ast }=& \overline{E}-\frac{r{q}^{\ast }}{\epsilon },\text{mbox}if{E}^{\ast }\ge r{q}^{\ast },\end{array}\{\begin{array}{rl}{E}^{\ast }=& g^{-1}\left(\frac{1}{p}\left(1-(1-q^{\ast }{)}^{\frac{1}{q^{\ast }z}}\right)\right),\\ E^{\ast }=& \overline{E},\text{mbox}if{E}^{\ast }<r{q}^{\ast }.\end{array}$$

Full text

Simple model of complex dynamics of activity patterns in developing networks of neuronal cultures

Ivan Y. Tyukin1,3,4,*, Dmitriy Iudin1,2, Feodor Iudin1, Tatiana Tyukina4, Victor Kazantsev1,2, Irina Muhina1, Alexander N. Gorban1,4

**1 Nizhny Novgorod State University, 23 Gagarin ave., 603950, Nizhny Novgorod, Russia

2 Institute of Applied Physics of RAS, 46 Uljanov str., 603950, Nizhny Novgorod, Russia

3 Saint-Petersburg State Electrotechnical University (LETI), prof. Popova str. 5, 197376, Saint-Petersburg, Russia

4 University of Leicester, Leicester LE1 7RH, United Kingdom

**

*** I.Tyukin@le.ac.uk**

Abstract

Living neuronal networks in dissociated neuronal cultures are widely known for their ability to generate highly robust spatiotemporal activity patterns in various experimental conditions. These include neuronal avalanches satisfying the power scaling law and thereby exemplifying self-organized criticality in living systems. A crucial question is how these patterns can be explained and modeled in a way that is biologically meaningful, mathematically tractable and yet broad enough to account for neuronal heterogeneity and complexity. Here we propose a simple model which may offer an answer to this question. Our derivations are based on just few phenomenological observations concerning input-output behavior of an isolated neuron. A distinctive feature of the model is that at the simplest level of description it comprises of only two variables, a network activity variable and an exogenous variable corresponding to energy needed to sustain the activity and modulate the efficacy of signal transmission. Strikingly, this simple model is already capable of explaining emergence of network spikes and bursts in developing neuronal cultures. The model behavior and predictions are supported by empirical observations and published experimental evidence on cultured neurons behavior exposed to oxygen and energy deprivation. At the larger, network scale, introduction of the energy-dependent regulatory mechanism enables the network to balance on the edge of the network percolation transition. Network activity in this state shows population bursts satisfying the scaling avalanche conditions. This network state is self-sustainable and represents a balance between global network-wide processes and spontaneous activity of individual elements.

Author Summary

We propose a model of spiking activity in living developing neuronal networks in dissociated neuronal cultures. Such neuronal cultures generate a broad range of complex yet highly robust spatiotemporal activity patterns in various experimental conditions. These patterns are often regarded as neuronal avalanches that satisfy the power scaling law and thereby exemplify self-organized criticality in living systems. In this work we approach an important question of how these patterns can be explained and modeled by just few simple equations and variables of which the phenomenological description make clear biological and physical sense. A distinctive feature of the model is that at the simplest level of description it comprises of only two variables, the network activity variable and an exogenous variable corresponding to energy needed to sustain the activity and modulate the efficacy of signal transmission. Despite its apparent simplicity the model is capable of explaining several common behavioral features of developing neuronal cultures, including emergence of network spikes and bursts as a function of days in development. We show that when the model is extended to the larger, network scale, introduction of the energy regulatory mechanism allows the system to balance on the edge of the network percolation transition. Network activity in this state shows population bursts satisfying the scaling avalanche conditions. This network state is self-sustainable and represents a balance between global network-wide processes and spontaneous activity of individual elements.

Introduction

Exploiting physics’ concepts in life sciences has long been reconginzed as a successful strategy for developing systematic understanding of complex phenomena observed in empirical data. One of the most striking and fashionable illustrations facilitating potential and power of this approach is the well-known concept of self-organized criticality (SOC) – the ability of systems to self-tune to a critical state. Initially proposed as a model for explaining how an abstract system can remain at a critical state in presence of perturbations [2, 20], the concept is now broadly used for describing biological neural networks (see e.g. [8, 6]). It has been shown in [23], [26] that networks whose structural evolution is linked with the node’s dynamcis can exhibit highly robust global SOC-like behavior. What is important is that this behavior can be maintained by simple local adjustment rules.

Examples of SOC-like behavior have been found in experimental studies of neuronal cultures [3, 4, 29]. The cultures grow autonomously and form synaptically coupled networks of living cells. After a period of initial growth and development the cultures start to generate spontaneous activity patterns in the form of population bursts. These bursts were shown to satisfy the power scaling law and hence are often referred as neuronal avalanches [3, 4].

Since then a number of mathematical models have been proposed to capture and analyse mechanisms involved in the generaaion of spontaneous burst of activity in neuronal networks. The spectrum of network’s features linked to emergence of persistent bursts includes, but is not limited to, e.g. re-wiring, delays [12, 13], frequency dependent and spike-timing dependent synaptic plasticity [39, 24, 25]. With regards to the mathematical frameworks describing neuronal avalanches, models of network’s growth [1] and stochastic networks [5] have been proposed.

It has been shown in [27] (see also references therein) that population spikes and bursts can be attributed to cell’s adaptation and short time plasticity mechanisms. The authors showed that population spikes similar to the ones observed in in-vitro cultures can occur in networks of excitatory model neurons with leaky integrate and fire dynamics. These neurons where subjected to Gaussian white noise and equipped with adaptation and short term plasticity mechanism. The network connectivity was all-to-all. The influence of neuronal connectivity on bursts was studied in [19]. It has been found that the number of synaptic connections per neuron, may play an important role for spiking and bursting activity in cultures. Additional links between connectivity development, firing activity homeostasis, and criticality are exposed in [37].

In this work we further contribute to the idea that several features of complex and critical behavior (e.g. the neuronal avalanches, super-bursts, periodic and chaotic spiking) observed in live neuronal cultures and networks can be explained by just few variables. These variables can be linked to local connectivity patterns (expressed e.g. by connection densities between cells) and neuronal activation dynamics.

We demonstrate that main critical transitions can be captured by a hierarchy of simple models. Starting from elementary phenomenological description of neural firing, we present a simple $2D$ mean-field model that is capable of a broad range of behaviors conistent with the ones observed in cultures. Remarkably, the network connectivity parameter appears to be a natural bifurcation parameter of the model; it regulates emergence of activity spikes from an initially silent mode to occasional network bursts which, in turn, eventually develop into a whole-network activation. This is in good agreement with previous works on SOC phenomena in developing neuronal networks [26], [37]. In contrast to [26], [37], however, we consider the problem from a different angle. Instead of focusing on the activity-connectivity interplay we consider and analyze the system’s dynamics in the activity-energy plane for various values of connectivity. This allows additional modelling capabilities, including e.g. the analysis of how oxygen and energy deprivation affects the network’s behavior and dynamics. Moving further to multi-agent model reveals emergence of neuronal avalanches showing power law scaling.

The manuscript is organized as follows. In  Results we present ingredients of the model at three different levels of phenomenological detail. We begin with a simple percolation-based geometric model describing the evolution of cells’ connectivity is presented. The model allows to accommodate biologically relevant features such as axons and dendrites; it also enables to replicate directional connectivity that is inherent to living systems including neuronal cultures. The model analysis reveals that sharp changes in the overall clustering and connectivity of the evolving network in both directed and undirected settings is determined by a single parameter describing average connection density in the network. The analysis is followed by expressing dynamics of neuronal activity by a mean-field approximation. We show that the corresponding single-dimensional model does not explain network spikes and bursts frequently observed in developing cultures. This limitation, however, can be resolved if neural activation is linked to an additional exogenous regulatory energy variable. Introduction of the latter variable needs an additional comment. It behaves as a sort of “energy” or a resource. Its physical nature, nevertheless, may or may not be associated with a specific type of the physical energy. A prototype of such energy, the notion of adaptation energy, was introduced by Selye in his analysis of physiological adaption [34], [35] and was successfully employed in modelling of various complex phenomena [17], [18]. The extended model is shown to be able to reproduce periodic spiking, irregular dynamics, and population bursts observed in live cultures. What is important is that dynamic regimes exhibited by the model can be regulated by just a single parameter corresponding to network connectivity. Next we provide results of large-scale simulation of evolving network of agents of which the activation probability depends on their current energy level. The network dynamics in this state shows population bursts satisfying the scaling avalanche conditions. This network state is self-sustainable and represents energetic balance between global network-wide processes and spontaneous activity of individual elements. The results are compared with empirical data. Section  Materials and Methods describes details of electrophysiological measurements and experiments, and Section Conclusion concludes the paper.

Results

Model

Geometric model

We start with geometrical arrangement of the network elements. Consider a network of $N$ neurons whose spatial coordinates are randomly and uniformly distributed in the unit square. Each individual neuron is described by two basic elements. The first is the region of reception of inputs represented by a circle of a given radius $R$. The circle models neuron’s ability to sense input signals from other neurons, and is referred to as the dendrite region (in biology, dendrite is an input). The second element is an axon (in biology, output), which in our model is simulated by a straight segment of length $H$ (on the mature stage of the network development $H>R$) and whose end point is acting as a transmitter of the neuron’s signal. If this point reaches out to the dendrite region of another neuron, a connection is established between these neurons [1]. There are three different ways that yield geometrical coupling or connectivity of the network elements:

Case 1: cells without axons, i.e. $H=0$. In this case $N$ circles with radius $R$ are randomly and uniformly distributed in the unit square. If a circle A overlaps with a circle B, and circle B is connected with a circle C, then A is connected with C. Thus, a path between two distant cells can be defined as a chain of overlapping circles joining these cells. Emergence of large groups of connected elements in this network can be analyzed within the framework of standard circle percolation problem. Let $n$ be the cells density defined e.g. as the number of circles’ centers in a unit area. According to [33, 36], emergence of large groups of interconnected cells, the percolation transition, in a set of randomly distributed circles can be characterised by the mean number of centers that fall within a circle of radius $R$:

[TABLE]

In particular, there exists a critical concentration $B=B_{c}$ at which two arbitrary circles become connected with high probability. Thus percolation occurs and a large cluster of connected circles appears. In contrast with typical thermal phase transitions where a transition between two phases occurs at a critical temperature, percolation transition relates to distribution and topology of clusters corresponding to the values of $B$ in a neighborhood of $B_{c}$. At low values of $B$ only small clusters of overlapping circles exist. When the concentration $B$ increases the average size of the clusters increases too. At the critical concentration value, $B=B_{c}$, a large cluster appears; it includes groups of cells that are close to the opposing boundaries of the original square. This cluster is called spanning cluster or percolating cluster. In the thermodynamic limit, i.e. in the infinite size system limit the spanning cluster is called infinite cluster. For scalar problem the value of $B_{c}\approx 1.1$.

Case 2. Cells have axons, $H>0$, and axons are allowed to transmit signals in both directions. Each neuron can be represented as an undirected pair of head- and tail-circles both having radius $R$. When the head-circle or the tail-circle of an neuron overlaps with the head- or the tail-circles of another neuron we consider these neurons connected. Despite obvious difference of this setting from the previously considered one in that we are to operate with dipoles here rather than with just circles as in Case 1, the problem remains within a class of scalar percolation, albeit for dipoles of circles not just a single circle.

Case 3. Cells have axons, $H>0$, and these axon can only transmit signals along a straight line that determines direction of connectivity for a given cell. The coupling direction from soma to synaptic terminal has isotropic distribution, and hence each neuron could be represented as a directed pair of head- and tail-circles both having radius $R$. Vectors linking centers of the head- and tail-circles are allowed to have arbitrary direction. Their lengths, $H$, however, are fixed. When the tail-circle of a neuron overlaps with a head-circle of another neuron the pair is considered as connected. In contrast with two other ways of establishing neuronal connectivity considered above this is the most realistic scenario. It is no longer within the scope of simple scalar circle percolation framework but is a vector percolation problem.

The three cases are illustrated with Fig. 1. Fig. 2 shows dependence of the percolation threshold parameter, $B_{c}$, on the ratio $H/R$.

In accordance with he definition of $B$ in (1) the cell’s concentration variable $B$ can be related to the expected number of neurons that are connected to a randomly chosen query neuron in the system. The latter quantity is known as the average network coordination number $z$. As can be seen from Fig. 2 the value of the coordination number $z$ that corresponds to $B_{c}$ is always bounded from above. $B_{c}\leq 1.4$ for all configurations considered so far.

The above analysis reveals that networks with coordination numbers exceeding these critical values are likely to form a spanning cluster that is capable of connecting opposite edges of the system [21, 22]. Thus, intuitively, one can argue that signals initiated by spontaneous activation of neurons in the cluster can spark waves of activity through the whole network. Such conclusion is based on the restrictive assumption that neurons always elicit spikes in response to a spike on their input. Not only this assumption does not necessarily hold true, but also the above geometrical model alone does not explain the wealth of excitation propagation phenomena observed in cultures. To account for more plausible situations as well as to possibly increase explanatory power of the model two additional variables are introduced: one is the maximal probability $p$ of neuronal activation in response to incoming spike, and the other is an exogenous “resource” variable $E$ determining if a neuron has enough energy to elicit a spike. The consequences of adding these two variables are discussed in the next sections.

Mean-field dynamic model of neuronal excitation

We begin with a simple mean-field approximation of the dynamics of neural excitation in the system. Consider a connected network of neuronal cells. Let $z$ be the expected coordination number, i.e. the expected number of neighbors of a randomly chosen query cell. Suppose that at a time instance $t$ some neurons in the network are excited. Let $q_{t}$ denote the number of these neurons relative to the total number of cells the whole network. If the value of $z$ is sufficiently large then the number of excited neurons among all neighbors of a given neuron can be estimated as $zq_{t}$. Let $p$ be the probability of neuronal activation in response to activation of at least one neuron from its nearest neighbours. Moreover, we suppose that all excitatory signals are independent. Thus the probability that a given neuron is activated equals to $1-(1-p)^{zq_{t}}$, and hence the expected proportion of all excited neurons at the time step $t+1$ is:

[TABLE]

The range of dynamics which model (2) is capable to reproduce is summarized in Proposition 1.

Proposition 1

Consider (2) and suppose that $p\in(0,1)$, $z>0$ be constants. Then the interval $[0,1]$ is forward-invariant, all forward orbits $q_{t}$ are monotone functions $t$, and the point map (2) has only fixed points as attractors. Furthermore

if $-z\log(1-p)\leq 1$ then the map has only one fixed point, $q^{\ast}_{1}=0$, and it is an attractor; 2. 2.

if $-z\log(1-p)>1$ then the map has only two fixed points with $q^{\ast}_{1}=0$ being a repeller and the other one, $q^{\ast}_{2}\in(0,1)$, a stable attractor.

Proof of Proposition 1 is provided in the Appendix.

According to Proposition 1, asymptotic mean-field dynamics (3) of the network is a steady activity at an equilibrium in the entire domain of the model’s feasible parameters: $p\in(0,1)$, $z>0$. Moreover, all transients are monotone trajectories. Emergence of the unique non-zero asymptotic activity is fully determined by the values of the connectivity parameter, $z$, and the probability of neural activation, $p$. Critical values of these parameters, e.g. the critical connectivity, $\overline{z}_{1}$, at which the transition occurs, satisfies:

[TABLE]

For $\overline{z}_{1}\gg 1$, this relation is approximately reciprocal: $p\simeq 1/\overline{z}_{1}$. Indeed, expressing $p=1-\exp\left(-\frac{1}{\overline{z}_{1}}\right)$ from (3), and expanding $\exp\left(-\frac{1}{\overline{z}_{1}}\right)$ as a power series with respect to $1/\overline{z}_{1}$ one obtains:

[TABLE]

The value of the stable equilibrium, $q^{\ast}_{2}$, can be estimated as follows. At the steady state we have that $1-q^{\ast}_{2}=(1-p)^{zq^{\ast}_{2}}$. Hence $(1-q^{\ast}_{2})^{1/q^{\ast}_{2}}=(1-p)^{z}$. It is well known that $(1-x)e^{-1}<(1-x)^{1/x}$ for all $x\in(0,1)$. Thus

[TABLE]

Observe that the larger is the value of the connectivity parameter, $z$, the higher is the level of the mean-field network activity.

Whilst model (3) is consistent with the very basic observation that increasing the network’s overall connectivity may lead to emergence of a self-sustained activity in the network, the model’s explanatory capability is limited. The model does not explain widely-reported richness of the dynamics in live neuronal cultures, including emergence of spontaneous activity bursts and irregular and seemingly chaotic spikes.

This limitation is not surprising since (2) is a crude approximation of the network’s dynamics. Model (2) does not account for a broad spectrum of biological mechanisms involved in spike generation and assumes that the neuron’s ability to produce spikes depends exclusively on stimulation. A possible way to overcome this unrealistic assumption is to explicitly account for these missing mechanisms. To keep the model simple, we account for joint effect of these mechanisms by adding a single energy-like variable $E_{t}$ to (2). The new variable determines the neuron/network’s ability to produce a spike depending on the amount of resources or “energy” available. Generic models of this type have been proposed and analyzed in [16, 17] in the context of adaption to stress and external environmental factors. These models have been shown to capture periodic and irregular behavior in multi-agent systems [18].

Here we extend the original phenomenological mean-field model (2) as follows:

[TABLE]

where

[TABLE]

and $\mathcal{H}$ is the Heaviside function. In (4) $p$ is the maximal probability of neuronal activation, $z>0$ is the coordination number, $E_{t}$ is the exogenous phenomenological “energy resource” variable; $r>0$ is the energy cost of neuronal activation, $\underline{E}>0$ and $w>0$ are parameters that determine the minimal activation probability and the energy activation threshold, $\overline{E}>\underline{E}$ is the energy recovery value, $\varepsilon\in(0,1)$ is the energy relaxation parameter. General shape of the function $\sigma(\cdot,\underline{E})$ in the energy-dependent synaptic efficacy component, $p\sigma(\cdot,\underline{E})$, is shown in Fig. 3.

Mean-field model (4) of the network dynamics inherits phenomenological transparency of (3). It does, however, account for generic constraints of spike-generation through the new energy variable $E_{t}$ and energy-modulated synaptic efficacy $p\sigma(\cdot,\underline{E})$. Despite retaining simplicity, the model produces remarkably rich dynamics. Its equilibria, however, can still be described by just few parameters as follows from Proposition 2 below.

Proposition 2

Consider (4) with $p\in(0,1)$. Then the domain

[TABLE]

is forward-invariant. In addition,

If $-z\log\left(1-p\sigma\left(\overline{E},\underline{E}\right)\right)\leq 1$ then (4) has only one fixed point:

[TABLE]

This fixed point is an attractor when $-z\log\left(1-p\sigma\left(\overline{E},\underline{E}\right)\right)<1$. 2. 2.

If $-z\log\left(1-p\sigma\left(\overline{E},\underline{E}\right)\right)>1$ and ${\overline{E}}/{r}\geq 1$ then (4) has at most two fixed points: fixed point (5) and, if exists, an additional fixed point $(q^{\ast}_{2},E^{\ast}_{2})$:

[TABLE]

The fixed point $(q^{\ast}_{1},E^{\ast}_{1})$ a repeller. 3. 3.

If $-z\log\left(1-p\sigma\left(\overline{E},\underline{E}\right)\right)>1$ and ${\overline{E}}/{r}<1$ then (4) has at most three fixed points. The fixed point $(q^{\ast}_{1},E^{\ast}_{1})$ specified by (5) and, possibly, additional two fixed points: the fixed point $(q^{\ast}_{2},E^{\ast}_{2})$ specified by (6) and a fixed point

[TABLE]

The fixed point $(q^{\ast}_{1},E^{\ast}_{1})$ is a repeller, and $(q^{\ast}_{3},E^{\ast}_{3})$, if exists, is a stable attractor.

Proof of Proposition 2 is provided in the Appendix.

An illustration showing relationships between parameters of the model and emergence of the three different equilibria described in Proposition 2 is provided in Fig. 4. The equilibria are shown as white circles. Green lines show curves

[TABLE]

as functions of $q^{\ast}>0$. According to the first equation of (4), all equilibria of the model with $q^{\ast}\neq 0$ must belong to these curves (see also the proof of Proposition 2 in Appendix). Depending on the value of $z$, the curves move up and down, and intersect with line segments (shown as red solid lines in Fig. 4):

[TABLE]

and

[TABLE]

These intersections correspond to equilibria (6) and (7), respectively.

The equilibria persist over intervals of $z$, and the greatest lower bounds of these intervals (critical values of $z$) are:

[TABLE]

We also note (see the proof of Proposition 2) that equilibria (6) are always above or on the line $E=rq$ (dashed orange line in Fig. 4), whereas equilibria (7) are to be below this line.

According to Propositions 1, 2, models (2) and (4) share some similarity. For $z$ sufficiently small all orbits are attracted to a single equilibrium. At this equilibrium, the systems are silent. When $z$ increases and eventually exceeds the first critical value (eq. (3) for (2) and $\overline{z}_{2}$ for (4)), the silent equilibrium becomes a repeller and the systems start to exhibit non-zero activity. However, further increases of $z$ trigger drastically different dynamics in these models.

All orbits of model (2) with $q_{0}\neq 0$, as ensured by Proposition 1, converge monotonically to a single non-zero steady state regardless of how large the values of $z$ become. The spectrum of orbits in model (4) is different. Our numerical experiments demonstrated that, in addition equilibria, the model is capable of generating periodic orbits too. Moreover, for a broad range of parameters it produced complicated and apparently chaotic motions. Examples of these complicated motions are shown in Fig. 5.

Observe that the model parameters corresponding to the trajectories in Fig. 5 satisfy statement 2) of Proposition 2. In this case, at most two equilibria may exist. As we can see from Fig. 5, these equilibria (fixed points (5) and (6)) are not attracting the orbits, and trajectories appear to be chaotic with some apparent intermittency.

In order to gain additional insight into the model’s dynamics, we numerically explored asymptotic regimes of (4) for varying values of $z$, $\underline{E}$, and $r$. Other parameter were as follows: $\varepsilon=0.05$, $w=1.5$, $\overline{E}=4$, $p=0.1$. Outcomes of these experiments are summarized in Figs. 6 – 8 (see Materials and Methods for details of the steps taken to produce these figures). In these experiments, the values of $\underline{E}$ were chosen from a uniform equispaced grid of $21$ points in the interval $[1,3]$. This grid is shown as grey dashed horizontal lines in Figs. 6 – 8. Parameter $z$ was varying adaptively (increments ranged from $0.1$ in the intervals $(0,50]$ and $(200,300]$ to $5$ in the interval $(50,200]$). For these values of model parameters, we assessed the type of the model’s asymptotic dynamics and mapped these onto relevant parametric regions. These regions are shown with different colour in Figs. 6 – 8.

According to Figs. 6–8, development and evolution of the dynamics of (4) follows a robust pattern. For a fixed value of $\underline{E}\in[1,3]$ and $z$ small, trajectories of the model converge to a unique attractor (fixed point 5). This attractor corresponds to the system’s state in which no elements/neurons are excited. When the value of $z$ increases and exceeds $\overline{z}_{2}$ (specified in (9) and shown as blue dashed lines in 6–8), equilibrium (5) becomes a repeller and a second equilibrium emerges (fixed point (6). Numerical evaluation of the eigenvalues of the Jacobian at this equilibrium showed that it is locally asymptotically stable. Further increases of $z$ lead to that fixed point (6) loses stability thorough the Neimark-Sacker bifurcation, and an attracting periodic orbit emerges. The boundary of this transition is depicted as green dashed lines in Figs. 6–8. As we keep increasing the values of $z$, non-trivial and complex dynamics eventually occur (red area in Figs. 6–8). Complex orbits and behavior persist over intervals of values of $z$. Notice that some of the complex trajectories shown in panel $C$, 8 eventually converge to the stable equilibrium specified by (7). This is an empirical manifestation of slow relaxations and critical delays in model (4) [14], [15]. For $z$ sufficiently large, these complex orbits disappear and reduce to periodic orbits, Figs. 6, 7, or mere equilibria, Fig. 8.

Mean-field bursting dynamics, shown e.g. in panels D and E in Figs. 6, 7, resembles that of the population bursts observed in live neuronal evolving cultures. An important factor in successful replication of this behavior was the energy variable, $E_{t}$, coupled with the energy-dependent activation probability $p\sigma(E_{t},\underline{E})$. The mean-field model, however, does not capture spatial effects and as such is only a rough approximation of activity propagation in neuronal cultures. In the next section we extend the proposed mean-filed model (4) to a multi-agent network with randomized energy-dependent activation and numerically assess relevant parameters of its dynamics, including distributions of sizes and durations of firing avalanches.

Multi-agent model of neuronal excitation

As a natural extension of (4) we considered a connected network of $N$ neurons of which the activity is governed by exogenous energy variables. The network’s topology is defined by its adjacency matrix, $C$, whose elements $c_{ij}$ are:

[TABLE]

No links from a node to itself are permitted, but cycles are allowed. For simplicity, all links in the network have been assigned equal weights of which the value was assumed to be $1$. For the given adjacency matrix $C$, we determined the average number of inputs, $\langle{N}_{\mathrm{in}}\rangle$, and the average number of outputs $\langle{N}_{\mathrm{out}}\rangle$

[TABLE]

Each $i$-th node in the network is described by two variables: the activity variable, $q_{i,t}\in\mathbb{R}$, and the energy variable, $E_{i,t}\in\mathbb{R}$. Dynamics of these variables is defined as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

The variables $q_{i,t}$ take values in the set $\{0,1\}$, and $E_{i,t}$ are in the interval $[0,\overline{E}]$. The function $\sigma(\cdot,\underline{E})$ is as in (4).

Phenomenological motivation for the dynamics of individual nodes in model (10) is similar to that of the mean-field model, (4). There are, however, several key differences. The evolution of variables $q_{i,t}$ and $E_{i,t}$ in (10) depends explicitly on the network topology and activity of the node’s neighboring cells. The energy balance equation, the second equation in (10), accounts for the costs of transmitting active signals at the neuron’s input (term ${r_{2}\sum_{j=1}^{N}c_{ji}q_{j,t}}{\langle{N}_{\mathrm{in}}\rangle}^{-1}$) and generating activity signals on the neuron’s output (term ${r_{1}\sum_{j=1}^{N}c_{ij}}{\langle{N}_{\mathrm{out}}\rangle}^{-1}$). If the node’s energy level is insufficient to trigger a spike, $E_{i,t}\leq E^{\mathrm{fire}}_{i,t}$, then no spikes will generated. The latter property is difficult to fully capture at the level of the mean-field approximation, as low sub-threshold values of the bulk energy do not necessarily imply absence of activity at the level of individual neurons (cf. Proposition 2, alternative $3$, and Fig. 4).

In our numerical experiments we focused on fully connected networks for which $c_{ij}=1-\delta_{ij}$, where $\delta_{ij}$ is the Kronecker’s delta. We also observed that adding a fraction of inhibitory connections does not altar the network dynamics qualitatively. These simplifications are consistent with the approaches used in earlier works [26]. The model parameters where set as follows:

[TABLE]

and parameters $r_{1}$ and $r_{2}$ varied in the intervals $[1,1.5]$ and $[4,6]$, respectively.

We simulated forward orbits of model (10) for various initial conditions and parameter values, and determined sizes and durations of avalanches of firing events. In our experiments, the avalanches were defined as events corresponding to the intervals $T_{j}=[t_{j},t_{j+1}]$ of the network nonzero firing activity such that $\sum_{i=1}^{N}q_{i,t}>0$ for all $t\in T_{j}$ and $\sum_{i=1}^{N}q_{i,t_{j-1}}=\sum_{i=1}^{N}q_{i,t_{j+1}+1}=0$. Each orbit was simulated for $10^{6}$ time steps, with $q_{i,0}=0$, $i=1,\dots,N$ and $E_{i,0}$, $i=1,\dots,N$ chosen randomly in the interval $[0.5\underline{E},\overline{E}]$. For each orbit, we gathered statistics of sizes and durations of the observed avalanches. A brief summary of these experiments is shown in Fig. 9 and Fig. 10.

Fig. 9 presents size and duration histograms as a function of parameters $r_{1}$ and $r_{2}$ in (10). As we can see from this figure, energy feedback has the capacity to inhibit system-size events in networks with identical graph topologies. Parameters of this feedback may affect the values of size and duration exponents. In particular, we observed that increading the value of $\varepsilon$ facilitates occurances of system-size events and eventually pushes the system into the super-critical state.

Fig. 10 shows statistics and estimated exponents of avalanches for $r_{2}=5.5$, $r_{1}=1.1$.

The estimated exponents are close to those reported for live neuronal cultures [4], [29], [10]. Observe that the size and duration curves have noticeable humps (cf. [10], Fig. 2), albeit different and less prominent exponents than those reported in [10].

Activity bursts appear to be synchronized with the peaks of the energy function. This bears some similarity with network models in which connections change according to rate-dependent synaptic plasticity [39]. Analogous behavior was observed in the dynamics of mean-field approximation (4), Figs. 6–8.

Comparison with empirical data

The model presented above shows a range of behaviors controlled by just few parameters wich may be explicitly linked with relevant physical quantities. Plots presented in Figs. 6 – 8 show evolution of activity patterns in the network as a function of connectivity parameter $z$. It would therefore be nice to see how these observations relate to patterns observed in actual cultures. Fig. 11 shows typical activity of live developing neuronal cultures over time.

We can see from Fig. 11 that, depending on days in development, networks of living cultures exhibit irregular spiking, bursts, and (nearly) periodic spiking. The mean-field dynamics of the activity variable, $q_{t}$, in our model shows qualitatively similar regimes. Plots shown in Figs. 6– 8, however, have not been subjected to any threshold subtraction. Fig. 11, bottom row, shows behavior of the model with threshold subtraction for growing values of $z$. The value of threshold, $\theta$, is set to $1/0.6745$ of the median of $q_{t}$ over the interval of simulation $[0,4000]$. Increases in values of $z$ model development of synaptic connections over time. Evolution of the variable $q_{t}$ resembles that of the culture activity. Spikes of $q_{t}$ in the leftmost plot in are narrow and do not correspond to the network spikes. When the value of $z$ is increased their width grows and they start to appear in packs (middle picture); this corresponds to network bursts shown in the first row. When the value of $z$ is increased further the bursts start to break (the rightmost plot), and when $z$ reaches even larger values no high-amplitude spikes are generated in the model (not shown in the figure). The latter phenomenon was not observed in live cultures. Network connectivity in cultures, however, does not grow indefinitely. Therefore network activity for extreme values of $z$ maye be discarded as that outside of the model validity zone.

Last but not least is that behavior of the model is consistent with the effects of oxygen deprivation in cultures. Fig. 12 shows network activity during and after acute but short oxygen deprivation (see also [41]). In the latter experiments primary cultures of hippocampal cells were subjected to $10$min of acute oxygen deprivation in the $21$st DIV (see Methods). Reoxygenation after short-term hypoxia rapidly restores energy deficit and neuronal ATP levels and increases the release of glutamate. Glutamate is a major excitatory neurotransmitter in mammalian central nervous system. Glutamatergic neurons form the main excitatory system in the brain and play a pivotal role in many neurophysiological functions [43]. Excessive release of glutamate is homeostatic response to the hypoxia-induced network silence [30]. On the one hand, it is aimed at restoration of network activity. On the other hand, it over-activates ATP-dependent ion pumps and changes calcium homeostasis. This in turn leads to a cascade of intracellular events causing neuronal degeneration or excitotoxicity [9], [31].

In order to see how the proposed model might capture hypoxia as well as the sequence of complex biological of changes related to oxygen deprivation the following experiments have been performed. The model, eq (4), with $r=1.5$, $z=170$, $\varepsilon=0.05$, $w=1.5$, $\underline{E}=2$, $\overline{E}=4$, $p=0.1$ was run for $1500$ steps. Then for the next $1000$ steps the value of $\overline{E}$ was set to $0.4$ and then restored to the nominal level of $4$ for $t>2500$. This corresponds to energy deficit caused by acute oxygen deprivation. At $t=2500$, however, the value of $p$ was increased to $0.15$ to account for glutamate release, and then dropped to the level $p=0.07$ in the interval $[2700,5000]$ to emulate glutamate-induced suppression [28]. For $t>5000$ the value of $p$ was made to decrease linearly to account for degenerative processes triggered by oxygen deprivation. Model behaviour as well as the evolution of $p$ and $\overline{E}$ over time are shown in Fig. 12. Overall evolution patterns that the model produces qualtiatively coincide with empirical observation. Notice also that our empirical data showed increase in the level of activity after $2$h of hypoxia. Qualitatively very similar behavior has been reproduced by the model. In this regime, the value of neural excitability parameter $p$ was lower than in the nominal pre-hypoxic state to simulate the combined influence of glutamate-induced suppression. Yet, the network produced larger number of spikes. This can be attributed to network effects and the interplay between activity-energy variables. We also note that other experiments reported e.g. in [11] show that cultures exposed to hypoxia may show reduced activity during oxygen deprivation and after reoxygenation for $2$h (Fig. 5 in [11]). The experimental protocol, however, in [11] differed from ours in that the hypoxia was induced by exposing the cultures to $3\%$ $O_{2}$ for two hours, and that the age of the cultures was $11$ and $18$ DIV. Variability of development and neurons might also have attributed to differences in observed behavior.

Materials and Methods

Cell cultures

Hippocampal cells were dissociated from embryonic mice 57L/6J (on embryonic day $18$) and plated with a high initial density of approximately $9000$ cells/mm2 on microelectrode arrays (Alpha MED Science, Japan) pre-treated with the adhesion promoting molecule polyethyleneimine (Sigma P3143) according the protocol [40]. The cells were cultured under constant conditions of $35.5^{\circ}$C, $5\%$ CO2 and $95\%$ air at saturating humidity in a cell culture incubator (MCO-18AIC, Sanyo) in a medium containing Neurobasal medium (Invitrogen 21103-049) with $2\%$ B27 (Invitrogen 17504-044), $1$ mM L-glutamine (Invitrogen 25030-024) and $0.4\%$ fetal calf serum (PanEco 055) without any antibiotics or antimycotics. Glial growth was not suppressed because glial cells are essential for long-term culture maintenance. One half of the medium was changed every $2$ days.

Electrophysiological methods

Extracellular potentials were recorded simultaneously through $64$ planar indium tin-oxide (ITO) platinum black electrodes with the integrated MED64 system (Alpha MED Science, Japan). MEAs were $8\times 8$ $(64)$ with a $50$ $\mu$m $\times$ $50$ $\mu$m electrode size and a $150\mu m$ spacing, and the sampling rate was $20$ kHz/channel. All of the signal analyses and statistics were performed using custom made software (Matlab$\circledR$).

Spike detection

Detection of recorded extracellular spikes was based on threshold calculation using the signal median: $T=N_{s}\sigma$, $\sigma={\mathrm{median}}\left(|x|/0.6745\right)$, where $x$ is the bandpass-filtered ($0.3-8$ KHz) data signal, $\sigma$ is an estimate of the standard deviation of the signal with zero spikes, and $N_{s}$ is a spike detection coefficient determining the detection threshold [32]. The standard deviation of Gaussian noise is equal to the median of the absolute values of the signal divided by $0.6745$. $N_{s}=4$ was used for all data, resulting in a reliable detection of spikes with amplitudes greater than $20$ $\mu$V. The minimal interspike interval was set to $1$ ms. Detected spikes were plotted in raster diagrams.

To detect small network bursts, we calculated the total spiking rate (TSR) accounting for the total number of spikes from all electrodes within 50 ms time bins. A fast appearance of a large number of spikes over all of the electrodes in a small ($50$ ms) time bin was used as a criterion for small burst appearance (see [32] for more details).

Superbursts in the electrical activity recorded from the multielectrode arrays were detected as follows. First, we defined a Gaussian function with an effective width equal to $10$ s. Next, we iteratively moved that function from the beginning of the recording to the end in $10$ ms time steps and calculated a cross-correlation of the function with the TSR. The resulting cross-correlation indicated how much of the synchronized activity (bursts) was recorded in the $10$ s window. The procedure was identical to the one reported in [40]. To detect superbursts in the spiking activity, we applied threshold detection, in which the threshold was estimated as the spiking superburst detection accuracy coefficient multiplied by the standard deviation of the calculated cross-correlation. The superburst detection accuracy coefficient was found empirically and was equal to $0.4$. All time points that crossed the threshold were defined as the beginnings and endings of the superbursts.

Construction of Figures 6 – 8

To construct the figures, we run $3\cdot 10^{5}$ iterations of model (4) from $5$ random initial conditions for each relevant set of parameters $\underline{E}$, $z$, $r$, whilst keeping the values of $w,p,\overline{E}$, and $\varepsilon$ constant ($w=1.5$, $p=0.1$, $\overline{E}=4$, $\varepsilon=0.05$). The values of parameter $\underline{E}$ where chosen in the equispaced grid of $21$ points in the interval $[1,3]$. The values of $z$ were varying adaptively. In the intervals $(0,50]$ and $(200,300]$ these values were taken from equispaced grids with distances between the grid’s nodes being equal to $0.1$. In the interval $(50,200]$ these distances were set to $5$. For each set of model parameters and each initial condition, last $2\cdot 10^{4}$ points in each run have been recorded. Points corresponding to orbits from different initial conditions have been collated together (color-coded), plotted in $(q_{t},E_{t})$ space and stored as .gif files at [38]. This resulted in circa $32000$ images for $r=1$ and $r=1.5$, and circa $10000$ images for $r=10$. The resulting figures were visually inspected and classified into orbits converging to a) single equilibrium (5), b) single equilibrium (6) c) single periodic orbit, d) complicated sets like the ones shown in Fig. 5, and e) multiple attractors. These were color-coded and mapped onto the relevant parametric domains.

Hypoxia modelling

Hypoxia modelling was conducted by replacing the normoxic culture medium with a medium containing low oxygen (0.37mL/L) for 10 min on day $21$ of culture development in vitro according to the previously described protocol [42].

Conclusion

In summary, we have proposed a simple network model explaining burst generation in living culture networks. A distinct feature of our model is presence of a dynamic exogenous energy variable and neuronal activation probability that is made dependent on the energy, like in general models of physiological adaptation [18]. We showed that introduction of these modifications already enables to explain evolution of cultures from resting state to population bursts, at least in the mean-field approximation. In accordance to the model, emergence of bursts and spikes is regulated by just few parameters that correspond to network connectivity and efficacy of synaptic transmission. We also note that our energy-based model is complementary to more traditional connectivity-focused approaches [26].

In this particular study, when comparing empirical data with model behavior, the number of days in development has been related to network connectivity. We note, however, that the latter in a broader biological context can depend on various external factors such as e.g. stress [7]. The proposed model hence might be able to predict qualitatively the effect of stress and adaptation to stress in neuronal activity.

Large-scale multi-agent simulations demonstrated that these additional variables are capable of governing the network’s dynamical state and maintaining it at the edge of percolation transition, depending on parameters. The feedback acts as as a mechanism for controlling and maintaining metabolic homeostasis; this enables communication between nodes across the whole network and at the same time prevents network’s overload due to excessive propagation of the activity. The energy feedback explicitly suppresses excitation in individual neurons (locally) by disabling high frequency local spike generation. In the other words, it can be viewed as a realization of frequency dependent synaptic plasticity [39].

Despite the model qualitatively explains certain phenomena observed in neuronal cultures, the model is simplistic. It does not account for varying strengths of synaptic efficacy as well as for various plasticity mechanisms and their time scales. Accounting for these is the subject of our future work.

Ethics statement

C57Bl/6j pregnant mice (E18) were euthanized via cervical dislocation according the recommendations in the Guide for the Care and Use of Laboratory Animals of the Ministry of Health of the Russian Federation. The protocol was approved by the Committee on the Ethics of Animal Experiments of the Privolzhsky Research Medical University (Protocol Number: 09-02.2016). All efforts were made to minimize suffering.

Acknowledgments

The work supported by the Ministry of education and science of Russia (Project No. 14.Y26.31.0022).

Appendix. Proofs of technical statements

Proof of Proposition 1

Notice that $a=(1-p)^{z}\in(0,1)$ for all $p\in(0,1),z>0$. Thus $1-a^{q_{t}}\in[0,1]$ for all $q_{t}\in[0,1]$, and forward invariance of $[0,1]$ follows. The right-hand side of (2) is continuous and strictly monotone with respect to $q_{t}$ on $(0,1]$, with $q_{t}=0$ being an equilibrium. Hence all forward orbits of this map, i.e. $q_{t},q_{t+1},q_{t+2},\dots$ are monotone, and map (2) has only fixed points as attractors. Furthermore, the right-hand side of (2) is strictly concave, which implies that the number of fixed points is at most two.

If the value of $p$ is such that the right derivative of

[TABLE]

$q_{t}=0$ is less or equal to $1$ then strict concavity of $f(\cdot)$ implies that $f(q_{t})<q_{t}$ for all $q_{t}\in(0,1]$. Hence (2) has only one fixed point, $q_{t}=0$. The corresponding condition is $-z\log(1-p)\leq 1$. This fixed point is attracting: $\lim_{t\rightarrow\infty}q_{t}=0$. If $-z\log(1-p)>1$ then the trivial equilibrium $q_{t}=0$ becomes a repeller and the second fixed point $q_{t}=\tilde{q}$, $\tilde{q}\in(0,1)$ appears. At this point, the line $y=q$ and the curve $y=f(q)=1-(1-p)^{zq}$ intersect transversely. Indeed, if this is not the case then there is a point $q^{\prime}>0$ such that $1-(1-p)^{zq^{\prime}}=0$ (see Fig. 13, left panel). The latter, however, is impossible as $(1-p)\in(0,1)$. Moreover, at the point of this intersection, the slope of the curve $y=f(q)=1-(1-p)^{zq}$ is always strictly smaller than one (see Fig. 13, right panel). In order to see this, recall that the following must hold at $q=\tilde{q}$:

[TABLE]

If $df/dq(\tilde{q})\geq 1$ then strict concavity of $f$ implies that $df/dq(q)>1$ for all $q\in(0,\tilde{q}]$. This, however, contradicts to that $\tilde{q}=f(\tilde{q})$. Hence the slope of the function $f(\cdot)$ at the second fixed point, $\tilde{q}$, is strictly smaller than one, and the fixed point is locally exponentially stable. $\square$

Proof of Proposition 2

Forward invariance. Similar to the proof of Proposition 1, we observe that $0<p\sigma(E_{t},\underline{E})<1$ for all $E_{t}$. Hence $a=(1-p\sigma(E_{t},\underline{E}))^{z}\in(0,1)$ and $1-a^{q_{t}}\in[0,1]$ for all $q_{t}\in[0,1]$. This implies that $q_{t}\in[0,1]$ for all $t$ regardless of the values of $E_{t}$. Let $E_{t}\in[0,\overline{E}]$ and $E_{t}<rq_{t}$. In this case, $E_{t+1}=(1-\varepsilon)E_{t}+\varepsilon\overline{E}\in[0,\overline{E}]$. If $E_{t}\in[0,\overline{E}]$ and $E_{t}\geq rq_{t}$, then

[TABLE]

[TABLE]

Thus the domain $\{(q,E)\ |q\in[0,1],E\in[0,\overline{E}]\}$ is forward-invariant.

Statement 1. Observe, that $q^{\ast}_{1}=0$, $E^{\ast}_{1}=\overline{E}$ is always an equilibrium of (4). At this equilibrium, $\overline{E}=E^{\ast}_{1}\geq rq_{1}^{\ast}=0$, and hence the Jacobian $J(q^{\ast}_{1},E^{\ast}_{1})$ of the right-hand side of (4) at this equilibrium is:

[TABLE]

It is clear that the eigenvalues of $J(0,\overline{E})$ are $\lambda_{1}=-z\log(1-p\sigma(\overline{E},\underline{E}))$, $\lambda_{2}=1-\varepsilon$. Therefore, the fixed point $(q_{1}^{\ast},E_{1}^{\ast})$ is an attractor when $-z\log(1-p\sigma(\overline{E},\underline{E}))<1$ and is a repeller when $-z\log(1-p\sigma(\overline{E},\underline{E}))>1$. This proves statement 1 of the proposition.

Statements 2,3. Let $(q^{\ast},E^{\ast})$ with $q^{\ast}\neq 0$ be another equilibrium of (4). All such equilibria of (4) must satisfy

[TABLE]

Depending on the sign of $E^{\ast}-rq^{\ast}$, the above system splits into the following two cases:

[TABLE]

Let $g(\cdot)=\sigma(\cdot,\underline{E})$. Note that the function $g(\cdot):\mathbb{R}\rightarrow(0,1)$ is continuous and strictly increasing for $w>0$. Hence $g^{-1}(\cdot):(0,1)\rightarrow\mathbb{R}$ exists, and is continuous and strictly increasing too. Moreover, points $(q^{\ast},E^{\ast})$ satisfying (11) should satisfy conditions below (and vice-versa):

[TABLE]

Consider the functions

[TABLE]

Let $\Omega$ be the domain of the definition of the function $h(\cdot)$. If a solution of (12) exists with $q^{\ast}\in(0,1)$ then the intersection

[TABLE]

must be non-empty. The function $f(\cdot)$ is continuous on $(0,1)$, and its first derivative,

[TABLE]

is strictly negative in $(0,1)$. This implies that the function $h(\cdot)$ is strictly increasing in $(0,\overline{q})$.

Consider the case when $E^{\ast}\geq rq^{\ast}$, (left system in (12)). Given that $E^{\ast}=\overline{E}-rq^{\ast}/\varepsilon$, as a function of $q^{\ast}$, is strictly decreasing on $(0,1)$, and $h(\cdot)$ is strictly increasing on $(0,\overline{q})$, there must be at most one equilibrium of (4) in $(0,\min\{1,\overline{q}\})\subset(0,1)$. If such equilibrium exists then it must satisfy

[TABLE]

for some $q^{\ast}\in(0,1)$. This, however, is possible only if the the derivative of the right-hand side of (13) at $q^{\ast}=0$ is larger than $1$. The corresponding condition is

[TABLE]

Consider the case when $E^{\ast}<rq^{\ast}$ (the right system in (12)). Equilibria corresponding to this alternative must satisfy $q^{\ast}>\overline{E}/r$. This, however, is possible only if $\overline{E}/r<1$. The alternative condition, $\overline{E}/r\geq 1$, therefore implies that $E^{\ast}\geq rq^{\ast}$. The latter observation completes the proof of statement 2.

Let $\overline{E}/r<1$ (or $E^{\ast}<rq^{\ast}$). Given that the function $h(\cdot)$ is strictly monotone on $(0,\min\{1,\overline{q}\})$, only one equilibrium with $E^{\ast}=\overline{E}$ may exist. At this equilibrium,

[TABLE]

and

[TABLE]

The latter inequality is due to that $1-(1-p\sigma(\overline{E},\underline{E}))^{q^{\ast}z}$ is strictly concave with respect to $q^{\ast}$ on $(0,1)$ (see the proof of Proposition 1). Consider the Jacobian $J(q^{\ast},\overline{E})$:

[TABLE]

where $[*]$ stands for the corresponding entry of the Jacobian matrix. The absolute values of its eigenvalues are clearly less than $1$, and hence the fixed point is a stable attractor. $\square$.