Abstract

Phasic dopamine (DA) release is related to reward processing and addiction. The prevailing view posits that this DA release originates from synchronous bursts of DA neuron groups, rather than individual neuron bursts. However, the mechanism by which diverse excitatory inputs synergistically induce synchronous bursts remains unclear. In biophysically realistic networks with complex structure, the responses of functionally connected DA neurons to various excitatory inputs are examined. Activating NMDA receptors alone results in asynchronous bursts, while co-activating with muscarinic receptors significantly enhances burst synchronization. The synchronization trends display qualitatively similar characteristics across all tested topological networks, indicating that these synchronous bursts are universal. Research on a dual-node network reveals that inhibitory couplings, specifically the inward rectifying K+ currents activated by G protein linked to D2 receptors (D2-GIRK currents), participate in inducing synchronous bursts. A detailed analysis of decoupled DA neurons shows that these synchronous bursts are induced by transitioning bursts from integrator-like to resonator-like behavior, a process dependent on sufficient intracellular Ca2+ accumulation. NMDA receptors directly supply Ca2+, while muscarinic receptors indirectly provide Ca2+ by enhancing calcium-activated, nonspecific cation (CAN) current, causing depolarization, and activating L-type calcium channels. Therefore, simultaneous activation of both receptors is more effective in achieving the required Ca2+ accumulation than activating either receptor alone. These findings elucidate the mechanism by which diverse excitatory inputs work together to induce synchronous bursts, providing new insights into their inductions and regulations, potentially advancing our understanding of the physiological diversity of phasic DA releases and their addictive abnormalities.

1 Introduction

The ventral tegmental area (VTA) DA system is vital for reward-based neural activities, driving adaptive behaviors essential for survival and reproduction (Berke, 2018). Abnormalities in this system have been associated with conditions such as addiction. Traditionally, it was believed that changes in individual DA neuron electrical activity over time encode DA release, with tonic release associated with low-frequency spike firings and phasic release linked to high-frequency burst firings. However, this view is increasingly being challenged (Morales and Margolis, 2017), with experimental data supporting both homogeneity in individual DA neuron burst firing patterns (Eshel et al., 2016) and diversity in phasic releases in target regions (Parker et al., 2016). Several studies have reported that the synchronization of DA neuron populations is also involved in phasic release induction (Joshua et al., 2009; Li et al., 2011). As a result, an emerging perspective posits that spatial organization (synchronicity) of bursts is crucial for encoding robust phasic release (Beeler and Dreyer, 2019; Liu et al., 2021). Nonetheless, the precise mechanisms remain unclear.

Synchronization in neural networks depends on network structure, connections between neurons, and properties of neurons. Considering the impact of network structure on synchronization, DA networks exhibit stable complex network features, such as small-world (SW) properties, and cannot be regulated through D2 receptors (Miller et al., 2021). However, the role of network structure in synchronous bursts induced by excitatory inputs is still not fully understood.

Regarding the impact of connections between DA neurons upon synchronization, direct connections like gap junctions and synaptic links are relatively rare, whereas functional connections are prevalent. Activating D2 receptors by DA generates an inhibitory postsynaptic current, namely the D2-GIRK current (Beckstead et al., 2004). Theoretical studies show that inhibitory couplings typically lead to anti-phase (Li et al., 2019) or clustered synchronous bursts (Kim and Lim, 2019), which are less effective in producing phasic release compared to in-phase synchronous bursts. However, experimental evidence reveals that increased D2 receptors paradoxically facilitate phasic release (Kita et al., 2007; Jones and Fordahl, 2021). Theoretical studies on spiking networks reveal that sufficient inhibitory couplings can foster the synchronization of resonators (Izhikevich, 2010; Chik et al., 2004). This implies that strengthening inhibitory couplings through D2 receptors might similarly promote the synchronization of resonator bursts in DA neuronal networks.

