Cross-scale excitability in networks of quadratic integrate-and-fire neurons

From the action potentials of neurons and cardiac cells to the amplification of calcium signals in oocytes, excitability is a hallmark of many biological signalling processes. In recent years, excitability in single cells has been related to multiple-timescale dynamics through canards, special solutions which determine the effective thresholds of the all-or-none responses. However, the emergence of excitability in large populations remains an open problem. Here, we show that the mechanism of excitability in large networks and mean-field descriptions of coupled quadratic integrate-and-fire (QIF) cells mirrors that of the individual components. We initially exploit the Ott-Antonsen ansatz to derive low-dimensional dynamics for the coupled network and use it to describe the structure of canards via slow periodic forcing. We demonstrate that the thresholds for onset and offset of population firing can be found in the same way as those of the single cell. We combine theoretical analysis and numerical computations to develop a novel and comprehensive framework for excitability in large populations, applicable not only to models amenable to Ott-Antonsen reduction, but also to networks without a closed-form mean-field limit, in particular sparse networks.


Introduction
Excitability is a fundamental all-or-none property of many living cells including neurons. For simplicity, we describe neuronal excitability, but this notion extends beyond membrane biophysics. It manifests itself by a very nonlinear response to a sufficiently strong external input, leading to the emission of an action potential before going back to a rest state, whereas any weaker input has no effect on the cell other than a small fluctuation of the membrane potential around its equilibrium value. The concept of excitability is well known to biologists, in particular through the existence of a non-observable boundary in the response space marking the abrupt transition from rest to spike. However, while the geometry of single-cell excitability is well understood [1], the idea of population excitability (for example in a network of coupled neurons) has been far less studied. What makes a population respond normally (such as in working memory tasks [2]) or abnormally (such as in seizures and other pathologies [3]) is a critical question in neuroscience.
The canonical class I excitable systems, such as the QIF neuron models [11], do not intrinsically possess multiple timescales. Nevertheless slow periodic forcing can bring out a bursting rhythm in theta neurons and the threshold to bursting dynamics is again formed by canard solutions [12]. Networks of QIF neurons are capable of generating similar bursting rhythms upon periodic input [13,14], so the question of network excitability in relation to threshold comes naturally in this context.
In this article, we provide a novel approach to the question of population excitability by showing that the geometry of excitability at the microscopic level scales up to large networks, involving similar key objects related to the slow-fast nature of the system. Excitability in a single neuron is determined by studying how a model responds to small perturbations in the external input or initial conditions: an excitable's system all-or-none response is mapped in parameter and phase space, by carefully selecting perturbations so that a boundary between markedly different behaviour reveals itself [4,15,16]. A major obstacle towards transferring such study from single neuron to population excitability is that it is often unclear how to translate the small perturbations used for single neurons to network and mean field level. In this paper we show that such step is possible and natural for QIF networks, the canonical class I neurons, and that canard solutions are the fabric of excitability thresholds, from single neurons to mean-field descriptions. Dissecting this transition is important, as excitability is a fundamental aspect of the cortex providing a substrate for the propagation of information from one location to another: at the single cell level (via action potentials) and at network level, via propagating waves (for example, see the recent reviews [17,18] for network models).
Initially, we build on the results by Montbrió et al. [13] (extending previous work by Ott and Antonsen [19]) in the case of dense (all-to-all) networks of QIF neurons, with randomly distributed constant inputs following a heavy-tail (Lorentzian) distribution. For this network, several groups have studied the existence of a simple mean-field limit [13,[20][21][22][23][24]. We show that excitability in large networks of this type is organised via canards in the very same way that it is at the mean-field limit, and we showcase these results computationally, by exhibiting an accurate approximation of the network threshold for networks of size N = 10 5 .
What is more, we extend this approach to a much wider class of networks of QIF neurons, including: networks with heterogeneous weights, sparsely-connected networks, networks with electrical as well as chemical synaptic coupling, networks with asymmetrically-spiking neurons, and multi-population networks that combine any of the above features.
The geometry of excitability in all such QIF neural networks beautifully persists across scales, in great generality. For large-enough networks (and up to the mean-field limit) this persistence reveals itself once we consider the correct macroscopic variables, namely the firing rate and the mean membrane potential. This is in contrast to the single-neuron level where the slow-fast variables endowed with this excitable geometry only encompass the membrane potential.
We show that systems of the type above, at any scale, support a continuous route from nonbursting to bursting solutions, upon slow (external) periodic forcing. This continuous route visits canard solutions, which form an interface for excitable transitions, from down network states (neural population silent phase) up to network bursting -as observed in [13] but without explanation of the threshold transition-as well as for the dual transitions from up network states (neural population tonic firing) down to network bursting which was not been reported before and involves the same canard geometry and dynamics.
As stated this geometry emerges through a slow periodic forcing. While the best known slow modulatory rhythms in cortex are theta oscillations with frequencies of 5 to 8 Hz, these are not the only slow rhythms that modulate the activity of the cortex. There is a global oscillation that organizes the activity of cortex operating at less than 1Hz and causes a transition from the down state (inactive) to the up state (active). This slow oscillation is present during sleep and anaesthesia (see [25] for a recent review, and also [26,27]). In addition, there exist slow "delta" oscillations whose frequency ranges from 0.5 − 4 Hz. These have been shown to modulate gamma oscillations in the whisker barrel cortex of awake mice and are coupled to respiration ( [28]). These oscillations modulate the excitability allowing gamma and other high frequency oscillations.
In the paper we follow a didactic approach, whereby calculations are performed initially on the original mean field QIF network derived by Montbrió, Pazó, and Roxin in [13], to which we add a synaptic variable as the same group considered in [21]. We henceforth denote this model as the MPR network. With this example we develop intuition and all the technical ingredients to describe the canard population thresholds. We then show how to extend this approach to more general cases.
The paper is organised as follows: in section 1 we introduce the MPR model, present excitability and routes to bursting at single cell and network level; in section 2 we present the mathematical tools to study excitability through folded-saddle canards, and we use them to interpret network excitability, which we showcase numerically in networks of 10 5 neurons; in section 3, we explain how this approach naturally extends to QIF networks in great generality; in section 4 we demonstrate that the same continuous routes to bursting exist in sparse networks, in the absence of an exact mean field limit for the network; we conclude in section 5.

QIF Network model
We study a network of N all-to-all coupled QIF neurons. The ith neuron has membrane potential V i , synaptic variable s i , and is subject to both a background current η i , and an external, zero-mean current I(t) = A sin(εt), leading to for 1 ⩽ i ⩽ N. We refer to the sum K i = η i + I(t), as the (external) input to the ith cell. The ODEs above hold between two consecutive firing times: these are the finite, computable times at which a membrane potential in the network diverges to +1. Each time this condition is met by V i we: (i) stop the simulation, (ii) reset V i to −1, (iii) send a spike to all synapses, whose values are instantaneously incremented by an amount 1/N; (iv) restart the simulation of system (1) from these updated initial conditions. A voltage threshold at +1 is clearly non physical, but one typically addresses this problem in two ways. Firstly, it is possible to use the transformation V i = tan θ i /2 to obtain one then computes a firing event when θ i crosses the value π/2 from below. Secondly, voltages at V t = 1, V r = −1 may be replaced by finite, large values V t = −V r . Below we shall present two types of numerical simulations: for N = 1, we use the above transform, hence thresholds are attained at V ! 1; for N > 1, we use the second strategy instead. The simulation for N = 1 in θ presents no appreciable difference to the one with N = 1 with finite threshold and reset. This system is ideally suited to study excitability across scales because: (i) we can analyse and compare single cell-dynamics, N = 1, network dynamics N � 1, and mean-field dynamics N ! 1; (ii) one can switch from constant input currents (ε = 0) to slowly-varying, oscillatory currents (0 < ε � 1) in order to uncover transitions between various cellular regimes.

Single-neuron excitability
Let us set N = 1, ε = 0, and examine a QIF neuron with self-coupled synapse [29], subject to a constant input current K 1 = η 1 , as in Fig 1a and 1e. When 0 < τ s � 1 and Jτ s is sufficiently large, the cell supports two coexisting attracting states: an equilibrium (down state), and a periodic solution with tonic firing (up state), separated by an intermediate unstable equilibrium. The equilibria belong to a curve which folds when the input K 1 is null; periodic solutions collide with unstable equilibria at a homoclinic bifurcation, when K 1 = h. In passing we note that in Fig 1a, 1b, 1e and 1f are only schematics of bifurcation diagrams, without units or scales, which is why they feature "states" on the ordinates. This is in contrast to Fig 1c, 1d, 1g and 1h, which show voltages of computed trajectories.
When ε = 0 and η 1 2 (h, 0) in the intrinsically bistable regime, Fig 1a, initial conditions determine whether the voltage is attracted to the down or to the up state; the threshold is given by the middle unstable state. When η 1 > 0, the only attractor is the periodic solution, hence the cell is intrinsically tonic, Fig 1e. We are interested in how the cell (and later on, the network) transitions from the rest state to the repetitive firing state, and how the two states are concatenated together to form a bursting state. In a standard QIF model without synapses (J = 0), there is no bistability (the up state is to the right of the fold): this changes some of the waveforms supported by the cell, but not the mechanisms we aim to describe, namely the bursting transitions between down and up states [12].
To study these transitions, one sets ε > 0 small and hence examines slow forcing [1,4], as sketched in Fig 1b and 1f. This causes the input to oscillate around the mean value η 1 (see ellipses on the horizontal axes). The exact dynamics of the system depends on A, the amplitude of the input oscillations, and on the sign of η 1 . By varying these two parameters, one can construct a great variety of solutions where up and down states alternate. Some trajectories stand out, in that they signal the onset or termination of a phase. With reference to Fig 1b, smallamplitude forcing with average η 1 2 (h, 0) causes the cell to oscillate around its rest state; these are subthreshold oscillations, which stick to the down state at all times; upon increasing the amplitude A (see ellipse around η 1 ), the trajectories reach a turning point, when η 1 + A � 0, that is, for A � −η 1 ; near this value, there are trajectories which follow the branch of unstable states for increasingly longer times, before jumping to the down state, or to the up state (segments 1-2). A temporal profile of solution jumping down is given in Fig 1c, obtained for A = 0.20318, ε = 0.01 (green); a profile of a solution jumping up is also in Fig 1c, for A to 0.20319 (purple curve). The narrow region of parameter space between A = 0.20318 and A = 0.20319 contains an entire family of solutions, which spend O(1) times near the repelling branch of equilibria of the system with ε = 0. This is surprising: when ε is small, one would expect the ε = 0 analysis to play a role, and therefore trajectories to be repelled exponentially fast from the unstable branch; on the contrary, here we consider trajectories that stay close to the unstable branch for long times. Orbits with such a feature, like the ones labelled with segments 1-2 (and 3-4) in Fig 1 are called canards. These counter-intuitive canard solutions (marked with colors and numbers in Fig 1b), constitute a computable interface between subthreshold oscillations (down-down orbits) and Sketch (not to scale) of the bifurcation diagram of steady states (curve) and periodic solutions (cylinder) of the single QIF neuron subject to a constant input K 1 = η 1 (ε = 0). A stable quiescent state (down state) coexists with a stable tonic firing solution (up state), separated by an unstable equilibrium (dashed curve). A homoclinic bifurcation is present when K 1 = h. In this intrinsically bistable regime (K 1 = η 1 2 (h, 0)), the cell selects the up or down state depending on initial conditions. (b): When 0 < ε � 1, K 1 (t) = η 1 + A sin(εt) becomes a slowly varying quantity, oscillating around the value of η 1 (ellipses on the K 1 axes) with amplitude A, and transitions between the up and the down phases become possible. The onset between phases is determined by a family of canard solutions 1-2 (see text); in the bistable regime they appear in the down-down (green), and down-up (purple) transitions. In the latter the variables (V 1 , K 1 ) are used, and are superimposed on the curve of equilibria of the ε = 0 system (grey parabola), providing evidence of canard behaviour (1)(2).
The fact that canards occur in such a small range of parameter space makes one wonder whether they are detectable and useful in nature. As we shall see below, canard solutions determine the effective thresholds of the all-or-none responses in a model; they act as excitability thresholds, whose biological relevance is well established [4,15,30], even though isolating such threshold may be challenging in experiments. To observe a canard in experiments, one would need access to phase-plane information similar to Fig 1b and 1f, but most electrophysiological experiments record only one observable time trace, like membrane potential, which does not allow for a phase plane interpretation. Nevertheless, time traces can display canard signatures, as we show in Fig 1c and 1g, and similar traits are found in experiments [31], in particular when discussing delayed onset of firing [32]. In addition, analog circuits engineered to reproduce excitability display a clear canard behaviour as they provide access to experimental phase plane data [33].
We then move to the intrinsically tonic regime, for η 1 > 0, Fig 1e-1h. We find up-up and up-down orbits with canards in Fig 1f-1h, whereas it is possible to prove that down-down and down-up canards cannot exist, that is, the transition at the fold is a jump. This is why we have segments of type 3-4, but not of type 1-2, in this scenario. In Fig 1g and 1h show many spikes, most of which are greyed-out for enhancing visibility of the canard segments. In passing we note that segments of type 3-4 are also present in the bistable scenario, but we do not discuss them in the single cell, for the sake of brevity.
The orbits described above capture excitability transitions at a single-cell level. Developing a mathematical understanding of these special orbits is crucial, because canards act as basin boundaries between different cellular responses: biophysical and idealised single-cell models support generically continuous canard-mediated transitions [7,9], as we will exemplify in a moment. The transitions are brutally sharp, but can be continuous, even though they may appear discontinuous upon running simulations. The main contribution of this paper is to show that this scenario also occurs generically and robustly in networks of type-I neurons, of which QIF are universal prototypes.
In addition, the canard-mediated transition from non-bursting to bursting states in networks of QIF neurons transfers across scales: in an isolated QIF cell, as well as in a networks of QIF cells, there exists a continuous route from non-bursting to bursting states, and this path is made of solutions with canard segments such as the ones seen in   We label the scenario above as a up-to-down route to bursting. Other routes to bursting are also possible (from down to up states, for instance) when η varies. We do not pursue a classification of bursting routes for a single QIF cell, but we will do so for the mean-field network, after showing that continuous routes to bursting persist across scales, when N ! 1.

Network excitability
Let us now consider the network (1), together with reset conditions, subject to random background currents: η i are taken from the Lorentzian distribution with density gðZÞ ¼ D=ðpðZ À � ZÞ 2 þ pD 2 Þ, hence the network is heterogeneous, with some neurons in the bistable regime, and others in the tonic regime. However, the centre of the distribution � Z will turn out to play a role. If � Z < 0 (� Z > 0) we say that distribution, or the network, is predominantly bistable (tonic). For N ! 1, there is a well-known mean-field limit [13,21] for the coupled system: where r, v, s are the mean firing rate, mean membrane potential, and mean synaptic input, respectively. The external input is now given by KðtÞ ¼ � Z þ IðtÞ. Recall that at microscopic level the background current is constant but heterogeneous (η i , sampled from a Lorentzian distribution), whereas in the mean field limit it is constant and homogeneous (equal to � Z). Note that alternative Ott-Antonsen QIF reductions of QIF networks use amplitude and phase of a complex-valued order parameter in place of mean voltage and rate [20,23]. The order-parameter and rate-voltage descriptions are related through a conformal mapping [13].
It is important to remark that the mean-field limit introduces a new quantity, the population firing rate, r [21] that is not a part of the single or finite system of equations. It emerges in the limit as N ! 1 of the microscopic model (1). As seen in Fig 3, when ε = 0, hence KðtÞ � � Z, the equilibria of the system lie on an S-shaped curve. More precisely, it can be shown that the curve has no folds or 2 folds ( [34], Eq 12 and Fig 2(c)). Henceforth we shall assume that J is sufficiently large to guarantee the existence of 2 folds, which must occur for negative values of K.
In the down state all neurons are close to rest (quiescent network state), whereas the up state corresponds to asynchronous network tonic firing, an averaged version of the tonic state in Fig 1e-1f. The up state can be a stable focus away from the fold and have complex eigenvalues. Between these two stable fixed points is an unstable (saddle) point that serves as a separatrix between the two stable states. Remarkably, when 0 < ε � 1, the geometry of excitability still persists in this macroscopic description, and transitions are now determined by the distribution peak at � Z. We now show that the orbits of the mean field and of the network directly parallel the transitions of the single neuron model, involving the same canard types. The qualitative difference between the two behaviors is more striking in the time traces of the two curves, shown in Fig 4a, where the transition to the up state is accompanied by a burst while for the smaller input, there is just a subthreshold oscillation. The reason for the burst is that the up state is a stable spiral for a range of input values. When A � À � Z, as expected, the network jumps (without canards) from asynchronous to synchronous firing, following the upper branch of the S-shaped curve in Figs 3a and 4b. In conclusion, the trajectories displayed in Fig 4a and 4b are the mean-field equivalents of the single cell ones in the bistable regime, shown in Fig 1c and 1d, respectively. They undergo similar transitions from down-down to down-up states, with canards being the threshold.
The tonic case (� Z ¼ 5) is simulated in Fig 4c and 4d. A similar situation occurs in this case, but the dynamics revolves around the upper fold F + , and features canards in up-up and updown orbits.

Folded-saddle canard behavior across scales
We now derive three central results of the paper: firstly, we characterise the mean-field transitions described above, valid at the ODE level (system (2)), using standard methods from Geometric Singular Perturbation Theory [35]; secondly, we use this characterisation to infer the existence of canard behavior at the network level (system (1)) with N = 10 5 neurons, and provide numerical evidence of this phenomenon; thirdly, we explore canard-mediated routes to bursting at the network level, a feature that persists from the single neuron level. The latter results are remarkable and novel, as canard behavior in large networks is greatly unexplored, in particular for systems with resets and random data.
As anticipated, we initially present the theory for the network described above, and then adapt these results to more general networks. To study the behavior of (2) for small ε > 0, we extend the system with two ODEs describing the oscillatory dynamics of K(t), namely K 0 ¼ εQ; Note that in [13,24], the mean-field limit model assumes instantaneous synaptic processing, which amounts to taking τ s = 0 and replacing s by r in the equations for r and v. However, assuming 0 < τ s � 1, that is, a fast synapse, does not change anything about the threshold analysis below while making it more general. Now that we have expressed the slow dynamics of the current I using a second-order harmonic equation, we rescale time in Eqs (2) and (3) so as to parametrise them by the slow time τ = t/ε, as was done in [36], and obtain To shed further light onto the transition from low-rate (down) states to high-rate (up) states of the mean-field, it is key to consider the slow limit of (4) as the forcing speed ε tends to 0. Therefore, we set ε = 0 in (4) and obtain the three algebraic constraints In the ε ! 0 limit, the system variables (r, v, s, K, Q) evolve on a three-dimensional manifold in R 5 , the so-called critical manifold, given by The subscript 0 in S 0 refers to the fact that this manifold is found by setting ε = 0 in (4). The transitions discussed in this paper occur when S 0 is folded, and it is the union of two attracting and one repelling submanifolds. These conditions occur generically: for common choices of parameters J, Δ and � Z (see [13], Fig 1(a)), one finds that S 0 has two loci of folds (two lines of folds), F + and F − , corresponding to the set {Dψ(v) ≔ ψ 0 (v) = 0}. A projection of the manifold S 0 onto the (K, v) plane is visible in Fig 4b and 4d, where the fold lines project onto points F ± ; compare with Fig 5 where S 0 is projected onto the (K, Q, v) space and the fold lines are fully visible.
The ε ! 0 limit introduced above corresponds to a differential-algebraic problem referred to as the slow subsystem of the original equation, and in the present case it reduces to The algebraic constraint in (6) hides the dynamics of v in this slow limit. To reveal it, we then differentiate the constraint with respect to time, and obtain the following set of ODEs One can relate system (7) to the canards shown in Fig 5: we are considering the ε = 0 dynamics, hence we focus on the black curves on S 0 , and we take the one in Fig 5a as an example. It would appear that system (7) breaks down at points where ψ 0 (v) = 0, that is, along the fold set F ± . Along such folds the first equation reduces to Q = 0 hence (7) is undefined at points along the folds where Q 6 ¼ 0. However, an inspection of Fig 5a shows that the flow is well defined at one specific point, which is called a folded singularity and is marked as (fs). In fact, there are two other points on F + where the trajectory seems to cross the fold. However, these are not associated with canard dynamics, and the flow is not passing through the fold, as the system (7) is singular at those points.
It is around this point that canard solutions are born, because the trajectory passes through (fs) from an attracting to a repelling sheet of S 0 . We also note that for ε � 1 this behaviour We superimpose them onto the (grey) surface of ε = 0 mean-field equilibria, S 0 , and fold lines F ± (also shown in Fig 4b and 4d). On the surface are visible the folded saddle singularity (fs) and its associated canards (black). The networks follows the canard orbits predicted by the mean-field theory remarkably well. persists (green curve in Fig 5a): slow-fast theory predicts [37] that canards at ε = 0 survive for small-enough ε > 0 and, as we will see below, they organise the excitable structure of QIF networks. We therefore investigate the passage through (fs) more precisely: intuitively, at (fs) ψ 0 (v) = 0 and Q = 0, so that the quotient Q/ψ 0 (v) stays finite, and the flow of (7) is well defined. We formalise this step by: (i) desingularising (7) with a time rescaling, (ii) identifying (fs) as an equilibrium point of the desingularised problem, (iii) classifying the type of equilibrium, which in turn determines the type of canard in the original system (as done in [36]).
In step (i), we desingularise system (7), and rescale time by a factor −ψ 0 (v), which eliminates the prefactor to _ v, and regularises the problem, leading to the desingularised reduced system (DRS). v 0 ¼ Q; An important subtlety is that the time rescaling by −ψ 0 (v(t)) depends on the state variables v, hence the time orientation depends on the position on S 0 . In fact, the time rescaling transforming the reduced system (7) into the desingularised reduced system (8) is such that the flow in both systems has the same orientation on the attracting sheets of S 0 but opposite orientation on its repelling sheet.
In step (iii) we study the linear stability of these equilibria, which is determined by the Jacobian matrix ; whose eigenvalues are given by l ¼ � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À c 00 ðv � Þcðv � Þ q : The equilibrium (v�, 0) is therefore either a saddle or a center: in the former case the point (fs) is called a folded saddle, and gives rise to folded-saddle canards; in the latter case (fs) is a folded centre; it is known that this singularity does not give rise to canards, but rather to a discontinuous transition. A quick calculation reveals that: (1) In the bistable regime, when � Z < 0, there are two (fs) points, one on F + and one on F − , both of which are folded saddles; both folded saddles are visible (in projection) as F + , F − in Fig 4b, whereas only the latter is shown in Fig 5a and 5b. (2) In the tonic regime, when � Z > 0, the singularity on F + is still a folded saddle, whereas the one on F − is a folded centre, and both are visible in projection in Fig 4d; we show only the folded saddle as (fs) in Fig 5c and 5d. In passing, we note that these findings remain valid for a much larger class of networks and mean field limits, as we will show below. Fig 5 displays the dynamics of system (7) (ε = 0, in black) superimposed on the dynamics of (1) with N = 10 5 (ε > 0, in color), for both subthreshold and suprathreshold forcing . Fig 5a  shows, in green, the full network orbit with N = 10 5 for A slightly below threshold. A canard segment is visible, where the trajectory hugs the black curve on the fold before falling back to the down state. Fig 5b shows the same projection, but for slightly larger A; in this case, the trajectory makes the jump to the up state before falling back down. The folded-saddle canards predicted by mean-field theory (in black) are in striking agreement with the behaviour of the network, and the correspondence with the similarly coloured mean-field transitions in Figs 4a, 4b, 1c and 1d is remarkable. These solutions define analogue scenarios to the down-down and down-up states for the single neuron model, but they are exhibited at the level of the network. Interestingly, at single cell level we have a bistable neuron (η 1 < 0), while at the network level some neurons will be tonic, but the network is predominantly bistable (� Z < 0). Similar considerations apply to the predominantly tonic network scenario, showcased in Fig 5c and 5d, and corresponding to Figs 4c, 4d, 1g and 1h.

Continuous routes to bursting
We will now show that for any value of � Z, the network robustly supports a continuous route from a non-bursting state to a bursting state upon increasing the input amplitude; this transition involves a canard explosion, as seen in Fig 2 for the single cell.
For sufficiently large coupling values J, the critical manifold S 0 is S-shaped, with folds occurring at η + < η − < 0. As depicted in Fig 6a, there exist 4 different scenarios depending solely on Z in each region are � Z ¼ À 6:5, � Z ¼ À 5, � Z ¼ À 3:5 and � Z ¼ À 2, respectively. Panels (b-d) display both the solution branches obtained by varying the forcing amplitude A for a given value of � Z (top), and a selection of solutions on the branch, plotted in the phase plane (K, v) on top of S 0 (bottom). As observed, a continuous branch of solutions bridging from the non-bursting regime to the bursting regime always exists, regardless of the � Z value. This branch connects in parameter space a down-down non-bursting solution to a down-up bursting one for � Z < Z 0 (panels b,c), or an up-up non-bursting solution to an up-down bursting one for � Z > Z 0 (panels d,e). Additionally, for Z þ < � Z < Z À , another solution branch exists, starting at low A amplitude and which does not connect to the bursting regime; this branch contains upup (resp. down-down) solutions for Z þ < � Z < Z 0 (resp. Z 0 < � Z < Z À ). Parameter values are: Δ = 1, J = 15, τ s = 0.02, ε = 0.05; � Z and A as indicated in the panels.
https://doi.org/10.1371/journal.pcbi.1010569.g006 the value of � Z with respect to the fold values η ± and the midpoint η 0 = (η + + η − )/2 between them. This is due to the symmetry of the forcing amplitude around � Z. Case I: � Z < Z þ . Since the forcing oscillates around � Z < Z þ , when a low-amplitude forcing is switched on, the network can only oscillate near the bottom branch of S 0 (orbit 1 in Fig 6b). Upon increasing the amplitude of the forcing, the orbit reaches the fold and hugs the middle unstable branch of S 0 (orbit 2 in Fig 6b). A continuum of orbits with canard segments are now visited by the network: in this transition to bursting visible in the (A, Δr) diagram in Fig 6b, the branch of solutions is almost vertical, as A varies in a tiny region of parameter space. As we ascend the branch, we pass from a down-to-down solution (green, labelled as 2) to a down-toup solution (purple, labelled 3) with canard segments. Past the vertical branch, we obtain a fully developed bursting solution (4). The solution branch in Fig 6b reveals a continuous down-to-up route from non bursting to bursting network solutions. All solutions with canard segments in this case are of down-to-up or down-to-down type, as witnessed by the green and purple coloring.
Case II: Z þ < � Z < Z 0 . Since we enter the bistable region of S 0 , we can switch on the forcing near two starting points, on the lower and upper branch, respectively. If one starts from the lower branch, the same considerations as in Case I are valid, and we have a continuous downto-up route (see green branch in Fig 6c). Starting from the upper branch, low amplitude forcing generate state hovering near the upper branch (solution a in Fig 6c)). When the forcing amplitude increases the orbit grows a canard segment from the upper fold F + (solution b). However, due to the proximity between � Z and F + , the branch can only grow canards up to the solution labelled c, and folds back onto itself while displaying solutions d-f. Therefore in Case II there is a continuous down-to-up transition (and no up-to-down) transition to bursting.
Case III: Z 0 < � Z < Z À . This scenario is the mirror image of Case II. The network possesses continuous down-to-up transition to bursting, but no down-to-up transition, which is interrupted (Fig 6c). The continuous transition is explained in case IV below.
Case IV: � Z > Z À . When the forcing is small, the solution can only stay near the upper branch of S 0 (solution 1 in Fig 6e). We can still transition continuously from this non-bursting solution to a bursting solution (solution 4) via canards that start near F + . This case mirrors Case I, but involves canards of up-to-up and up-to-down type.

