Biological emergent properties in non-spiking neural networks

: A central goal of neuroscience is to understand the way nervous systems work to produce behavior. Experimental measurements in freely moving animals ( e.g. in the C. elegans worm) suggest that ON- and OFF-states in non-spiking nervous tissues underlie many physiological behaviors. Such states are deﬁned by the collective activity of non-spiking neurons with correlated up- and down-states of their membrane potentials. How these network states emerge from the intrinsic neuron dynamics and their couplings remains unclear. In this paper, we develop a rigorous mathematical framework for better understanding their emergence. To that end, we use a recent simple phenomenological model capable of reproducing the experimental behavior of non-spiking neurons. The analysis of the stationary points and the bifurcation dynamics of this model are performed. Then, we give mathematical conditions to monitor the impact of network activity on intrinsic neuron properties. From then on, we highlight that ON- and OFF-states in non-spiking coupled neurons could be a consequence of bistable synaptic inputs, and not of intrinsic neuron dynamics. In other words, the apparent up- and down-states in the neuron’s bimodal voltage distribution do not necessarily result from an intrinsic bistability of the cell. Rather, these states could be driven by bistable presynaptic neurons, ubiquitous in non-spiking nervous tissues, which dictate their behaviors to their postsynaptic ones.


Introduction
Spiking neurons and spiking neural networks are the most widely studied cases in computational neuroscience due to their ubiquity in neural processes [1][2][3][4][5][6]. However, such neurons are not the only information processing component of the nervous system. Indeed, numerous nervous tissues in both vertebrate and invertebrate animals function normally without triggering spikes. Spiking neurons driven by bistable neurons which dictate their behaviors to their postsynaptic ones through bistable synaptic inputs. The generality of our study gives a rigorous framework for better understanding how the behavior of coupled non-spiking neurons depends on their intrinsic dynamics and their couplings. In other words, it allows to gain insight into how intrinsic dynamics and couplings of single non-spiking neurons impact neural network behaviors.

Amplitude Duration Waveform
Graded potential Action potential and bistable (AFD) non-spiking neurons. The evolution of the membrane potential has been recorded in response to a series of current injections, in the space of 5 seconds, starting from −15 pA and increasing to 35 pA by 5 pA increments. In all cases, the higher the current intensity, the higher the amplitude of the membrane potential. The experimental data have been reproduced from [33] and [17] with the consent of the authors.
The remainder of this paper is organized as follows. Section 2 presents the new non-spiking phenomenological model recently developed in [32]. The analysis of its stationary points and its bifurcation dynamics is carried out in Section 3. Section 4 presents elements for modeling synaptic connections in non-spiking neural networks, with a particular focus on C. elegans. Then, we study both mathematically and numerically how the intrinsic dynamics of a simple presynaptic neuron impact the behavior of a postsynaptic neuron.

