{'fraction_non_alphanumeric': 0.033295130623879635, 'fraction_numerical': 0.005554086220576567, 'mean_word_length': 4.728956228956229, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Spiking Neural Dynamical Drift-Diffusion Model on Collective Decision Making with Self-Organized Criticality

This article proposes a novel collective decision making scheme to solve the multi-agent drift-diffusion-model problem with the help of spiking neural networks. The exponential integrate-and-fire model is used here to capture the individual dynamics of each agent in the system, and we name this new model as Exponential Decision Making (EDM) model. We demonstrate analytically and experimentally that the gating variable for instantaneous activation follows Boltzmann probability distribution, and the collective system reaches meta-stable critical states under the Markov chain premises. With mean field analysis, we derive the global criticality from local dynamics and achieve a power law distribution. Critical behavior of EDM exhibits the convergence dynamics of Boltzmann distribution, and we conclude that the EDM model inherits the property of self-organized criticality, that the system will eventually evolve toward criticality.


I. INTRODUCTION
Self-Organized Criticality (SOC), a ground breaking achievement of statistical physics, has gained growing interest in neural firing and brain activity in recent years [1]. Bak's hypothesis [2] and recent studies [1], [3] suggest that criticality is evolutionarily chosen for optimal computational and fast reactionary purposes, and that the brain is always balanced precariously at the critical point.
Such critical dynamics emerge during the phase transition between randomness (sub-critical) and order (super-critical), and usually follow power-law distributed spatial and temporal properties [3]. A dynamic network system with the SOC behavior then has the potential to be spatial and/or temporal scale free [4] and to fast switch between phases and attain optimal computational capability. This offers a possible approach to model a decision making process.
The systems that exhibit SOC behavior are usually high dimensional and slowly driven, with nonlinearity properties [2], [4]. To this end, Brochini et al [5] have discussed the phase transitions and SOC in stochastic spiking neural networks. Also, Bogacz et al [6] demonstrated that standard Drift Diffusion Model (DDM) can be used for stochastic spiking dynamics and they relate DDM to a highly interactive "pooled inhibition" model. However, to authors' best knowledge, although the SOC has been recognized as a fundamental property of neural systems [1], there has yet to be a decision making model capitalizing the SOC property. In this paper, we incorporate the nonlinearity of Exponential Integrate-and-Fire (EIF) model to replace the stochastic spiking scheme in DDM proposed by [6]. Introduced by Fourcaud-Trocme et al [7], the nonlinear EIF model is experimentally verified to be able to accurately capture the response properties. It will be demonstrated that neural sampling and mean field branching can be derived with the Boltzmann distribution. The proposed Exponential Decision Making (EDM) model therefore reaches a set of absorbing states, and the corresponding global criticality follows power law distribution, thus attaining the SOC behavior.
The contributions of this paper are summarized as follows.
1) To the best of authors' knowledge, this paper is the first published work on modeling decision making processes with the SOC property. 2) A collective decision making model, i.e., EDM, is proposed to implement the DDM methodology on EIF spiking neurons.

3) A probability inference scheme on EIF sampling is
proposed, which extends an existing leaky integrateand-fire sampling method. 4) Mean field analysis of the connectivity of EDM is given, which exhibits global criticality. 5) Detailed analysis is given to reveal SOC behavior of the EDM model under criticality conditions. This paper is organized as follows. In Section II, literature review is given and the necessary background concepts are introduced. Then, in Section III, the dynamics of each agent as well as the EIF neural sampling are discussed. In Section IV, the collective behavior of the network system is analyzed. Section V provides convergence analysis and simulation results. Section VI concludes the paper and looks into future work.

