In the context of a two-choice decision task, at zero coherence, the stimuli to the two selective pools are similar and hard to distinguish, so the system will make decisions randomly in a “coin tossing” manner. Therefore, it can be anticipated that there exists two symmetric attractors in favor of one of the two distinct choices (the left or right motion direction) in the binary decision system. When the parameters are specified as w+ = 1.61, μ = 58 and c = 0, we obtained a bistable attractor landscape (shown in Figure 2). The normalized steady state probability distribution P_ss(r₁, r₂) is firstly quantified by collecting the statistics of the system state in the “decision space” constructed by the selective populations averaged firing rate r₁ and r₂, and then the potential landscape can be mapped out by U(r₁, r₂) = −lnP_ss(r₁, r₂) (Sasai and Wolynes, 2003; Li and Wang, 2013a, 2014a) where U(r₁, r₂) is the dimensionless potential (see section 2 for details). Here, the potential U(r₁, r₂) does not correspond to the real energy in physical system, but a reflection of steady state probability distribution. In Figure 2A, the blue regions represent higher probability or lower potential, while the yellow regions indicate lower probability or higher potential. The landscape displays two symmetric basins of attraction. The two attractors with higher activity for one neural population (winner) and lower activity for the other (loser) are identified as distinct decision states, DS1 and DS2.

The decision making process corresponds to the transition from a spontaneous state where both S1 and S2 fire with similar low rates, to a decision state (DS1 or DS2) where they compete with each other in a winner-take-all manner, in response to the motion input. The attractor landscape provides an intuitive explanation for the decision process. Starting from an initial resting state, the system will evolve dynamically to arrive at a decision state. This process can be pictured as the system follows a path downhill to the bottom of the nearest valley which has the minimal potential energy (DS1 or DS2).

In response to the visual stimulus, subjects will make an initial decision. However, they continue to accumulate noisy evidences (Resulaj et al., 2009) or monitor the correctness of the previous decision (Cavanagh and Frank, 2014) after a choice has been made. Then the decision-maker may or may not change the mind, i.e., the initial decision is made to subsequently either reverse or reaffirm. Since the decision-making is endowed with randomness from the external stimulus and spontaneous activity fluctuations, the stability of an attractor state in landscape is only guaranteed up to a limited time-scale. In the case of changing minds, the noise will terminate the self-sustained pattern of activity and drive the system to escape away from the initial attractor and switch to the other although the stimulus input remains unchanged.

The neural system operates far from equilibrium. So we quantified the non-equilibrium probability flux map (denoted by the white arrows in Figure 2B) to measure the extent of the violation of detailed balance. The flux is calculated from the dynamical stochastic trajectories of neural activity (the time window used to calculated firing rate is 50ms) over a long time (see section 2) (Battle et al., 2016). We can see that the force from the curl flux drives the system away from one local attractor and guides the transition to the other. To see the effects of using different time windows for simulations, we showed typical trajectories for using different time windows (Supplementary Figure 1). We also estimated the landscape and probabilistic flux when the time window used to calculate firing rate is 20ms (Supplementary Figure 2). A major difference for different time windows is that for smaller time windows there are larger fluctuations (Supplementary Figure 1), which is also reflected by the landscape results showing that the landscape using 20ms as time window displays a little more variations (compare Figure 2B and Supplementary Figure 2). As for the curl flux, different time windows lead to qualitatively similar results (compare Figure 2B and Supplementary Figure 2), i.e., the magnitude of flux is larger in the region close to the basins, and the flux has certain curl direction.