Concerning the last category, akin to spikes, most bursts generally exhibit either a Saddle-Node (SN) or an Andronov-Hopf (AH) bifurcation structure of equilibria (Izhikevich, 2010, 2000), allowing them to function as integrators or resonators. Theoretical studies on spikes show that integrators and resonators respond differently to brief inhibitory perturbations. Integrators usually experience phase delays, whereas resonators often display both phase delays and phase advances (Izhikevich, 2010). These distinct responses result in significantly different synchronization properties, with resonator networks synchronizing more effectively than integrator networks. The properties of DA neurons are regulated by diverse external excitatory inputs through muscarinic and NMDA receptors (Wickham et al., 2015; Galaj et al., 2022). Studies on these bursts (Chen et al., 2022) reveal that, in the parameter plane defined by the intracellular Ca2+ concentration (Ca) and the gating variable z of SK current, the state-space trajectory of Ca and z dynamics intersects the SN bifurcation curve of equilibria during burst initiation. This finding aligns with earlier studies (Yu and Canavier, 2015; Oster et al., 2015), suggesting that these bursts function as integrators. However, as Ca and z increase, the bifurcation curve transitions from the SN to the AH curve via a Bogdanov-Takens (BT) bifurcation (Chen et al., 2022). This implies that the trajectory may intersect the AH curve during burst initiation, potentially allowing excitatory inputs to switch bursts from integrator to resonator behavior, thereby enhancing synchronization of the DA network.

Although activation of muscarine (Zhang et al., 2005) and NMDA receptors (Chergui et al., 1993) can both trigger bursts, co-infusion of carbachol, a muscarinic receptor agonist, with AP5, an NMDA receptor antagonist, into the VTA eliminated phasic release (Spanos et al., 2019), suggesting burst asynchronization; however, co-administration of carbachol and NMDA into the VTA restored phasic release (Spanos et al., 2019), indicating burst synchronization. These results imply that activating muscarinic and NMDA receptors separately tends to induce difficult-to-synchronize integrator bursts, whereas their simultaneous activation results in easy-to-synchronize resonator bursts, suggesting a synergistic effect between these receptors on inducing resonator bursts, though the exact mechanism is yet to be elucidated.

VTA DA neurons express TRPC channels, which mediate the CAN current (Klipec et al., 2016). Modulating the expression and function of these channels helps regulate DA neuron firing, and in turn, influences animal behavior (Wang et al., 2024). As G protein-coupled receptors, muscarinic receptors have the ability to modulate TRPC channel expression and function—a capability that has been confirmed in hippocampal pyramidal neurons (Tai et al., 2010). To substantiate speculations above, a functional network with SW characteristics is created, comprising 50 DA neurons bidirectionally connected by D2-GIRK currents. Each neuron features an ICAN modulated by muscarinic receptors and an INMDA mediated by NMDA receptors (Chen et al., 2022). To examine the responses of DA neuron populations to diverse excitatory inputs, the spatiotemporal dynamics of networks are observed after separately and simultaneously activating NMDA and muscarinic receptors. To assess the impact of diverse excitatory inputs on synchronous bursts, synchronization index curves as a function of NMDA maximal conductance are compared before and after muscarinic receptor activation. To evaluate the influence of network structure on these bursts, synchronization index distribution diagrams across networks with different structures are compared. To investigate the role of inhibitory couplings in synchronous burst induction, D2-GIRK currents are evaluated across different network dynamics within a dual-node DA network. To understand how the local properties of DA neurons affect synchronous burst induction, a phase-plane analysis is performed on bursts of all conditions within a reduced model of decoupled DA neurons. To elucidate how muscarinic and NMDA receptors collaborate to transition bursts from integrator to resonator behavior, detailed fast-slow analyses are conducted on all electrical activities within a full model of decoupled DA neuron. To investigate the stability of the synchronization manifold induced by excitatory inputs, synchronization stability analyses of the systems under finite small perturbations are conducted.