Presentation of the phenomenological non-spiking model
Non-spiking neurons are ubiquitous in many nervous tissues. For this reason, a recent simple model has been proposed in [32] as an alternative to conductance-based models (a.k.a. Hodgkin-Huxley type models), known to be computationally expensive and difficult to treat mathematically. In this section, we briefly present this simple non-spiking neuron model, derived from the bifurcations' study of conductance-based models (CBMs) of non-spiking neurons [17,34].
Let V represent the membrane potential of a neuron. The neuron 1-D model takes the general form where I describes an injection current, τ the constant time for which V reaches its equilibrium value V * , and V → f (V) a cubic function which reads as The parameters a, b, c, and d are dimensionless. These are estimated to reproduce the steady-state current of the neuron, noted I ∞ , shown to determine the neuro-computational features of non-spiking neurons [17]. The biological meaning of the steady-state current is the following: in typical voltageclamp experiments, the membrane potential is clamped at several values V H (H for hold) for which the resulting currents are measured ( Figure 2.A, left). Asymptotic values (t → ∞) of those currents, depending only on V H , are called steady-state currents ( Figure 2.A, right). As a representative example, the experimental steady-state currents of the RIM, AIY and AFD neurons are represented in Figure 9. The stationary points V * of the model (2.1) satisfy f (V * ) = I so that the function V → f (V), in particular its shape, determines the neuro-computational features of the non-spiking model (2.1). Indeed, two typical shapes of the curve f exist, monotonic and N-shaped ( Figure 2.B), each involving two fundamental different behaviors of the neuron: • As shown in Figure 2.B (left), a model with a monotonic function f has only one equilibrium for any value of I. Such a model displays a near-linear behavior characterized by smooth depolarizations or hyperpolarizations from the resting potential, such as the RIM and AIY neurons ( Figure 1.B, left). • As shown in Figure 2.B (right), a model with a N-shaped function f has one, two, or three equilibria depending on the value of I. Such a model displays a bistable behavior characterized by voltage jumps between the resting potential and a depolarized potential, such as the AFD neuron ( Figure 1.B, right).
In a previous work [32], the parameters of the model (2.1) were estimated using the differential evolution algorithm [35] to fit the experimental steady-state current of the RIM, AIY and AFD neurons V c 1 * and V c 2 * correspond to equilibria for a current injection I = c 1 and I = c 2 respectively. (Right) N-shaped steady-state current. The number of equilibria of the system depends on the value of I. For the sake of readibility, we highlight equilibria only for I = c 1 , noted V c 1 1 * , V c 1 2 * and V c 1 3 * . displayed in Figure 9. Further, we reproduced the mean of the steady-state current (Figure 3.A) in order to capture the experimental trial-to-trial variability of non-spiking neuron responses, which is known to be larger than that of spiking neurons [36]. Finally, Figure 3.B shows the resulting voltage dynamics of the model against an example of experimental voltage traces. Again, the model is not intended to reproduce the specific shown traces, but to capture the trial-to-trial variability of the neuron responses. As explained in [32], although the model slightly diverges from the experimental voltage traces, it preserves the main qualitative features of the neuron dynamics. Indeed, the RIM model is depolarized or hyperpolarized in a smooth manner under current clamp, while the AIY model is more sensitive to hyperpolarization than depolarization inputs with transition point around −30 mV, as in the experimental membrane potential behaviors. Likewise, the voltage jump of the AFD model to its up-state occurs for the same values of current injections as in the biological neuron.

Mathematical analysis of the non-spiking model
This section aims at studying the phenomenological non-spiking model (2.1) presented in Section 2. The first part analyses the stationary points of the model, while the second studies its bifurcation dynamics. Moreover, we apply the model with the parameters of the RIM, AIY and AFD neurons ( Table 2), representative of known types of behaviors of non-spiking cells: near-linear and bistable. All numerical simulations of the models presented in this article are performed with the function ode of the Scilab software. As a first step, some biological hypotheses are formulated and used throughout the manuscript based on experimental studies.
Biological hypotheses throughout the manuscript. Biological cells have limited tolerance to electrical stimulation. This is even more true for the C. elegans neurons due to their small size. Beyond a certain threshold, the stimulus bursts the cell. Therefore, we define the range of the neuron functioning, noted I, as follows: where m = −15 pA and M = 35 pA as in [33].