Networks with heterogeneous currents
Let us now consider generalisations of the MPR network (1) and mean-field limit (2), for which the excitability scenario described above still holds. We will discuss the generalisations only at the level of the mean field, and refer to existing literature for descriptions of the corresponding microscopic networks.
The starting point is the following generalisation of the MPR mean field (2) where: I(t) = A sin(εt) is a slow zero-mean periodic forcing as before; Δ, J, � Z and τ s are parameters as before; Γ, g, and a are additional parameters. This generalisation encompasses a variety of exact mean-field limits of QIF networks including: 1. The MPR network with heterogeneous background currents η i sampled using a Cauchy distribution with peak at � Z and half-width at half-maximum (HWHM) Δ. To recover this model from (9) [13,22] one sets Γ = 0, g = 0, τ s = 0.
3. The MPR network with first-order fast or slow synapses [14,21] characteristic time τ s , which corresponds to Γ = 0, g = 0, τ s 6 ¼ 0. Note that the analysis done above on (1) assumed 0 < τ s � 1 (fast synapse), however it is still valid for τ s = O(1) (slow synapse), because the folded-saddle structure only requires the existence of a slow periodic forcing (3). This dynamics also persists with second-order synapses.
5. The modified QIF network studied in [38] with electrical coupling and asymmetric spikes, which differs from the previous case only by a 6 ¼ 1.
Before further extending the networks analysable with the proposed formalism, let us rewrite (9) using the generalised coefficientsGðG; gÞ ¼ G=p À g,JðJ; g; aÞ ¼ J þ g ln a, yielding It is apparent that one can analyse the slow-fast structure of this system in the exact same way as we have done for system (4), with a cubic-shaped critical manifold given by where ψ is as in (5) withJ instead of J. It is apparent that the same folded-saddle dynamics organise the excitable structure of the corresponding mean-field model and can be observed in associated large-enough generalised QIF networks. Instead of pursuing this analysis, we introduce first a further generalisation, namely we consider the exact mean-field limit of p synaptically coupled populations of QIF networks, where we suppose, for simplicity, that only one population (the k-th one) receives a slow external periodic forcing. The coupled equations read for i = 1, � � �, p, where δ is the Kronecker symbol, and ρ i , ν i , σ i are the functions defined in system (10) for population-specific choices of parametersG i , Δ i , � Z i , and (τ s ) i . The critical manifold of system (11) is defined by the algebraic constraints where Hence, the critical manifold S 0 can be compactly written as The expression for S 0 contains p independent algebraic conditions in R pþ2 , therefore the critical manifold is indeed a surface, which is consistent with the fact that system (11) has two slow variables. As a consequence, one can write the reduced system associated with (11) in the form The system above mirrors the constrained system (6) in the single-population MPR network. Proceeding like in the single-population case, we differentiate the algebraic constraints with respect to time and project the resulting limiting system onto the (v k , Q)-plane to obtain where v 1 , � � �, v p satisfy the algebraic constraints defining the critical manifold, that is, the first p equations in (14). We recognise in system (15) the same form as (7) for the one-population mean-field limit. We note that the starting system has p + 2 equations, but the reduced system (15) has only 2 equations, and it is singular at the fold set of the kth population system, given by the condition @ v k C k ðv 1 ; � � � ; v k� ; � � � ; v p Þ ¼ 0. The 2-dimensional system (15) is singular along that fold, and it can be desingularised as in the one-population problem. The folded-saddle and folded-centre classification carries through in this case. Hence we can conclude that the same canard-induced excitability scenario appears in the generic p-population case described above.

