To investigate how excitatory inputs trigger DA neuron populations, a bidirectionally connected SW network is developed by randomly rewiring the connections of a ring lattice with a probability of p = 0.1. The construction process includes the following steps: (1) Start with a ring of N = 50 DA neurons. (2) Connect each neuron to its 4 nearest neighbors (Miller et al., 2021). (3) Reconnect each connection to a randomly selected neuron with a probability of p. And our model of networked DA neurons (Equation 1) is as follow:

with i,j = 1,2,…,N representing the DA cell indices, and X_i denoting the state vector of the ith DA neuron. Xi.=F⁢(Xi,IiN⁢M⁢D⁢A,IiC⁢A⁢N) describes the local dynamics of the ith DA neuron, with IiN⁢M⁢D⁢A and IiC⁢A⁢N representing the diverse excitatory inputs to DA neuron i. The coupling relationship among DA neurons is captured by the adjacency matrix A = {a_ij}, with a_ij = a_ji = 1 if DA neurons i and j are functionally interconnected, and a_ij = 0 otherwise. Meanwhile, the diagonal elements of A are set as a_ii = 0.

When DA neuron j is functionally interconnected with neuron i and becomes activated, it releases DA, which binds to D2 receptors on DA neuron i. These D2 receptors regulate GIRK channels through a biochemical cascade involving G proteins, resulting in a delayed and slow D2-GIRK current. The coupling function H(X_j) is described by

with V_j denoting the membrane potential of DA neuron j, τ = 50 ms denoting the lag constant (Beckstead et al., 2004), and Θ_s = −20mV representing the coupling threshold (Otomo et al., 2020). And the coupling current Ii⁢jG⁢I⁢R⁢K is activated by switching H(V_j(t−τ)) from 0 to 1, as described by Tian et al. (2022):

The dynamics of the gating variable r_i (Equation 2) are as follows:

with the functional equations for r_i given by

In line with previous studies (Leone et al., 2015), the coupling current Ii⁢jG⁢I⁢R⁢K is constrained to ensure a consistent total coupling current to each DA neuron: Ii⁢jG⁢I⁢R⁢K=IiG⁢I⁢R⁢K/κi⁢n (with κ_in representing the in-degree of DA neuron i). This constraint reflects biophysical reality, as experimental evidence shows that homeostatic plasticity mechanisms preserve network stability by preventing excessive or insufficient connections in DA neurons (Friedman et al., 2014; Zhang et al., 2019). In this study, the maximal conductance of the total coupling current received by each DA neuron is set at g¯iG⁢I⁢R⁢K=0.005⁢m⁢S/c⁢m2 (Courtney and Ford, 2014).

The local dynamics of DA neurons (Equations 3–9) are derived from previous research (Chen et al., 2022) as follows:

where, for the ith DA neuron, V_i is the membrane potential, h_i is the inactivation variable for Na⁺ current (IiN⁢a); n_i is the activation variable for delayed-rectifier K⁺ current (IiD⁢R); dl_i and fl_i are the activation variable and inactivation variable for L-type Ca²⁺ current (IiC⁢a⁢L), respectively; z_i is the activation variable for SK current (IiS⁢K); Ca_i is the intracellular Ca²⁺ concentration. Iil is the leak current, IiN⁢a⁢P is the persistent Na⁺ current, IiK is the generic persistent K⁺ current, IiN⁢M⁢D⁢A is the NMDA current, and IiC⁢A⁢N is the muscarinic receptor-modulated CAN current.

All currents are represented by chord-conductance equations, as shown below:

And functional equations for all gating variables are

All parameter values are given in Table 1. The values of g¯C⁢A⁢N and g¯N⁢M⁢D⁢A are uniformly distributed over the intervals [0.9, 2.5] mS/cm² and [0, 0.025] mS/cm², respectively.

To evaluate burst synchronization, the order parameter R is utilized as follows (Sun and Xue, 2018):

with N = 50 representing the total number of networked DA neurons, and φ_j(t) denoting the burst phase of the jth DA neuron at the time t, this can be expressed as