Numerical results demonstrate that activating NMDA receptors alone shifts the DA network from a resting state to a burst asynchronization state. However, activating muscarinic receptors significantly improves the synchronization of NMDA receptor-induced bursts. These synchronous bursts are universal, as they exhibit qualitatively similar synchronization trends across all tested topologies. Further research on a dual-node DA network shows that inhibitory couplings participate in inducing synchronous bursts. Analyses on decoupled DA neurons reveal that transitioning to resonator bursts is crucial for synchronous burst induction, a process dependent on sufficient intracellular Ca2+ accumulation, which is more effectively achieved by simultaneously activating muscarinic and NMDA receptors. Stability analysis conducted under finite perturbations confirms that the synchronization induced by excitatory inputs is stable. These findings elucidate how diverse excitatory inputs collectively trigger synchronized bursts in midbrain DA neurons.

2 Materials and methods2.1 The network model

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 Xi denoting the state vector of the ith DA neuron. describes the local dynamics of the ith DA neuron, with and representing the diverse excitatory inputs to DA neuron i. The coupling relationship among DA neurons is captured by the adjacency matrix A = , with aij = aji = 1 if DA neurons i and j are functionally interconnected, and aij = 0 otherwise. Meanwhile, the diagonal elements of A are set as aii = 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(Xj) is described by

with Vj 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 is activated by switching H(Vj(t−τ)) from 0 to 1, as described by Tian et al. (2022):

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

with the functional equations for ri given by