The dynamical transition paths corresponding to the changes of mind process (green line for DS1 to DS2 and megenta line for DS2 to DS1) are shown in Figure 2C. Note the fact that the forward and backward paths are irreversible, which is the consequence of the non-equilibrium curl flux (Wang et al., 2008; Li and Wang, 2014a; Yan and Wang, 2020). It is worth noting that the curl flux and noise play different roles in the decision making system. Here, the role of the none-zero flux is breaking detailed balance, which is also the cause for the irreversibility of forward and backward transition paths (the green path and the magenta path are not the same, Figure 2C). If the flux is zero, the green and magenta path will duplicate completely, which corresponds to an equilibrium case (Wang, 2015). To see whether this irreversibility is due to stochastic effects, we used different stochastic trajectories to calculate transition path (Supplementary Figure 3) and obtained consistent curl directions for the irreversible transition paths for both bistable and tristable system. We also showed single-trajectory examples, which support the irreversibility for the transition path for both bistable and tristable cases (Supplementary Figure 4). Therefore, flux is important for the state transitions. In terms of the role of noise, we see that although the landscape and flux keep the same, the single sample for transition path varies due the noise effects (Supplementary Figure 4). The noise also contributes to the state transition for the barrier crossing process.

Furthermore, we found that the number of stable attractors on landscape changes under different parameter combinations. Figure 3A shows the tristable attractor landscape of the decision making network for w+ = 1.66, μ = 16 and c = 0. Compared to Figure 2A, the attractor that emerges at the bottom left corner is identified as the spontaneous undecided state with low activities of both selective neural populations, indicated by SS on the landscape. Similar to the case of bistability, the property of curl flux and irreversibility of transition path also exist in the tristable system (Figures 3B,C), illustrating the inherent non-equilibrium property of neural system. It is worth noting that the kinetic switches between two decision states go through the SS state (green and megenta path in Figure 3C), indicating that the system will erase the former decision (back to the resting state) firstly and then make another decision (Pereira and Wang, 2014). So the SS state can be treated as an intermediate state.

The intermediate state has been observed in various biological process. For example, in epithelial–mesenchymal transition, stem cell differentiation and cancer development, the intermediate state cell types play crucial roles in cell fate decision system governed by corresponding gene regulatory networks (Lu et al., 2013; Li and Wang, 2014b; Li and Balazsi, 2018; Kang et al., 2019). For working memory, the presence of the intermediate state significantly enhances the flexibility of the system to a new stimulus without seriously reducing the robustness against distractors (Yan and Wang, 2020). For the decision making function in this work, our landscape picture provides some hints on the roles of the intermediate state on mind changing. The mind changing process occurs in a step-wise way with the existence of the SS state. The system will switch to SS state firstly and stay there for a while, and then depending on the new accumulated evidences, decide whether transit to the other decision state (change of minds) or return back to the initial decision state (reaffirm previous decision) (Figure 3C). This demonstrates that the intermediate state may increase plasticity and robustness of perceptual decision making as the system can switch back to the original decision state from intermediate state if the accumulated evidences are not enough for change of mind.

The stimuli strength applied to the two selective pools is a key factor in decision making tasks. In Figure 4, we quantitatively mapped out the potential landscape for the decision-making dynamics under different stimulus strength μ. The landscape displays three qualitatively different topographies (monostable, bistable and tristable) across the range of stimuli strength from 5 to 90Hz.

For small inputs, the spontaneous state with both selective populations firing at low rates, is the exclusively stable state (Figure 4A). The additional increase of stimulus strength induces the emergence of the two decision attractors (Figure 4B). The central basin of the spontaneous state becomes weaker for increasing stimulus strength (Figure 4C) until the system operates in a binary decision-making region (Figure 4D). When the selective inputs are sufficiently high, an intermediate state with high activity for both neural assemblies encoding the possible alternatives (Figure 4E), called "double-up" state, emerges, in line with both experimental evidences and computational results in delay response tasks (Shadlen and Newsome, 1996; Roitman and Shadlen, 2002; Huk and Shadlen, 2005; Wong and Wang, 2006; Martí et al., 2008; Albantakis and Deco, 2011; Yan et al., 2016; Yan and Wang, 2020). The “double-up” state reduces the barrier between decision states and the transitions are more likely to occur, facilitating changes of mind (Albantakis and Deco, 2011; Yan and Wang, 2020). The two symmetric decision states disappears simultaneously for the further increase of the stimulus strength, so the system loses its ability to compute a categorical choice (Figure 4F).