Canard transitions across scales in sparse networks
The slow-fast scenarios uncovered in the previous section are valid in a large variety of (all-toall coupled) QIF networks with exact mean-field limits. We now present evidence that the phenomenon persists in sparse networks, for which no exact mean field limit has been derived to date.
We present this extension for two main reasons: on one hand, we show that the mechanism discussed in the previous section extends further, to sparse networks; on the other, we want to emphasise that the availability of a mean-field description is not strictly necessary for the canard phenomenon, which is supported by generic network systems of QIFs with finite size.
In Section (1.3) we studied networks with exact mean fields: since an ODE description was available for the case N ! 1, we used these ODEs to predict the region in parameter space where canard dynamics occur, and to classify the folded singularities organising the transition from non bursting to bursting patterns. However, networks with finite size also support canard-mediated transitions, as evidenced in Figs 2 and 4, where the orbits are computed for a very small (N = 1) and a very large (N = 10 5 ) network, respectively.
Having a mean field description at our disposal is useful to pinpoint regions of parameter space where canards will occur, through the study of S 0 and its folded lines; also, large networks of neurons possess canard solutions that are almost indistinguishable from their mean field ones. However, canard mediated transitions are present (and can be documented) in finite-size networks, even when the mean field is inexact, or unavailable in closed form.
To substantiate this claim, we study a sparse network of N synaptically-coupled QIF neurons. For a similar network, a heuristic mean field description has been proposed, based on sparsity scaling arguments [22,39]. The heuristic mean field is in the form (10) and hence canards of folded saddle type are supported by this set of ODEs. However, this mean field is not the exact limit of a network of QIF neurons, and the extent to which the mean field approximates the finite-size network is also immaterial to find canards: as we shall see, both the heuristic mean field and the finite-size network have canard-mediated routes to bursting, even if the two models do not agree well in certain regions of parameter space.
We consider N synaptically-coupled QIF neurons of the following form arrive at the following candidate approximate mean field given by Differently from [22,39], the model above has heterogeneous currents, in addition to sparse, heterogeneous connectivity. Also we consider an excitatory neuronal population, as opposed to a inhibitory one. The inhibitory population considered in [22,39], with J = 1 and a term −Js in the v-equation, is not suitable for studying excitability and transition to bursting, as the critical manifold is not folded. In passing we note that inhibition prevents canard behaviour in this particular model, but not in others: the generalisations of the MPR network given in 3 is valid for coupled populations of excitatory and inhibitory networks, which generically support canard behaviour.
A bursting orbit for a network with N = 10 4 neurons is visible in Fig 7a, in the (v, K, r)space. This simulation is done for a network with sharply peaked current distribution (Δ � 1), which generates a particular bursting pattern, as we will now discuss.
Assuming that the heuristic mean-field description approximates the network simulation, one can reason as in a standard MPR network: the burst is due to a family of foci on the upper branch of an cubic-like critical manifold (visible in grey in the figure). The figure shows that the manifold S 0 of the candidate mean field captures well the geometry of the bursting orbit, and in particular it displays a canard segment along the repelling branch of S 0 .
A further inspection of S 0 reveals that, because Δ is small, both folds F ± of S 0 occur at vanishingly small values of r. We observe down-to-down and down-to-up orbits in the network as well as in the heuristic mean field (see Fig 7b), and this suggests the possibility of a continuous down-to-up route to bursting.
It is important to note that the presence of nearby down-down and down-up solutions, on its own, suffices to get a hint of the presence of canards; in this case, we also have an approximating heuristic mean-field description with a computable manifold S 0 , which clearly helps us finding values of parameters where the nearby down-down and down-up solution exist; the considerations that follow, however, hold for the finite-size network, even if we obliterate S 0 from Fig 7a and 7b, and from the discussion above.
To uncover a continuous route to bursting in the network, we compute several orbits of the system when A varies between A ; while A changes, we keep the connectivity matrix W constant, that is, we extract it once and reuse it thereafter. The results are given in Fig 7c-7f. In Fig 7d  we show the voltage profiles. As anticipated, the canard structure of this network is peculiar, because the currents are almost homogeneously distributed: the canard segment in these orbits stretches along v = 0; during this transition, however, the rate increases sharply, as seen in the raster plot Fig 7e and in the histograms in Fig 7f. This means that the excitability threshold for this network occurs for states at constant voltage and progressively large rate, unlike in the cases presented before. This is induced by the fact that the critical manifold S 0 is different between system (17) and system (4). The former is a perturbation of the latter by a term proportional to s, which gives rise to a very sharp fold followed by a middle branch along which v remains almost constant. This effect is visible on Fig 7a. In Fig 7c we compare the routes to bursting for the network and the candidate mean field. We find a good agreement in the subthreshold regime, as well as the presence of canard transitions in both systems, marked by quasi-vertical branch segments.
We also notice two types of discrepancies. Firstly, the value of A at the quasi-vertical segment differs slightly in the N = 10 4 and in the mean field; this is to be expected, given that the mean field is only heuristic. Solution types along the canard explosion, however, are very similar in the two systems. Secondly, we notice a discrepancy in the maximum voltage v, which is the solution measure chosen for the network, especially in the bursting regime. This type of discrepancy is not unexpected: in the finite-size network, each neuron has a voltage at most equal to the reset value, hence the maximum mean voltage is capped in the network; the heuristic mean field, on the other hand, supports solutions with a very large maximum voltage.
As anticipated, the discrepancy in Fig 7c is relevant if one wants to assess the accuracy of the mean field in the bursting regime, but is immaterial for the canard-mediated route to bursting: the blue route in Fig 7c, and the data in Fig 7d-7f show such route independently of the existence of S 0 , or of the accuracy of the heuristic mean field.