Notation and Preliminaries
Here we use a classical directed graph representation G = (V, E, A) with nonempty finite number of nodes and edges. Specifically, V is the set of nodes, E is the set of directed edges, and A = [a ij ] is an adjacency matrix with weights a ij > 0 if, i, j ∈ E, an edge from node i to node j. Note that the assumed graph is simple, i.e., i, i / ∈ E, ∀i, with no multiple edges going in the same direction between the same pair of nodes and no self loop. In this case, the diagonal elements of A are zero. In addition, the Laplacian matrix of G is denoted by L.

II. BACKGROUNDS AND RELATED INFORMATION
A. Self-Organized Criticality Self-Organized Criticality describes a self-tuned internal interactions that show critical dynamics in complex systems [2]. The interacting node groups are called active sites while the nodes that are less sensitive to the input are called inactive sites. Usually in the sand pile model, each agent has their own steep slope, which represents the membrane potential of the spiking neurons. When a certain threshold is hit and the sand in that specific area is steep enough Z local > Z critical , the avalanches will be triggered, which follow a power-law distribution of 1/f noise. Plenz and Beggs [8] observed a similar pattern of avalanches in the cortical neural electrical activity, which was the first evidence that the brain functions at criticality.
For the spiking neural network sense, if an agents activates too many neighboring neurons (super-critical), it leads to the massive activation of the entire network, while if too few neurons are activated (sub-critical), propagation dies out too fast [6].
In our case, the local rigidity level is expressed by two terms: the firing threshold of each agent, and the correlation between each two agents in the same local active sites.

B. Boltzmann Machine
Boltzmann Machines (BM) is a special type of stochastic recurrent neural networks based on non-stochastic Hopfield nets. In recent years, BM's property of binary output attracts more attention in both the theoretical neuroscience and the high dimensional parallel stochastic computation [9], [10].
BM is proven to be efficient for the models with connectivity properly constrained, to be specific, machine learning and probability inference are two major applications. In this paper, we apply neural sampling and show that the probability function of the gating current variable for the activation term follows the Boltzmann functions.
The global energy function of the Boltzmann Machine is defined as: where E is the global energy, h ij is the connection strength between Unit j and Unit i, s i ∈ {0, 1} is the state of Unit i, and b i is the bias of Unit i.
As for the Boltzmann distribution, the probability that the ith unit is on is where T is the temperature of the system. The probability is calculated using only the information of the energy difference from the initial state, and the temperature of the current time. In our framework, we consider the term − E i /T as a logistic activation function as in [11], [12]. It has already been shown in [12] that some stochastic neurons sample from a Boltzmann distribution. Ideally, after a long running period and without further inputs, the probability of a global state will not be affected by other terms, i.e., time constants and conductance values in the EIF model. At this stage, the system is at its "thermal equilibrium", that converges to a low temperature distribution where the energy level hovers around a global minimum.
This feature presents a similar behavior as SOC if we consider the criticality as this thermal equilibrium that energy level fluctuates around. Also, the log-probabilities eventually becomes a linear term, which helps us simplify the exponential term in the EIF model. Further discussion will be given in later sections.
Moreover, the neural sampling technique in the later section incorporates the Boltzmann machine according to some local switching, with conditional probability integrated, the multivariable Boltzmann joint distribution has the form [13] where P m stands for the probability of state m, m is the energy at state m, k is a constant, and M is the total number of the states.