To track the landscape change and dynamical transition property of decision-making process, we calculated the typical kinetic transition paths averaged over many trials for continuously increasing (from 5 to 90Hz, magenta lines) and decreasing (from 90 to 5Hz, green lines) stimulus inputs, as displayed in Figures 5A–D. There exists four combinations and corresponding landscapes of the decision path since the forward and backward paths can pass through either the DS1 or the DS2 state due to the stochasticity. In Figures 5A,D, the backward path deviates away from the forward decision-making path due to the non-equilibrium flux force although the two paths pass through the same decision state (DS1 or DS2). Figures 5B,C show the forward and backward paths pass through distinct decision states. It can be seen that the system will spend some time staying at the spontaneous state (μ = 5, 16) before making a decision (μ = 40) for the forward paths while maintaining the "double-up" state for a period of time (μ = 90, 76) before making a decision (μ = 40) for the backward paths, implying the existence of hysteresis for state switches in biological neural system.

The energy landscape we obtained here also explains the delayed version of the RDM task well, which involves both the decision computation and working memory. The delayed RDM task additionally requires the subjects to withhold the choice in working memory across a delay period before responding to the saccadic eye movement (Shadlen and Newsome, 2001). Figure 5E shows the three layers of landscape corresponding to the targets, motion and delay in the delayed RDM task, respectively. The gray balls represent the instantaneous states of the system. The megenta and green lines are two possible decision paths due to the symmetry and fluctuations. Initially, the system rests at the spontaneous state until the visual inputs force the system switching into one of the two decision attractor states. The decision is stored in the DS1 or DS2 and this information will be retrieved to produce motor responses at the end of the delay period even though the stimulus is absent. This quantitative picture of landscape echos with the one-dimensional illustrative diagram of decision “landscape” proposed in Wong and Wang (2006).

To provide a more global picture of the system dynamics, we explored the change of landscape topography when the stimulus strength μ and recurrent connectivity in the two selective populations w+ vary (Figures 6A–I). We found that the number of stable states increases as the recurrent connectivity w+ increases (horizontal direction), while the stimulus strength μ increases (vertical direction) there is no apparent trend for the occurrence of multistability. This suggests that stronger self-activation strength may promote the occurrence of multistability.

The stability of attractor states is critical for decision-making network since the decisions can be changed easily for unstable decision state. Experimental and modeling works both suggested that the probability of changes of mind depends on the task difficulty, namely, the motion strength c (Resulaj et al., 2009; Albantakis and Deco, 2011). For biased stimulus input (c≠0), the symmetry of attractor landscape is broken (Supplementary Figures 5, 6 for the asymmetric landscapes). Correspondingly, we can calculate relative barrier height (RB) between pairs of local minima and mean first passage time (MFPT) to quantitatively measure the global stability of the neural network under different motion strength c.

For bistable system, the relative barrier height (RB₁) is defined in terms of the two basins of decision state. Based on landscape topography, we define U₁ and U₂ as the potential minimum of the DS1 and DS2, and U_saddle as the potential of the saddle point between the two states. The barrier heights quantifying the global stability are BH₁ = U_saddle − U₁ and BH₂ = U_saddle − U₂. Then RB₁ = BH₁ − BH₂, quantifying the relative stability of DS1 against DS2.

The MFPT represents the average kinetic time of the alternations between different attractor states, thus also describing the stability of attractors. With dynamical neural activity trajectories over a long time window, we can estimate the time taken for the system to switch from one attractor to another for the first time and get the first passage time (FPT). Here, an attractor region is roughly taken as a small ellipse centered on the local minimum (megenta and green dashed circles in Figure 2C), and the transition is finished once the stochastic trajectory enters the destination ellipse. Then the MFPT is defined as the average of FPT by samplings.