Stationary point analysis of the model
Let a, b, c and d be the parameters of the model (2.1). In the following, let A and B be defined by Moreover, let p and q be defined by and I 1 and I 2 as follows: Then, the following proposition holds.
where p and q are defined in (3.4). Moreover, letĨ be defined by (3.3), and I 1 and I 2 be defined by (3.5). Then: has a unique stationary point V * (I) for all I ∈ I given by • for I = I 1 and I = I 2 , the model (2.1) has two stationary points V 1 * (I) and V 2 * (I) given by • for all I ∈]I 1 ; I 2 [, the model (2.1) has three stationary points V 1 * (I), V 2 * (I) and V 3 * (I) given by Proof. The stationary points of the model (2.1) are given by solving Let us denote V = X − α/3 with α = b/a. Then, the equation (3.11) can be rewritten as where p and q are defined in (3.4). The number of roots of the third-degree polynomial (3.12) depends on the sign of the function ∆ defined in (3.7) which is its discriminant. Moreover, ∆ is a convex parabola, so that the number of stationary points of the model (2.1) is determined by the sign of the minimum of the function ∆. We can rewrite ∆ in the following form: where A and B are defined in (3.2), and C is defined in (3.6). Then, we have so that the minimum of the function ∆ is reached inĨ = B/2A. From the Cardan-Tartaglia method, if ∆(Ĩ) > 0 then the model (2.1) displays a unique stationary point given by the equation (3.8). Moreover, if ∆(Ĩ) ≤ 0, then the model (2.1) has one, two, or three stationary points depending on the value of I.
The critical values of I at which the number of stationary points changes are those satisfying ∆(I) = 0. These are given by I 1 and I 2 defined in (3.5). Therefore, for all I < I 1 and I > I 2 , ∆(I) > 0 so that we have only one stationary point given by (3.8); for I = I 1 or I = I 2 , ∆(I) = 0 so that we have two stationary points given by (3.9); for all I 1 < I < I 2 , ∆(I) < 0 so that we have three stationary points given by (3.10).
Biological interpretation. If ∆(Ĩ) > 0, then the model (2.1) displays a near-linear behavior defined by smooth depolarizations or hyperpolarizations from the resting potential ( Figure 1.B, left). If ∆(Ĩ) < 0, then the model (2.1) exhibits a bistable behavior defined by a nonlinear transition between the resting potential and a depolarized potential (Figure 1.B, right). In this case, I 1 and I 2 defined in (3.5) are the injection current thresholds from which the neuron jumps to its down-and up-state, respectively.
Application to the RIM, AIY and AFD neurons. Here we determine the number and the expression of the stationary points of the three neurons under study (RIM, AIY and AFD). To do this, we compute ∆(Ĩ), defined in (3.7), with the parameters a, b, c and d given in Table 2, determined in a previous work [32] to reproduce the qualitative behavior of these neurons. The result for each ones is summarized in Table 1 and shown in Figure 4. Table 1. Minimum of the function I → ∆(I) for each neuron.
We observe that ∆(Ĩ) > 0 for the RIM and AIY neurons so that the following corollary holds from Proposition 3.1. Moreover, ∆(Ĩ) < 0 for the AFD neuron. We compute I 1 and I 2 defined in (3.5) and we find I 1 = 2.62 and I 2 = 3.12. Then, we have the following corollary: • for I = I 1 and I = I 2 , AFD has two stationary points V 1 * (I) and V 2 * (I) given in (3.9); • for all I ∈]I 1 ; I 2 [, AFD has three stationary points V 1 * (I), V 2 * (I) and V 3 * (I) given in (3.10).
Stability of stationary points. Consider I such as F(V * (I), I) = 0 where F is the right-hand side of the model (2.1). The first-order derivative of the function (V, I) → F(V, I) with respect to V is given by Then, for each neuron, we compute from Cardan-Tartaglia formula the equilibria V * (I) for each I ∈ I by 0.01pA, as well as the value of the function (3.14) in V * (I) to determine its nature. As shown in Figure 5, the sign of the function (3.14) is always negative for RIM and AIY, showing that their unique equilibrium point is always stable. Regarding AFD, it depends on the value of I. For any I < I 1 and I > I 2 , the unique equilibrium point is stable. For I = I 1 and I = I 2 , we have one stable and one non-hyperbolic equilibrium point. For any I ∈]I 1 ; I 2 [, we have two stable equilibria and one unstable equilibrium point.  14) with respect to V is always negative for all I ∈ I, so that their unique equilibrium point is stable. For AFD, the sign is negative for I < I 1 , so for these values of injection currents, the single equilibrium point is stable. For I = I 1 and I = I 2 , we have two equilibria, one stable and one non-hyperbolic. For I ∈]I 1 , I 2 [, we have three equilibria, two stable and one unstable.

Bifurcation analysis of the model
The analysis of the stationary points carried out in the previous section shows the existence of two non-hyperbolic equilibrium points for bistable dynamics of the model (2.1). The following proposition gives information on the dynamics of the model (2.1) near these points. Then, a saddle-node bifurcation occurs at I = I * . Moreover, the normal form of the model (2.1) near the stationary point V * (I * ) is given by Proof. For the convenience of readability, we set τ = 1.
• We carry out the following change of variables: The model (2.1) in the new coordinates becomes Therefore, (0, 0) is a non-hyperbolic stationary point of the model (3.16).
• The third order term in y is negligible in 0 compared to the second order and first order terms in y. Then, we obtain By assumption, −3aV * − b 0, so that • Finally, we rescale the Eqs (3.18) by considering the change of variable Application to the AFD neuron. In the previous section, we highlighted two non-hyperbolic equilibrium points V * (I 1 ) and V * (I 2 ) in the case of the AFD neuron (see Figure 5 as an illustration), where I 1 and I 2 are still defined in (3.5). Indeed, ∂F(V * (I 1 ), I 1 ) ∂V = 0 and ∂F(V * (I 2 ), where F is given in (3.14). From Proposition 3.2, the dynamics of the model (2.1) near the nonhyperbolic equilibrium V * (I 1 ) is given by the following normal form: whereas near V * (I 2 ), the normal form reads as Finally, we obtain the bifurcation diagram displayed in Figure 6.

Non-spiking neural network
Experimental measurements of the non-spiking neuronal activity in freely C. elegans worms show that ON and OFF network states underlie physiological behaviors [28][29][30]. This section aims at understanding how such correlated up-and down-states of the membrane potential of the neurons can emerge from the intrinsic neuron dynamics and their couplings. To that end, we first present the C. elegans synaptic connections and their modeling, before determining mathematical conditions informing whether the ON-and OFF-states result from the intrinsic neuron dynamics or the synaptic activity.

C. elegans synaptic connection
Anatomical reconstruction of the C. elegans nervous system reveals both chemical and electrical synaptic connections [37]. Synaptic transmission at chemical synapses occurs through the release of neurotransmitters similar to those in more complex vertebrates, such as GABA [38], glutamate [39], acetylcholine [40], and some amines (serotonine, dopamine, etc.) [41]. Regarding C. elegans gap junctions [42], i.e. electrical synapses, numerous works demonstrate their functional equivalence to mammalian gap junctions (we refer to [43] for a recent review). Chemical and eletrical synapses often coexist in same target neurons.
From a general viewpoint, neurons can transmit information through various modes of neurotransmitter release. In the case of C. elegans, experimental biological studies show that neurons release neurotransmitters in a graded manner [44][45][46][47][48]. In addition, there is evidence of a graded release in Ascaris lumbricoides [49,50], another extensively studied nematode related to C. elegans. More precisely, the graded synaptic connection is defined as follows: Definition 4.1 (Graded synaptic connection [45]). Graded synaptic connection is one in which action potentials are required neither for transmitter release nor for postsynaptic signalling, and in which incremental changes in membrane potential in the soma of the presynaptic neuron cause incremental changes in membrane potential in the soma of the postsynaptic neuron.
Graded release of neurotransmitters in C. elegans has been firstly modeled in [51] based on data from Ascaris lumbricoides [49,50] adapted to the particular case of C. elegans neurons. In the following, we will refine some parameter values based on recent experimental evidences from C. elegans.

Modeling of synaptic connections
Each instance of the models composing the network follow the general equation • V i denotes the membrane potential of a neuron i; • E i j the reversal potential for the synaptic conductance, set to 0mV for excitatory synapses and −48mV for inhibitory synapses [51]; • V j a presynaptic membrane potential; • g gap i j the constant conductance of the gap junction. Recent experimental works [52,53] directly recorded trans-junctional currents between C. elegans neurons, from AVA to VA5, giving a parameter value to about 0.4nS. Therefore, g gap i j will be set to 0.4; • g i j the synaptic conductance of the postsynaptic membrane which is gated according to presynaptic potential V j of a sigmoidal manner. Therefore, t → g i j (t) is described by the equation where τ i j is the synaptic constant time and V → g i j∞ (V) is a sigmoidal function relating preand postsynaptic membrane potential. In the C. elegans worm, τ i j 1, so that the activation of synaptic conductances can be considered as instantaneous [51]. Therefore, the synaptic conductances depend only on the presynaptic membrane potential V j , i.e. g i j (t) = g i j∞ (V j ) Finally, the expression of g i j∞ is given by where g i j is the maximal postsynaptic conductance, V i j 1/2 satisfies g i j∞ (V i j 1/2 ) = g i j /2, and V slope the slope factor where smaller values of V slope lead to a sharper g i j∞ .
In the following section, we carry out a mathematical analysis of the dynamics of coupled neurons based on the formalism described in this section.

Mathematical analysis of coupled non-spiking neurons
Let N be the number of neurons in the network, and i ∈ {1, . . . , N} a given neuron. Based on the previous section, the dynamics of the membrane potential of the neuron i is given by . Therefore, the system (4.2) can be written as follows: Stationary point analysis of the coupled models. From (4.3), the stationary points of the membrane potential of the neuron i verify where V j and I ext i are considered as parameters of the system of the neuron i. Then, based on the previous section, the discriminant of the third-degree polynomial (4.4) is given by where C i and D i are defined in (4.3). Then, the following proposition holds.
Proposition 4.1. Let ∆ i be defined in (4.5), p i and q i be defined in (4.6), and V j and I ext i be the inputs considered as parameters. Then: • if ∆ i > 0, the membrane potential of the neuron i has a unique stationary point V * i (V j , I ext i ) given by • if ∆ i = 0, the membrane potential of the neuron i exhibits two stationary points V 1 * i (V j , I ext i ) and • if ∆ i < 0, the membrane potential of the neuron i has three stationary points V 1 * i (V j , I ext i ), .

Biological interpretation.
Beyond informing us about the stationary points of a coupled neuron, Proposition 4.1 tells us about how the coupling impacts the instrinsic dynamics of the coupled neurons. Indeed, if ∆ i > 0, then the neuron i governed by the equation (4.2) displays a near-linear intrinsic dynamics, while it exhibits an intrinsic bistable one if ∆ i < 0. Application on a canonical example of two coupled neurons. Here we implement a non-spiking neural network of two coupled neurons through a gap junction and an excitatory chemical synapse (Figure 7.A). We consider the AFD bistable neuron as the presynaptic sensory neuron, noted A, and the RIM near-linear neuron as the postsynaptic neuron, noted R. It is worth noting that the behavior of the AFD and RIM neurons represents the two qualitative types of non-spiking behavior (bistable and near-linear), and are representative of two major quantitative classes of C. elegans neurons [33].
The mathematical expression of the network is given by The dynamics of the membrane potential V R of the RIM neuron shows up-and down-states when coupled with a bistable presynaptic neuron (Figure 7.B, left), as previously reported in [32]. These up-and down-states are apparent in the bimodal distribution of the membrane potential (Figure 7.B, right), computed using the distplot function of the Seaborn library in Python. This is in agreement with its bistable intracellular calcium dynamics when exposed to odor stimuli [30]. This shows that the connectivity strongly affects the behavior of RIM: the membrane potential of a near-linear neuron can have bimodal dynamics when coupled with a bistable neuron. The question that then arises is whether these up-and down-states result from the intrinsic bistability of the cell, or whether they reflect bistable synaptic inputs. In other words, we want to understand whether the connectivity has modified the intrinsic nature of the RIM neuron, from a near-linear type to a bistable type. To address this issue, we use the Proposition 4.1. As a first step, the equation on the postsynaptic RIM neuron is rewritten in the form (4.3): As the chemical synapse is excitatory, we have E RA = 0 so that RA V A Therefore, we plot in Figure 8.A the discriminant ∆ R given in (4.5). We observe that ∆ R > 0 for any value of the presynaptic membrane potential V A . According to Proposition 4.1, the RIM postsynaptic neuron exhibits a unique stationnary point. This allows to conclude that the nature of the intrinsic dynamics of this coupled neuron remains near-linear despite the bimodality of its membrane potential. Finally, Figure 8.B shows the bistable dynamics of the synaptic activity which confers a bimodal membrane potential dynamics to the RIM neuron.
The same analysis has been conducted with the AIY neuron as the postsynaptic neuron, and AFD as the presynaptic neuron. As can be seen in Figure 10, the same results are obtained: the apparent bimodality of the AIY voltages is not the result of an intrinsic change in its nature, but of bistable synaptic inputs.

Conclusion and discussion
Experimental observations through calcium imaging in freely moving animals (e.g. in the C. elegans worm) show that non-spiking neural networks exhibit ON-and OFF-states. This paper aimed at understanding how such states could emerge from the intrinsic dynamics of neurons and their couplings. To that end, we used a recently developed simple model of non-spiking neurons from which a mathematical analysis of its dynamics were carried out. In particular, we determined explicit expressions of the stationary points of the model depending on its type, near-linear or bistable. In addition, the bifurcation dynamics of the model were studied, with the determination of its normal form. Then, we determined mathematical conditions for determining whether the ON-and OFF-states in non-spiking neural networks result from the intrinsic neuron dynamics or the synaptic activity. On this basis, we highlighted that the bimodality of the membrane potential of coupled neurons could be a consequence of bistable synaptic inputs, and not of intrinsically bistable neuron dynamics contrary to what one might have expected. Further, we suggested that up-and down-states in nonspiking neural networks could be driven by bistable neurons which dictate their behaviors to their postsynaptic neurons. Such a configuration, with bistable neurons sending electrical signals to nearlinear (and bistable) postsynaptic neurons, is ubiquitous in non-spiking neural networks. Indeed, many sensory, inter-and motorneurons in C. elegans are bistable such as ASER [9], RMD [54], AWC [55], ASH [56], AFD [33], AWA [33], AIA [57] and numerous non identified neurons [9]. In the same way, the bipolar [58] and horizontal [59] cells of the retina are bistable and send eletrical signals to ganglion and amacrine cells respectively [8]. Therefore, our study gives a rigorous framework for better understanding how the behavior of coupled non-spiking neurons depends on their intrinsic dynamics and their couplings. In particular, it allows to gain insight into how intrinsic dynamics and couplings of single non-spiking neurons impact neural network behaviors.
On the importance of distinguishing between bimodality and bistability. Here we stress the importance of distinguishing between the terms bistable and bimodal. Bistability results from the intrinsic properties of the neuron, whereas a bimodal distribution of the neuron's voltages does not necessarily. A bistable neuron inherently displays a bimodal distribution of its voltages, but a neuron can exhibit a bimodality of its voltages when it is coupled without being intrinsically bistable. Indeed, we showed in this paper that the RIM near-linear coupled neuron can display a bimodal distribution of its voltages without being transformed into a bistable neuron.
The implications on the network signal processing could be radically different whether the bimodality of the RIM voltages are induced by intrinsic bistability or not. Indeed, a modification of the near-linear type of RIM into a bistable type would endow it with new computational characteristics. In this regard, bistable neurons enhance the short-term memory capabilities of excitatory recurrent networks [60], stabilize the persistent neural activity observed in working memory systems [61], and reduce the effects of noise and distraction stimuli [57,62]. Therefore, if we want to understand the way nervous systems work to produce specific behaviors, we then must understand how synaptic activities influence the dynamics of neurons and their intrinsic nature. We argue that the current paper, through Proposition 4.1, makes this possible in the case of non-spiking networks.
Perspective. In this paper, we determined mathematical conditions to study the influence of synaptic activities on the neuron intrinsic dynamics. We applied them to canonical examples of two coupled neurons, representative of widespread configurations in non-spiking networks. However, an application on larger scale networks could be implemented in the future, for example on C. elegans networks for which we have the fully mapped connectome (i.e. the wiring diagram of its nervous system) [21,37,63].
In particular, a biological network of four pair of neurons displaying correlated up-and down-states has been well-studied in a biological viewpoint [30]. This network underlies a chemotaxis behavior, characterized by C. elegans reversal (reorientation) responses to odors. Modeling this network could be of significant importance to deal with many issues. For instance, a currently open question is the specific contribution of the gap junctions vs. chemical synapses in the emergence of the up-and downstates [30].