C. DDM
Drift Diffusion Model has been applied on Two-Alternative Forced-Choice (TAFC) task in an extensive amount of work (see [6], [14] for instance). The fact that DDM integrates the difference between two choices according to one or two threshold makes it possible to describe a decision making process in a spiking neural network.
In the pure DDM, the accumulation of the unbiased evidence has the form where dx represents the changes in difference over the time interval dt, g is the increase in evidence supporting the correct choice each time, w is the independent, identically distributed (i.i.d.) Wiener processes, and β is the standard deviation. The probability density P (x, t) is normally distributed with mean gt and standard deviation β √ t. Since the second term in (3) is represented by a standard Wiener process that describes the noise, it is common to consider dx in DDM to be the change in membrane potential within a certain amount of time [15].
Here we consider a network system with N agents. Each agent relates the DDM model to the non-stationary dynamics of the firing of an EIF spiking neural model. While the forced-response protocol is usually considered, we follow the free-response protocol, that each consecutive fires determine the range of the time interval. The common assumption made for this equation usually considers g > 0 to support the first choice, and g < 0 for the other [6]. The term g can either be a constant for inactive nodes, or a function for active nodes that depends on membrane potential.
While (3) only describes the dynamics of a single DDM system, we need extra terms to capture the impact from neighbors. We have the following stochastic diffusion process with an initial condition x 0 : where α(t) is a measurable gain function that models the external input to accelerate the potential increment, and linear drifting term g(t) represents the dynamic drifting variable of the node itself. In this regime, we have transferred our model to a Ornstein-Uhlenbeck (OU) process, which is known to be the solution to the famous Fokker-Planck equation.
Here, to further simplify the model, we may eliminate the afterhyperpolarization, that is, let x 0 = 0.
For the model proposed above, it is possible to receive the spike generations with arbitrary shape, i.e., different spiking time intervals and different incremental speed of membrane potential. With a proper defined activation function, which will be discussed in Section III, the behavior of each single stochastic diffusion process can be bounded. Before further discussing into the individual dynamics and their boundedness properties, we need to look into the specific local dynamics by applying a most commonly used neuron model.

D. Generalized Exponential Integrate-and-Fire (EIF) Model
Exponential integrate-and-fire (EIF) is a well developed biological neuron model introduced by Fourcaud-Trocme et al [7] as an extension of the standard leaky integrate-andfire model. As concluded in several studies [11], [12], EIF is a suitable simple model for very large scale network simulations. For the generalized EIF [12], arbitrary spike shapes are allowed and gated currents usually reach a steadystate with nonlinear voltage activation function.
The EIF model holds a nonlinearity property consisting of a linear leakage term combined with an exponential activation term, which follows a simple RC-circuit dynamics before V , the membrane potential, reaches a set threshold V T . After reaching the threshold, it can be considered that the neuron has fired, and its membrane potential is then set to a resting voltage, V R , approximately −60 mV [11], [12] or 0 mV [15] by different assumptions.
The dynamics of the membrane potential is given by In this equation, C is the membrane capacitance, V T is the membrane potential threshold, ∆ T is the sharpness of action potential initiation, or slope factor, V L is the leak reversal potential, is the conductance, and I ion is input current. While I ion only represents synaptic current in Fourcaud's model, here we have extended the ionic current by summing up input current I neib from neighbors with connectivity, external noise current I noise that integrates i.i.d. Wiener process, and synaptic input I syn that incorporates the drifting term, serving as a bias. We have I ion = I neib +I noise +I syn , where the term I neib takes identical form as leakage current in the first term of (5), while I syn = syn Γ(V − V syn ) represents the slow voltage activated current with a gating variable Γ = Γ(t). The term Γ ∞ = lim t→∞ Γ(t) can be used to describe the instantaneous activation.
Here, we alter the usually constant conductance , and change it to a function of V and t. Multiplying dt to both sides of the (5), we now have It is clear that most terms in (6) have very similar forms to those in (4). As most current terms do not need to be altered to fit in (4), the conductance term can change over time and become a function. Functions and α are sometimes interchangeable. However, the exponential term can be tricky to work around, and we will talk about it in the later part after done discussing absorbing state.
Henceforth, for simplicity purpose, we call this EIF and DDM combined model as Exponential Decision Model, or EDM for short.

III. DECISION MAKING DYNAMICS
The commonly discussed decision making process is an adaptive behavior that makes use of a series of external input variables and then leads to an optimal or sub-optimal choice of action over other competing alternatives.
We begin by discussing this optimal decision rule. There are two thresholds z i in DDM model, with the same magnitude but different signs to represent different choices. In EDM, we consider this choice to be optimal if the threshold of correct choice is reached, or sub-optimal, if our expectation value E ends up with the same sign as correct threshold but with smaller magnitude. Now we start solving the OU process described in (4). Lemma 1: The solution of the collective decision making system in (4) is given by which has a similar form as in [16], with the updated expectation Proof: Let φ(t) be a fundamental solution matrix. Also, let Y be Then Y has the stochastic differential equation which implies Furthermore, combining all above together, it is straightforward to conclude (7).