We study whether the landscape can quantify the transition time. Figure 7A shows that for bistable system, RB₁ increases as we increase c for c = 0.5%, 1%, 1.6%, 2.4%, 3.2%. The positive value of RB₁ indicates that the two stable states are no longer symmetric and the DS1 has a deeper basin of attraction than DS2 due to the biased external inputs. When the motion strength gets stronger (larger c), the time taken to switch from DS1 to DS2 increases with higher barrier (Figure 7B), indicating that when the task is easier (larger c) the change of mind needs more time, i.e., the change of mind (the right decision becomes wrong decision in current case) is harder.

We further ask how the model predictions from landscape fit quantitative experimental data. We use tristable system as an example to describe the decision process. Here, the relative barrier height (RB_SS) between SS and DS1 has the similar definition as RB₁. Figure 7C shows that for tristable system, RB_SS increases as we decrease c for c = 12.8%, 6.4%, 3.2%, 0%, indicating that it is harder to make a decision since the barrier is higher for the more difficult task. The transition from SS to DS1 represents the decision process, thus related to the behavioral reaction time (RT). The positive correlation between RB_SS and the experimentally observed RT (Roitman and Shadlen, 2002) suggests that the speed of decisions decreases in harder decision tasks (smaller c) with higher barrier (Figure 7D). Intuitively, this means that for a more difficult task, it is harder to make decision, which takes longer reaction time. Therefore, the landscape results explain the experimental data well and provide an intuitive and quantitative way to understand the decision making function.

Of note, the linear correlation between the logarithm of MFPT and RT and the barrier height of landscape is not a perfect fit. In fact, this is due to another important factor, the curl flux. For a non-equilibrium system, the landscape may explain a major part of the transition rate (barrier crossing from a potential valley). However, the landscape does not solely determine the transition rate as in an equilibrium system, i.e., the flux will contribute to the transition dynamics for barrier crossing (Feng et al., 2014). We should notice that due to the approximation nature of this approach, at very small or very large time scales, it might be less precise for the estimation of the relationship between MFPT and barrier height.

The biologically plausible model for binary decision making was first introduced in Wang (2002), illustrated in Figure 1. The model is a fully connected network, composed of N neurons with N_E = 1600 for excitatory pyramidal cells and N_I = 400 for inhibitory interneurons, which is consistent with the observed proportions of the pyramidal neurons and interneurons in the cerebral cortex (Abeles, 1991). Two distinct populations in excitatory neurons (S1 and S2) respond to the two visual stimuli, respectively, with N_S1 = N_S2 = fN_E and f = 0.15, and the remaining (1−2f)N_E non-selective neurons (NS) do not respond to either of the stimuli. All neurons receive a background Poisson input of v_ext = 2.4kHz. This can be viewed as originating from 800 excitatory connections from external neurons firing at 3Hz per neuron, which is consistent with the resting activity observed in the cerebral cortex (Rolls et al., 1998).

Table 1 is the connection weight matrix between these four populations, where minus sign represents inhibition effect. The synaptic weights between neurons are prescribed according to the Hebbian rule and remain fixed during the simulation. Therefore, inside the selective populations, the synapses are potentiated and the weight w+ is larger than 1. Between selective populations and from nonselective population to selective ones, the synaptic weight is w₋ = 1 − f(w₊ − 1)/(1 − f) so that the overall recurrent excitatory synaptic drive in the spontaneous state remains constant when altering w+ (Amit and Brunel, 1997). w₋ < 1 indicates the synaptic repression. The remaining weights are 1.

Herein, we consider the leakage integrate-and-fire (LIF) model to describe both pyramidal cells and interneurons (Tuckwell, 1988). The membrane potential of a neuron, V(t), can be described by a capacitance-voltage (CV) equation when it is less than a given voltage threshold V_th,

Here, C_m is the capacitance of the neuron membrane, C_m = 0.5nF (0.2nF) for excitatory (inhibitory) neurons. g_L is the leakage conductance, g_L = 25nS (20nS) for excitatory (inhibitory) neurons. Each neuron has a leakage voltage V_L = −70mV and firing threshold V_th = −50mV. I_syn represents the total synaptic current flowing into the neuron.