where T_j,k represents the onset time of the kth burst of DA neuron j. A higher R value signifies greater burst synchronization. Specifically, R values below 0.4 indicate asynchronization, between 0.4 and 0.8 indicate moderate synchronization, between 0.8 and 0.99 indicate near synchronization, and between 0.99 and 1 indicate full synchronization (Xu et al., 2021).

Each decoupled DA neuron operates on a slow-fast system, where the (V, h, n, dl, fl)-equations constitute a 5-D fast “spiking” subsystem, and the (Ca, z)-equations make up a 2-D slow “bursting” subsystem. The CAN current, controlled by Ca, and the SK current, gated by z, together produce a slow voltage oscillation. When this oscillation exceeds the firing threshold (V_th) of the fast subsystem, a series of action potentials is initiated, namely a burst.

To explore the mechanism of synchronous burst induction, the decoupled DA neuron is reduced into a 3-D model comprising (V, Ca, z) by setting the fast gating variables (h, n, fl, dl) to their steady states. The resulting reduced model (Equations 10–12) is:

Phase-plane analysis of each burst is conducted within the reduced model. The equilibrium bifurcation structure of bursts is defined by the intersection of Ca nullclines and V nullclines at V_th. A tangential intersection indicates a SN bifurcation, characteristic of integrator bursts. A direct intersection, however, signifies an AH bifurcation, indicative of resonator bursts.

To explore how diverse excitatory inputs work together to induce resonator bursts, a fast-slow analysis of electrical activities within the full DA model is conducted, with special attention paid to the differences in state-space distributions of Ca and z dynamics during the initiation of different bursts.

To examine the stability of synchronization, how the system’s synchronization behaves is analyzed when subjected to finite perturbations. The specific idea is: once DA neuronal network has achieved a stable synchronization, small perturbations will be randomly introduced into a subset of DA cells, and then observe whether the system can recover and re-establish synchronization (Arenas et al., 2008).

In the two-node DA network, the occurrence time of the first spike within each burst is used to represent the burst occurrence time, denoted as t1(n) and t2(n) respectively. Then, the phase difference between bursts of coupled DA cells is calculated using the formula:

where T_burst represents the bursting period. The variation process of Ψ_n with respect to the number of bursting periods n after finite perturbations is used to characterize synchronization stability. If the Ψ_n−n curve decays monotonically or oscillatorily to the pre-perturbation level, this will intuitively confirm the stability of synchronization.

In complex networks composed of 50 DA neurons, the decay and recovery processes of the order parameter R(t) following finite perturbations are directly used to characterize the stability of synchronization. By comparing the relaxation time and recovery rate of R(t), the faster the recovery, the more stable the synchronization.

Numerical simulations are executed using Microsoft Visual C++, phase-plane analysis is performed with MATLAB, and fast-slow analysis is carried out using MatCont. Initial values for each DA neuron are randomly assigned, and differential equations are solved using the fourth-order Runge-Kutta method with a step size of 0.01 ms.

To investigate the impact of NMDA receptors on DA network’s behavior, muscarinic receptor-modulated CAN maximal conductance g¯C⁢A⁢N of each DA neuron is maintained at 0.9mS/cm², representing the state before muscarinic receptor activation. The spatiotemporal dynamics of the network are observed by adjusting NMDA maximal conductance g¯N⁢M⁢D⁢A within each DA neuron. With g¯N⁢M⁢D⁢A set to 0mS/cm², reflecting the state before NMDA receptor activation, the network dynamics after a 2.5 sec transient period are illustrated in Figure 1A, with all neurons at resting membrane potentials, indicating that DA network is in a resting state. Increasing g¯N⁢M⁢D⁢A to 0.015 mS/cm² to simulate NMDA receptor activation, as depicted in Figure 1B, results in burst firings in all DA neurons, but these bursts seem asynchronous.