A. Behaviors
Here we consider a very special but common network system.
Assumption 1: For a high dimensional dynamical network system with N agents described by (4), each unit integrates an inward stimulus α i and receives signals I neib and noise β i dw j i from local j neighbors. The dynamics of w j follows Correlated Brownian motions, with standard correlation col ∈ (−1, 1).
In Reference [17], the authors proposed a Itô consensus S.D.E formalized by N 2 -dimensional standard white noise, and expanded the right hand side terms in (3) to a matrix form with graph theory, i.e., α(t)Lx(t). Different from their encoded gain function, we now assume that all the state information is available to others.
Assumption 2: For the considered network system with N agents, the state of each agent is observable by others.
Extending (4), then we have the dynamic equation for each agent with neighbor's dynamics added where K is the total number of agents connecting the agent i, l ij are elements in the Laplacian matrix L, y ji denotes the observed membrane potential of jth agent by the ith agent.
In Reference [18], Charalambos et al have shown a method of achieving optimization problems under a reference probability measure by transferring continuous and discrete-time stochastic dynamic decision systems, via Girsanov's Measure Transformation. To this end, we have the following claim.
Proposition 1: The collective stochastic dynamic decision system with common team optimality can be transformed to the equivalent static optimization problem with independent distributed sequences. And under the reference probability space, states and observations are independent Brownian motions.
Consider a series of d inputs, X i (t) = [x i1 (t), ..., x id (t)], for N agents, each of which has 2 states The firing during τ ref is sometimes called absolute refractory period. Such series of inputs can be modeled as discrete decision making scheme, which we have applied the free-response paradigm upon, with adaptive prescribed time interval τ x for each agent in the network system. Equation (9) then becomes the following distributed protocol dX = L·α(X(t), t)· Y(t)−X(t) +g(t) dt+β(X(t), t)dw. (11)