When V(t) = V_th at t = t_k, the neuron will emit a spike and the membrane potential is reset at V_reset = −55mV for a refractory period τ_ref,

After then, V(t) is governed by the CV (Equation 1) again. Here, τ_ref = 2ms (1ms) for excitatory (inhibitory) neurons.

The synaptic model maps the spike trains of the presynaptic neuron to the postsynaptic current. For the fully connected neural network, the total postsynaptic current is the sum of the following four currents:

The first term is the external excitatory current, which is assumed to be exclusively mediated by AMPA receptor. The second and third terms are the recurrent excitatory currents mediated by AMPA and NMDA receptors. The last term is the inhibitory current mediated by GABA receptor. More specifically,

where V_E = 0mV, V_I = −70mV. For excitatory cells, the synaptic conductances for different channels are g_{ext, AMPA} = 2.1nS, g_{rec, AMPA} = 0.05nS, g_NMDA = 0.165nS, and g_GABA = 1.3nS; for inhibitory cells, g_{ext, AMPA} = 1.62nS, g_{rec, AMPA} = 0.04nS, g_NMDA = 0.13nS, and g_GABA = 1.0nS. w_j is the dimensionless synaptic weight and s_j is the gating variable, representing the fraction of open channels for different receptors. The sum over j denotes a sum over the synapses formed by presynaptic neurons j. Specially, the NMDA synaptic currents depend on both the membrane potential and the extracellular magnesium concentration ([Mg²⁺] = 1mM) (Jahr and Stevens, 1990). The gating variables are given by

where the decay time for AMPA, NMDA and GABA synapses are τ_AMPA = 2ms, τ_{NMDA, decay} = 100ms and τ_GABA = 5ms (Hestrin et al., 1990; Spruston et al., 1995; Salin and Prince, 1996; Xiang et al., 1998). The rise time for NMDA synapses are τ_{NMDA, rise} = 2ms (the rise times for AMPA and GABA are neglected because they are typically very short) and α = 0.5ms⁻¹ (Hestrin et al., 1990; Spruston et al., 1995). The sum over k represents the sum over all the spikes emitted by presynaptic neuron j at time tjk. For external AMPA currents, the spikes are generated by independently sampling Poisson process with rate v_ext = 2.4kHz from cell to cell.

By numerically simulating the neural network system for a long time, we can obtain the raster plot, i.e., neurons emit spikes at specific time points. However, since the spiking neural network we studied here is a high dimensional system consisting of thousands of interacting neurons, we focus on the macroscopic activity of population-averaged firing rate rather than microscopic neural spikes, which effectively reduced the dimension of the system dynamics to 4 (4 populations). For visualization, two more important dimensions, i.e., the population-averaged firing rates of two selective neural groups r1 and r2 are chosen as coordinates to form the "decision space" and then the probablistic distribution at steady state (energy landscape) is projected into the decision space by integrating other dimensions. r₁ and r₂ can be calculated by firstly counting the total spike numbers of a population in a time window of 50ms, which slides with a time step of 5ms, and then dividing it by the neuron number and the time window to get r₁ and r₂.

To visualize the probability distribution of the system state in the "decision space" constructed by the firing rate r₁ and r₂, we discrete the space into a collection of grids and collect the statistics for the system state falling into each grid. Finally, the potential landscape is mapped out by U(r₁, r₂) = −lnP_ss(r₁, r₂) (Sasai and Wolynes, 2003; Li and Wang, 2013a, 2014a), where P_ss(r₁, r₂) represents the normalized joint probability distribution at steady state and U(r₁, r₂) is the dimensionless potential.