To illustrate how the network dynamics vary with g¯N⁢M⁢D⁢A before muscarinic receptor activation, the order parameter R is plotted against g¯N⁢M⁢D⁢A in Figure 1C. As illustrated, increasing g¯N⁢M⁢D⁢A leads to a gradual rise in R, yet it never exceeds the synchronization threshold of 0.4. This indicates that DA network remains in a state of burst asynchronization throughout the g¯N⁢M⁢D⁢A range, demonstrating that activating NMDA receptors alone is insufficient to induce synchronous bursts.

Previous research showed that, in addition to NMDA receptors, muscarinic receptors are essential for natural reward stimuli to elicit phasic release (Wickham et al., 2015). This finding highlights the need to investigate the impact of muscarinic receptors on NMDA receptor-induced bursts. As depicted in Figure 1D, increasing the CAN maximal conductance g¯C⁢A⁢N to 1.9mS/cm² to mimic muscarinic receptor activation results in the re-emergence of asynchronous bursts prior to NMDA receptor activation (g¯N⁢M⁢D⁢A=0⁢m⁢S/c⁢m2). However, upon NMDA receptor activation (g¯N⁢M⁢D⁢A=0.015⁢m⁢S/c⁢m2), as depicted in Figure 1E, the bursts become more synchronized.

Following muscarinic receptor activation, the R-g¯N⁢M⁢D⁢A curve is explored again. As depicted in Figure 1F, R rises from below 0.4 to between 0.4 and 0.8 as g¯N⁢M⁢D⁢A increases. As per Xu et al. (2021), the g¯N⁢M⁢D⁢A value corresponding to the synchronization threshold R of 0.4 is designated as g¯N⁢M⁢D⁢A,0 (Figure 1F, magenta dot). For g¯N⁢M⁢D⁢A values smaller than g¯N⁢M⁢D⁢A,0, R remains below 0.4, indicating the DA network is in a state of burst asynchronization. However, for g¯N⁢M⁢D⁢A values greater than g¯N⁢M⁢D⁢A,0, R increases from below 0.4 to between 0.4 and 0.8, suggesting the DA network shifts to a state of moderate burst synchronization. These results highlight that muscarinic receptor-modulated CAN current promotes NMDA receptor-induced bursts synchronizing.

To globally understand the impact of muscarinic receptors on NMDA receptor-induced bursts, a detailed diagram of R distribution across a specific region in the two-dimensional parameter space (g¯C⁢A⁢N, g¯N⁢M⁢D⁢A) is scanned. As shown in Figure 2, when the CAN maximal conductance is low (g¯C⁢A⁢N<g¯C⁢A⁢N,c≈1.31), g¯N⁢M⁢D⁢A,0 is nonexistent, indicating the burst asynchronization region occupying the entire g¯N⁢M⁢D⁢A range. And when the CAN maximal conductance is high (g¯C⁢A⁢N>g¯C⁢A⁢N,c), g¯N⁢M⁢D⁢A,0 gradually decreases as g¯C⁢A⁢N increases (Figure 2, magenta dashed line). Therefore, the moderate burst synchronization region expands gradually to the left, occupying a larger portion of the g¯N⁢M⁢D⁢A range. This indicates that higher levels of muscarinic receptor-modulated CAN current significantly enhance the possibility of NMDA receptor-induced bursts synchronizing, indicating a complex interaction between two receptors in promoting synchronous bursts.

Similar R distribution diagrams are observed for other network structures. Figure 3A illustrates the R distribution for the Barabási-Albert (BA) network, which is constructed through the preferential attachment of newly added nodes to those with rich connections. The BA network’s size and number of links match those of the SW network depicted in Figure 2, but its degree distribution follows a power-law nature. As illustrated in Figure 3A, similar to Figure 2, in the high region of the CAN maximal conductance (g¯C⁢A⁢N>g¯C⁢A⁢N,c≈1.38), increasing g¯C⁢A⁢N progressively expands the moderate burst synchronization region to the left, covering a larger portion of the g¯N⁢M⁢D⁢A range; while in the low region of the CAN maximal conductance (g¯C⁢A⁢N<g¯C⁢A⁢N,c), the burst asynchronization region occupies the entire g¯N⁢M⁢D⁢A range. Figure 3B shows the R distribution diagram for the Erdös-Rényi (ER) random network, which is constructed by connecting the existing DA neurons with an equal probability. Again, it is seen that in the high region of the CAN maximal conductance (g¯C⁢A⁢N>g¯C⁢A⁢N,c≈1.42), the moderate burst synchronization region gradually expands to the left as g¯C⁢A⁢N increases, covering a larger portion of the g¯N⁢M⁢D⁢A range; while in the low region of the CAN maximal conductance (g¯C⁢A⁢N<g¯C⁢A⁢N,c), the burst asynchronization region occupies the entire g¯N⁢M⁢D⁢A range. Figures 2, 3 suggest that the main features of the R distribution diagram exhibit qualitatively similar trends across the tested topologies, suggesting the potential existence of a universal mechanism to explain all observed phenomena.