B. Sampling with EIF
Recall from Proposition 1 that, since the Brownian motions are indpendent of all other team decisions, it opens up the possibility for Markov chain related methods. For efficiency and flexibility purpose, Markov Chain Monte Carlo (MCMC) has been applied in sampling the spiking network of neurons [10].
In Reference [12], Richardson has shown the equilibrium value of a slow driven voltage-activated current gating vari- is an adaptive time constant characterized by different voltage values, and Γ ∞ is the equilibrium value. Then we have where Φ is the sigmoid function that digests a membrane potential function into a probability density function with range [0, 1].
Here, the equilibrium term Γ ∞ holds a very similar form to Boltzmann probability distribution as in (1). It is clear that (12) is a slowly varying function, which can be proved simply by applying the definition. Now we can use the property of the slow varying function to deal with the exponential term in (6). The Karamata representation theorem is one of the mostly used property of slow varying functions that transfers a function into a general exponential form. In our case, Γ ∞ is expressed as: Γ ∞ = exp (Γ) + Γ B ε(t) t dt , for some B > 0, where (Γ) is a bounded measurable function converging to a finite number, and ε(t) is a bounded measurable function converging to 0.
Here, since the exponential term is only a property of membrane potential increment that adds to nonlinearity, the potential accumulation of each agent does not affect collective network decision as much during the firing period (t − τ ref , t]. And for the absolute refractory period, the neuron model is guaranteed not to fire. For such a piece wise continuous function, if we only consider the time interval to be one firing, then the variable V can be bounded. Taking out the exponential term in (6), we have where t end is the end of the refractory time, V end is the membrane potenial at end of the refractory period, and V e is the membrane potential incremented by exponential term only. Since converges to a finite number and ε(t) converges to 0, the limit of the exponential term converges to a constant C during the refractory period. It can be thought of as entering an absorbing state where its behavior at infinity is very similar to the behavior of converging to infinity. Therefore, at the absolute refractory period, we can safely ignore this exponential term in (6), and treat it just as other linear terms. For simulation, we have used the same parameters as the work done by Barranca et al [11]. As shown in Fig. 2, this monotonically increasing activation function has a sigmoidal general form, with an upper bound 1 and lower bound 0. In some cases, as mentioned above, we may set the initial condition to be 0 mV to eliminate the afterhyperpolarization. For our simulation, we have removed such constraints to show that the dynamics of the membrane potential can be bounded by the activation function with arbitrary initial conditions. This is because the higher threshold in the EIF model always requires higher voltage to be breached, and the logistic activation function describes such positive correlation well enough.
In fact, very similar dynamics of activation function has been capitalized in the neural sampling frame work [10], [12], [19]. Thus, our collective decision making model can be thought of as a network consisting of N agents (or neurons) sampling from a probability distribution p using the stochastic dynamics carried from the Drift-Diffusion model. Proposition 2: The firing activity of the generalized EIF model, which represents each agent in a collective DDM network system, follows a Markov chain process.
With the information-coded signal from each DDM agent of the system, in other words, the firing information within the time interval (t − τ ref , t], the neural sampling follows conditional probability distribution, and most of time is Boltzmann distributions.
For each agent, we treat the collective behavior of connected nodes as a accelerator/damper. For instance, if the agent is surrounded by nodes with higher membrane potential, it receives more current than normal drifting, and vice versa.
IV. COLLECTIVE BEHAVIOR For most multi-agent dynamic network systems, it is common that new communication links can be established over two agents with no previous connection. In the previous sections, we have assumed that the connection is known at each state. Now we define the coupling and connection behavior among the agents in the system.

A. Coupling and Connectivity
Assumption 3: For the system described in this paper, agent i forms at most K outward links randomly at t = 0. When s i = 1, Agent i tries to establish new connections with new neighbors, for instance, connecting to Agent j with the coupling probability P ij depends on the voltage difference. When successful, the equal number of previous connections are lost according to a decoupling probability function Q. When s i = 0, Agent i will not actively modify its neighboring connections.
It is worth pointing out that both P and Q follow a sigmoid (or reverse sigmoid) relation Φ (or 1−Φ). For function P, the greater the difference in membrane potentials has, the higher the probability is. The situation is reversed for function Q, the greater the difference in membrane potentials has, the lower the probability is.
Then we have the following equations where K i is the number of current established connections of node i. In the case of probability values that are less than zero or greater than one, we simply set them to 0 and 1 respectively for mean field analysis. The negative probability is also provided in Fig. 4 . Different from the common equation for branching probability, our process is not a treelike process and the maximal connectivity is defined by K. Therefore they are Markovian processes with respect to the number of connected nodes.

B. Mean Field Analysis
In the mean field analysis, the individual drifting variable usually follows a distribution with the average. In our case, we assume that the expanded term follows the average distribution of g i /W , where W can be considered as the sum of the total drift. Here g i is still a gain function representing Recalling the term "active site" described in Section II-A, the local active site for each agent can be thought as the local field: F i = − N j a ji g j /W , where a ji is the elements in an adjacency matrix, then the global field has the averagē F = N i (F i /K) /N . Therefore, with regard to a sequence of d number of inputs, the probability of a single DDM having the state s i = 1, represented by a single generalized EIF neuron, has approximately the following probability with respect to the field.
where S is some partition function, b i is the bias term that supports the correct choice, γ is a thermodynamic beta in Boltzmann factor with the form γ = 1/ (k B F i ), and k B is a Boltzmann constant. Putting aside the coupling strength, the mean activity of this network is measured as N i=1 s i /N . However, for such a spiking neural model, it is not always practical to have normally distributed probability density function, and in fact, most of these processes are stochastic with certain thresholds or even highly constrained. That being said, we would have a stochastic Itô based integral for probability density for the mean activity: H = ∞ −∞ Φ(X)p(s i |X)dX.