Of note, a key issue is to decide when a stationary distribution has been reached. Theoretically, the steady state distribution need to be obtained as time t goes to infinity (or very large). Since the numerical simulation of spiking neural network is computationally intensive and time-consuming, the time length of the dynamical trajectory we can obtained is limited. To address this problem, we define relative Euclidean distance between two probability distributions as σ=∑ij(Pijt+500-Pijt)2(∑ijPt)2, where Pijt and Pijt+500 are probability distributions obtained by firing rate activity with time length t and t+500, respectively. So σ measures the deviation of distribution by prolonging the trajectory by 500s. If increasing time t does not significantly change this relative distance for the probability distribution between different time length (σ is less than a threshold, σ < 0.06% for bistable landscape in Figure 2), we consider that a steady state has been reached.

Biological systems, including neural circuits, are generally dissipative, exchanging energies or materials with the environment to perform functions (Lan et al., 2012). For a non-equilibrium system, the violation of detailed balance lies at the heart of its dynamics. Different from the equilibrium system whose dynamics is solely determined by the underlying energy landscape, the non-equilibrium system is also driven by the steady state probabilistic flux, which measures to what extent the system is out of equilibrium or the detailed balance is broken (Wang et al., 2008; Yan et al., 2013, 2016; Li and Wang, 2014a; Yan and Wang, 2020). For high-dimensional systems, for example, the neural network we studied here, it is challenging to quantify the non-equilibrium probabilistic flux from diffusion equation (Wang et al., 2008; Li and Wang, 2014a). Therefore, we employ an approach which is based on the fluctuating steady-state trajectories, to quantify the probability flux (Battle et al., 2016). As we discussed before, the reduced stochastic system trajectory evolves over time in a four-dimensional phase space. However, for the sake of simplicity and visualization, the two more important dimensions, r₁ and r₂, are utilized to estimate the probabilistic flux. To determine the probability flux of the non-equilibrium neural system, we discrete the subspace constructed by r₁ and r₂ in a coarse-grained way, i.e., the subspace is divided into N₁ × N₂ equally sized, rectangular boxes, each of which represents a discrete state α. Such a discrete state is a continuous set of microstates. Then the probability flux associated with state α is the following vector:

Here, x→α is the center position of the box related to state α. There exist four possible transitions since state α has two neighboring states in each direction. The rate wα-,α(r1) is the net rate of transitions into state α from the adjacent state α⁻ (i.e., the state with smaller r₁), while wα,α+(r1) represents the rate of transitions from α to α⁺ (the state with larger r₁). Similarly, wα-,α(r2) and wα,α+(r2) denote corresponding transitions rates between boxes arranged along the r₂ direction, respectively. Of note, each rate has a sign. For example, wα-,α(r1)<0 means there are more transitions from α to α⁻ than the reverse direction per unit time.

These rates can be estimated by the temporal trajectories of r₁ and r₂:

where t_total is the total simulation time and Nα,β(ri) (Nβ,α(ri)) is the number of transitions from state α (β) to state β (α) along the direction r_i. For the cases where the trajectories go from one box to a non-adjacent box in a single time-step, we perform a linear interpolation to capture all the transitions between adjacent boxes.

To quantify the transition paths between the steady-state attractor states on landscape, we firstly treat a small spherical area centered on the local minimum point as a attractor area (green and megenta dashed circles in Figures 2C, 3C), denoted by B₁, ⋯ , B_m where m is the total number of attractors. The transition path from attractor B_i to attractor B_j is defined as a set of trajectory points {X_t1, ⋯ , X_tn} starting from B_i and ending at B_j. Since the trajectory is noisy, the start point is defined as the last point in B_i before leaving B_i and the transition is finished once the trajectory wanders into the area of attractor B_j. Each transition path corresponds to a non-linear mapping φk:[0,tend]→ℝd. Suppose there are K transition paths from B_i to B_j, then the average transition path ψ_ij is defined as

Of note, since the noise should be large enough to drive the transitions between attractors, the time durations of each single transition trajectory is diverse. To deal with this problem, we firstly split the trajectories evenly into equal number of points and then average these points to form the average trajectory.