Conclusions
The geometry of excitability and transition to bursting behavior in single neurons and allied systems is governed by canard solutions, which act as thresholds and determine the response of the system to slow parametric changes. We have shown that this structure carries over in the mean-field limit of large populations of excitable cells, as well as in large finite systems.
In these cases, the average voltage of the population plays the same role as the voltage in the single cell, and a well defined rate emerges as a new macroscopic variable. If a separation of time scales exists between external input and voltage at the level of a single cell, such separation persists at network level, between the input and the coupled mean voltage and rate.
Two main results have been discovered at network level: (i) a large class of networks of QIF neurons subject to an external stimulus with Cauchy-distributed background currents undergoes a continuous transition to bursting (either up-to-down, or down-to-up) upon increasing the forcing amplitude; to the best of our knowledge, this statement holds for any QIF network amenable to the Ott-Antonsen reduction currently derived in the literature, that is, expressible as system (9). (ii) The canard-mediated route to bursting is present also in sparse networks, for which there is no exact limit; obtaining an approximate mean-field limit is convenient to compare network trajectories to low-dimensional manifolds, but it is not strictly necessary, and in fact continuous routes to bursting are present in small networks too.
Since QIF neurons are general representatives of type-I neurons, we expect that similar properties will survive in networks of more realistic cells, up to their mean field-limit, which need not be an ODE. Introduction of inhibitory networks could provide a connection between this work and the concept of balanced networks where, depending on the details of the connection strengths, the firing is either driven by the mean input (analogous to our tonic behavior) or by the fluctuations (analogous to the excitable case). Another open interesting research direction is to investigate heterogeneous networks whose distributions are asymmetric, thereby violating the assumption required for the Ott-Antonsen exact mean-field derivation: without a mean field, we could still be able to observe a contraction of the dynamics to a lowdimensional manifold, and a canard-mediated excitability threshold. There is currently no available theory for such cases, except for continuous coarse-grained networks [36,40], but a numerical exploration of the averaged voltages and synaptic variables may reveal an underlying low-dimensional structure, similarly to what has been found in spatially-extended neural field models.