V. CONVERGENCE ANALYSIS
There are two main evidences for systems presenting the SOC behavior [2], [4], the first is the power law distribution, and the second is the critical dynamics, which we consider as converging to absorbing states in the EDM model. In this section, we expand the results from Sections III and IV, and examine the global convergence behavior of the collective EDM model. We then provide both pieces of evidence to show the EDM system has the SOC behavior.

A. Global Criticality from Local Dynamics
Recall that in SOC, active nodes trigger self avalanches when a threshold is reached, and update the information of all connecting active nodes. Also, nodes in nearby inactive sites will be communicated to establish more connections if needed.
Moreover, each local active site has their own dynamics of reaching out to other active or inactive sites. Inspired by the work of Harris on the theory of branching processes [20], we use a branching parameter σ that captures the subsequent activity of connectivity triggering or dying-out [21]. The local branching ratio and global branching ratio have the form respectively. As discussed in [20], [21], the system exhibits criticality at σ = 1, and is sub-critical (super-critical) for σ < 1 (σ > 1).
In our simulation, sub-critical dynamics are characterized by low potential and rapidly decaying agent neuron firing distributions, while super-critical dynamics are characterized by high potential and slowly decaying firing activities. Critical dynamics are characterized by firing activity that follows power-law distributions.
As shown in Fig. 3, it can be easily recognized that the collective behaviors of the firing density in the EDM model follow a power law distribution, D(z) ∼ z Λ with different cluster size z and scaling factor Λ < 0, which is a primary evidence supporting the SOC behavior [4].

B. Absorbing States
Recalling the avalanches described in self-organized criticality, the system keeps tuning itself to one of many metastable states, which commonly have lifespans shorter than ground state and longer than excited states [2]. And without further inputs, the distribution reaches meta-stability, a very special energy well that is able to temporarily trap the system for a limited number of states.
This can be modeled as Boltzmann distribution with global energy level in a simulated annealing system from any initial conditions. With the results from (13) and (14), as well as the exponential property of the generalized EIF system, the Lyapunov based semistability [22] can be achieved with great potential, but due to the length restriction, this concept will not be discussed in this paper. Nevertheless, we are going to show that the system described in this paper converges to some absorbing states.
In mean field theory, as being discussed in [21], absorbing state becomes unstable when the probability of a node creating connection with neighbors is greater than 1/2. In  Fig. (4a), as the firing probability increases with the membrane potential, σ converges to 1 with proper connectivity constraints. Cases with different numbers of active and connected neighbors are shown in a network system with N = 10. Fig (4b) removes the constraint that cumulation of local branching ratio in each iteration caps at 1. It is clear that as the number of connected neighbors increases, the network system enters the active/loading phase first, and then evolves to the dissipation/absorbing phase. The system clearly shows SOC dynamics, that is σ = 1 at both minimal and maximal connectivity. our case, this can be thought as the coupling probability P > P critical = 1/2. Lemma 2: The attractor of the system is a set of discrete states.
Proof: If a non-conserving system, such as the drift and diffusion based EDM model described in this paper, has shown a temporary stable configuration after the avalanches, then the system is at least at a critical point. The critical and super-critical session are usually slow driving [21]. So there must be a drift load and a diffusion dissipation fluctuating to keep all the nodes in the system from either forming active sites, or staying quiescent completely.
Thus, if the system presents the thermodynamic behavior as Boltzmann distribution with simulated annealing, there can exist infinite numbers of infinitesimally varied absorbing states in thermodynamic limit. As for a finite number of total states, this absorbing phase becomes a set of discrete states.
This obeys another property of SOC, that is, the dynamical system with a critical point as an attractor, is able to keep itself at the critical point between two phases, which in our case, are the active phase and the absorbing phase.
Proposition 3: Meta-stable states generally hold more energy than the ground states, and less energy than the excited states. Therefore, SOC is the process of the EDM model losing global energy, and falling into a certain set of absorbing states, regardless of guaranteed stability.
Each agent in the EDM model is essentially drift diffusion terms taking input variable from EIF markup. It is clear that when the individual thresholds are reached, the agents will initiate the spike and send information-coded signals (current) to connected nodes or nodes with great probability of establishing connectivity. The exited states usually carry higher membrane potential than the incremental states, and then the fired node resets its membrane potential to V R and enters a refractory period. Therefore, the states during this absolute refractory period τ ref can be considered as comparable meta-stable states that trap the dynamics of each node for τ ref . And since we have proved Lemma 2, the system self-organizes itself, instead of fine-tuning, to a small region of absorbing criticality, or in another word, metastability.
Since the system is slow varying at the absorbing state, the fundamental solution is independent of t, that is, α(X(t), t)) = α(X). And (7) becomes At such slow varying, time independent absorbing states, it is natural to assume that the dimension of β is only 1.
Using the corollary in [16], φ(t) = exp t (t0) α(η)dη . Then we can further turn Theorem 1 to be The threshold presents the local rigidity level. Also, the convergence dynamics of these absorbing states show Boltzmann distribution as well, i.e., P i (V x |s = 0) = e −Xi/k BF / N j=1 e −Xi/k BF . As shown in Fig. 4, the branching pattern of the EDM model across multiple fires follows the SOC behavior and will eventually evolve to a certain set of absorbing states, known as the recurrence sets. Also, these stationary distributions of network states can be convergent from any initial state. Therefore, without the presence of a proper controller, the system fine-tunes itself, and then converges toū = lim t→∞ t−1 0 E(X(t))/t.