During the induction of synchronous bursts, the network’s structure and connections stay constant; only the external excitatory inputs affecting the DA cells change. How does muscarinic receptor-modulated CAN current facilitate NMDA receptor-induced bursts synchronizing? To address this question, a simplified network is employed, which consists of two DA neurons and is also functionally interconnected via D2-GIRK currents, as illustrated in Figure 4A. The impacts of NMDA and CAN currents on spatiotemporal dynamics and synchronization index are re-evaluated in this dual-node network. The spatiotemporal dynamics are plotted in Figure 4B. Initially, as depicted in Figure 4Ba, the network remains at a resting state before muscarinic and NMDA receptor activation. Separately activating either receptor type drives the network into a burst asynchronization state (Figures 4Bb, Bc). However, simultaneous activation of both receptors leads the network into a burst synchronization state (Figure 4Bd). These results replicate the phenomena that muscarinic receptor-modulated CAN current promotes NMDA receptor-induced bursts synchronizing.

The average minimum time difference between bursts of coupled DA neurons, denoted as Δt, is utilized to assess the synchronization of bursts in this dual-node network. A smaller Δt value signifies greater network synchronization. As depicted in Figure 4C, the Δt distribution diagram reproduces the main feature shown in Figures 2, 3. To be specific, in the high region of the CAN maximal conductance (g¯C⁢A⁢N>g¯C⁢A⁢N,c≈0.96), as g¯C⁢A⁢N increases, the burst synchronization region (Δt < 60ms) progressively expands to the left, occupying a larger portion of the g¯N⁢M⁢D⁢A range; while in the low region of the CAN maximal conductance (g¯C⁢A⁢N<g¯C⁢A⁢N,c), the burst asynchronization region (Δt > 60ms) occupies the entire g¯N⁢M⁢D⁢A range. These findings reaffirm that higher levels of muscarinic receptor-modulated CAN current significantly increase the possibility of NMDA receptor-induced bursts synchronizing. The striking similarities between Figures 2–4 render that the dual-node network is ideal for exploring the mechanism by which diverse excitatory inputs synergistically induce synchronous bursts.

The effect of connections between DA neurons on these synchronous bursts is then analyzed. In order to monitor the changes in the couplings between DA neurons under different network dynamics, due to the symmetry of the dual-node network structure, the D2-GIRK current affecting neuron 1 (I1G⁢I⁢R⁢K) is selected for observation. Assuming D2 receptors produce sufficient G protein, the gating variable r₁ of I1G⁢I⁢R⁢K depends solely on neuron 1’s membrane potential V₁. As depicted in Figure 4D, the steady activation curve for r₁ decreases gradually with increasing V₁. The analyses of I1G⁢I⁢R⁢K for different network dynamics are as follows. Prior to the activation of NMDA and muscarinic receptors, both DA neurons are in their resting state, with V₁ and V₂ at approximately −40mV, below the coupling threshold (Θ_S = −20mV). Therefore, H(V₂) is 0, resulting in I1G⁢I⁢R⁢K being 0pA, despite a relatively high r₁ (Figure 4Ea). After activation of NMDA or muscarinic receptors, burst firings occur, rendering V₂ to exceed the coupling threshold, thus H(V₂) equals 1. During this period, neuron 1 is in its inter-burst phase, with V₁ typically below −40mV, resulting in a relatively large r₁ and thus a high amplitude of I1G⁢I⁢R⁢K (Figures 4Eb, 4Ec). Simultaneous activation of NMDA and muscarinic receptors also triggers burst firings in neuron 2, making H(V₂) equal 1. However, neuron 1 is in its intra-burst period, with V₁ often above 0mV, leading to a comparatively low r₁ and thus a marginally reduced amplitude of I1G⁢I⁢R⁢K (Figure 4Ed).