In line with previous studies (Leone et al., 2015), the coupling current is constrained to ensure a consistent total coupling current to each DA neuron: (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 (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, Vi is the membrane potential, hi is the inactivation variable for Na+ current (); ni is the activation variable for delayed-rectifier K+ current (); dli and fli are the activation variable and inactivation variable for L-type Ca2+ current (), respectively; zi is the activation variable for SK current (); Cai is the intracellular Ca2+ concentration. is the leak current, is the persistent Na+ current, is the generic persistent K+ current, is the NMDA current, and 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 and are uniformly distributed over the intervals [0.9, 2.5] mS/cm2 and [0, 0.025] mS/cm2, respectively.

ParameterValueParameterValueParameterValueParameterValue250.0mS/cm23.5mS/cm2ECAN0mVC1.0μF5.0mS/cm20mS/cm2ENMDA0mVτz100.0ms0.4mS/cm20.005mS/cm2ε0.0025CMg0.5μM0.002mS/cm2ENa55.0mVκ10.3ΘS−20mV0.015mS/cm2EK−90.0mVκ22τ50.0ms0.075mS/cm2ECa100.0mVκSK0.3kin40.9mScm2El−50.0mVCabasal0.005

All parameters are directly adopted from Chen et al. (2022), except for the parameters marked in gray. The reference literature for the gray parameters is detailed in Table 2.

ParameterRangeReferences0.001∼0.01Courtney and Ford, 2014τ50∼60 msBeckstead et al., 2004κin3.97∼5.24Miller et al., 2021ΘS−20 mVOtomo et al., 2020

Source literatures for the gray parameters in Table 1.

2.2 Synchronization index

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 Tj,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).

2.3 Phase-plane analysis and fast-slow analysis

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 (Vth) 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 Vth. 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.

2.4 Stability analysis of synchronization with finite perturbations

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 and respectively. Then, the phase difference between bursts of coupled DA cells is calculated using the formula:

where Tburst 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.

3 Results3.1 Co-activation of NMDA and muscarinic receptors induce synchronous bursts

To investigate the impact of NMDA receptors on DA network’s behavior, muscarinic receptor-modulated CAN maximal conductance of each DA neuron is maintained at 0.9mS/cm2, representing the state before muscarinic receptor activation. The spatiotemporal dynamics of the network are observed by adjusting NMDA maximal conductance within each DA neuron. With set to 0mS/cm2, 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 to 0.015 mS/cm2 to simulate NMDA receptor activation, as depicted in Figure 1B, results in burst firings in all DA neurons, but these bursts seem asynchronous.

Synchronous bursts induced by simultaneous activation of NMDA and muscarinic receptors in a small-world (SW) DA network. Maintaining of each DA neuron at 0.9mS/cm2 to simulate conditions before muscarinic receptor activation, the spatiotemporal dynamics of DA network before () (A) and after () (B) NMDA receptor activation. (C) the variation of order parameter R as a function of before muscarinic receptor activation. Increasing to 1.9mS/cm2 to simulate conditions after muscarinic receptor activation, the network dynamics before () (D) and after () (E) NMDA receptor activation. (F)R- curve after muscarinic receptor activation. The dashed lines denote moderate synchronization threshold of R = 0.4, and the magenta dot denotes the corresponding to moderate synchronization threshold of R = 0.4, designated as .

To illustrate how the network dynamics vary with before muscarinic receptor activation, the order parameter R is plotted against in Figure 1C. As illustrated, increasing 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 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 to 1.9mS/cm2 to mimic muscarinic receptor activation results in the re-emergence of asynchronous bursts prior to NMDA receptor activation (). However, upon NMDA receptor activation (), as depicted in Figure 1E, the bursts become more synchronized.

Following muscarinic receptor activation, the curve is explored again. As depicted in Figure 1F, R rises from below 0.4 to between 0.4 and 0.8 as increases. As per Xu et al. (2021), the value corresponding to the synchronization threshold R of 0.4 is designated as (Figure 1F, magenta dot). For values smaller than , R remains below 0.4, indicating the DA network is in a state of burst asynchronization. However, for values greater than , 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 (, ) is scanned. As shown in Figure 2, when the CAN maximal conductance is low (), is nonexistent, indicating the burst asynchronization region occupying the entire range. And when the CAN maximal conductance is high (), gradually decreases as increases (Figure 2, magenta dashed line). Therefore, the moderate burst synchronization region expands gradually to the left, occupying a larger portion of the 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.

R distribution diagram on the parameter space (, ). When muscarinic receptor-modulated calcium-activated, nonspecific cation (CAN) maximal conductance is high (), with the increase of , monotonically decreases, demonstrating that the moderate burst synchronization region gradually expands to the left, covering a larger portion of the range. When CAN maximal conductance is low (), is nonexistent, indicating the burst asynchronization region occupying the entire range. The horizontal lines denote R- curves shown in Figures 1C, F, respectively. The white dashed line denotes the critical CAN maximal conductance . The magenta dashed curve denotes the boundary composed by .

3.2 Synchronous bursts induced by excitatory inputs are universal

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 (), increasing progressively expands the moderate burst synchronization region to the left, covering a larger portion of the range; while in the low region of the CAN maximal conductance (), the burst asynchronization region occupies the entire 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 (), the moderate burst synchronization region gradually expands to the left as increases, covering a larger portion of the range; while in the low region of the CAN maximal conductance (), the burst asynchronization region occupies the entire 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.

R distribution diagrams for the Barabási-Albert (BA) (A) and Erdös-Rényi (ER) (B) networks. The size and the number of links of the networks are the same to that of SW network used in Figure 2. For the BA (ER) network, the moderate burst synchronization region gradually expands to the left as increases from (), and the burst asynchronization region occupies the entire range in the region . The white dashed lines denote the critical CAN maximal conductance . The magenta dashed curves denote the boundaries composed by .

3.3 Inhibitory couplings participate in inducing synchronous bursts

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.

Synchronous burst induction in a dual-node DA network. (A) a dual-node DA network structure. (B) spatiotemporal dynamics. (C) synchronization index Δt distribution. Inserted dots represent the Δt values corresponding to the network dynamics shown in Figure 4B. The white dashed line denotes the critical CAN maximal conductance . (D) the steady activation curve r∞ for the gating variable r1 with respect to V. (E) the time evolution of r1 (Left) and the coupling current (Right) corresponding to the spatiotemporal dynamics illustrated in Figure 4B.

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,