VI. CONCLUSION AND FUTURE WORK A. Conclusion
In this work, we have proposed a collective decision making model with a specific type of spiking neurons, exponential integrate-and-fire (EIF). Our method is based on the well-known Two-Alternative Forced Choice Task solver-Drift-Diffusion model. We recognize that DDM and EIF share very common terms in their dynamic equations, and the exponential term can be ignored during the absolute refractory period. We have derived the probability of each agent's firing based on a Markov chain conditional premise. Then the mean field theory is used to approximate the global criticality from local dynamics.
Both analytically and experimentally, we have found out that the global branching ratio follows a power law distribution and the EDM system eventually evolves to a set of absorbing states, which are two main evidences suggesting the Self-Organized Criticality behavior. The activation function follows the Boltzmann state probability and the convergence dynamics of absorbing states follow Boltzmann distribution as well.

B. Future Work
At this point, we have set up a detailed model that is ready to be expanded from different aspects. For instance, since theoretical and experimental studies have demonstrated that critical systems are often optimizing computational capability, it is promising to suggest that the system with the SOC behavior is both robust and flexible to ensure homeostatic stability.
In fact, due to the nature of absorbing states and criticality property, any initial conditions of a spiking network decision making system can converge and/or fluctuate around a set of states, potentially semistability [22]. Therefore, the convergence property of such a model can be useful for fault prescreening and is in a way, robust to quantified uncertainties.
Besides the theoretical analysis, there are potential applications as well. Bak has demonstrated that both traffic dynamics and brain dynamics exhibit similar criticality [2]. This opens up the way to extend our model to formulate realworld applications. With specific problem solving scenarios, it is natural to extend the proposed multi-agent spiking neural decision making system to a parallel distributed process, and eventually leads to a neuromorphic chip development.
Since the nature of SOC is to organize the system between two phases, the fast switching capability is useful for sensitivity analysis. Furthermore, we will also look into fully stochastic neurons such as Galves-Löcherbach (GL) model to incorporate different activation profiles and to increase the candidate choices as well.