Theoretical studies on spiking networks have shown that sufficient inhibitory couplings, i.e., increases in inhibitory couplings, can improve the synchronization of resonators (Izhikevich, 2010; Chik et al., 2004). Our findings reveal that inhibitory coupling regulates the dynamics of bursting DA networks; However, compared to asynchronous states, its amplitude diminishes during synchronization, indicating that it is not the primary driver of synchronization. This prompts critical questions: If synchronization does not arise from inhibitory couplings, could it instead stem from the local dynamics of nodes? More specifically, might excitatory inputs drive a transition in burst modes—from integrators to resonators?

The analysis then explores the properties of decoupled DA neurons to determine if they function as integrators or resonators. This is achieved by studying the bifurcation structure of equilibria during burst initiation, which involves two steps: identifying the discharge threshold V_th of the fast subsystem and conducting a phase-plane analysis within the reduced DA model.

To identify the discharge threshold V_th of the fast subsystem, the slow subsystem is removed away from the full DA model. V_th is then calculated within the remaining fast subsystem by adjusting the external DC current I₀. The results reveal that V_th consistently stays around −38.7mV, regardless of variations in g¯N⁢M⁢D⁢A.

In the reduced DA model, conduct phase-plane analysis for each burst by setting z as the control parameter and calculating Ca nullclines and V nullclines. As shown in Figure 5, when z is large, the V nullclines (Figures 5A–C, blue curves) intersect the Ca nullclines (Figures 5A–C, black curves) at different voltages on both sides of V_th (Figures 5A–C, black dots). These intersections represent equilibria, with the right intersection being unstable and the left one stable. Therefore, the system resides at the left equilibrium, representing the resting state. A significant decrease in z results in a notable downward shift of the V nullclines (Figures 5A–C, green curves), preventing them from intersecting with the Ca nullclines in the region of interest. This change indicates the system has entered the firing state.

A moderate decrease in z leads to the intersection of Ca nullclines and V nullclines shifting to V_th, marking the critical transition from a resting state to a firing state. This intersection at V_th varies across different scenarios. When NMDA and muscarinic receptors are activated separately, the Ca nullclines are tangent to the V nullclines (Figures 5A, B, red curves). However, simultaneous activation of both receptors results in the Ca nullclines intersecting the V nullclines (Figure 5C, red curve).

At V_th, the tangency between the Ca nullclines and V nullclines signifies a SN bifurcation structure, confirming that bursts function as integrators. However, at V_th, their intersection denotes an AH bifurcation, demonstrating that bursts act as resonators. This reveals that separately activating two receptors generates challenging-to-synchronize integrator bursts, while their simultaneous activation yields easy-to-synchronize resonator bursts. In summary, these findings illustrate that diverse excitatory inputs collaborate to induce synchronous bursts by transitioning burst behavior from integrator to resonator mode.

To investigate how diverse excitatory inputs synergistically transition burst behavior to resonator mode, a fast-slow analysis within the full DA model is performed. For each electrical activity, the bifurcation curve of equilibria and the state-space distribution of Ca and z dynamics are plotted on the Ca-z plane. As shown in Figure 6, as Ca and z increase, the bifurcation curve of equilibria evolves from the SN curve (Figures 6A–D, blue curves) to the AH curve (Figures 6A–D, green curves) via a BT bifurcation (Figures 6A–D, magenta dots). These curves and bifurcations partition each Ca-z plane into a firing region (Figures 6A–D, shaded area) and a quiescent region (Figures 6A–D, unshaded area).

During examining the state-space distributions of Ca and z dynamics across different electrical activities, two key differences are identified. Firstly, there is a notable difference in form. Before activating NMDA and muscarinic receptors, the state-space distribution stabilizes at a point within the quiescent region of the Ca-z plane (Figure 6A, black dot). After activation, the distributions transition into trajectories oscillating between the quiescent and firing regions on the Ca-z plane (Figures 6B–D, black circles). A closer examination on each trajectory reveals that the burst terminates when Ca²⁺ accumulates enough to activate z of I_SK, leading to the trajectory to intersect the AH curve at a higher Ca_AH, thus entering the quiescent region and halting Ca²⁺ influx. The burst initiates when Ca drops sufficiently due to Ca²⁺ clearance, deactivating z, and resulting in the trajectory crossing the bifurcation curve at a lower Ca, thereby exiting the quiescent region and starting Ca²⁺ accumulation. Secondly, there is a vital difference in the bifurcation curves intersected by trajectories during burst initiation. When bursts are induced by the individual activation of two receptors, their trajectories intersect the SN curves (Figures 6B, C). However, when bursts are induced by the simultaneous activation of two receptors, the trajectory intersects the AH curve (Figure 6D).

To comprehend why trajectories intersect different bifurcation curves during burst initiation, it is essential to elucidate the relationship between Ca and z. As indicated in Equation (8), z is dependent on Ca, simplifying the slow subsystem to Ca alone. Previous studies (Chen et al., 2022) have shown that burst induction relies on a positive feedback loop between Ca and I_CAN. Activating NMDA receptors supports this loop by directly providing Ca²⁺. Activating muscarinic receptors boosts the loop by increasing I_CAN, causing depolarization and activating L-type Ca²⁺ channels, thereby indirectly supplying Ca²⁺. These mechanisms moderately boost Ca dynamics, respectively (Figures 6B, C). When both receptors are activated simultaneously, the combined direct and indirect Ca²⁺ significantly enhance Ca dynamics (Figure 6D). A moderate rise in Ca results in a corresponding moderate increase in z, while a significant rise in Ca leads to a substantial increase in z. Therefore, during burst initiation, a moderate increase in Ca and z dynamics typically causes the trajectory to intersect with the SN curve at a lower Ca_SN, while a significant increase in Ca and z dynamics leads to the trajectory intersecting with the AH curve at a higher Ca_AH.

In conclusion, sufficient Ca²⁺ accumulation leads to transitioning burst behavior from integrator to resonator mode. Simultaneous activation of both receptors is more effective in achieving the required Ca²⁺ accumulation compared to activating either receptor alone.

When all the oscillators within DA neuronal networks are set at the synchronization manifold, they will remain synchronized. Now the crucial question is whether the synchronization manifold is stable in the presence of small perturbations.

To address this question, in the two-node DA network, once the system has stabilized and achieved synchronization, neuron 1 is subjected to a current pulse I₀ with a duration of 25 ms. When I₀ = 2 pA, this represents an excitatory small perturbation (Figure 7Aa); when I₀ = −2 pA, it represents an inhibitory small perturbation (Figure 7Ab). The occurrence times of the bursts, t1(n) and t2(n), are carefully recorded, and the phase difference between bursts of the coupled DA cells, Ψ_n, is computed. As shown in Figure 7B, both excitatory and inhibitory perturbations induce changes in Ψ_n; however, after several bursting periods, these changes decay oscillatively back to the pre-perturbation level (Figure 7B). These results strongly confirm the stability of burst synchronization induced by excitatory inputs.

In complex networks composed of 50 DA neurons, after the systems reach moderate synchronized states, same current pulses, representing perturbations, are randomly introduced into 10 of the cells. The decay and recovery processes of the order parameter R(t) are then carefully monitored and recorded. As shown in Figure 8, overall, in all complex networks, R(t) exhibits a process of initial transient decay followed by gradual recovery. The subtle difference lies in that the relaxation time of R(t) in the SW network is significantly shorter than that in scale-free networks and random networks. These results further confirm that synchronization induced by excitatory inputs is stable, and such synchronization stability is more pronounced in SW networks.