Biophysical Modeling of Dopaminergic Denervation Landscapes in the Striatum Reveals New Therapeutic Strategy

Functional and spatial characterization of DA signaling in the healthy striatum

We began our investigation by modeling DA signaling in the fully innervated human striatum, specifically the putamen, which we defined as the healthy state (Dreyer et al., 2010; Dreyer, 2014). We simulated the firing of dopaminergic SNc neurons and described the extracellular DA concentration. For this, we employed a model describing DA in a subvolume of 10³ mm³ (see Extended Data 1). Given the estimated density of ∼0.1 dopaminergic axonal terminals per mm³ in the healthy striatum (Doucet et al., 1986; Dreyer et al., 2010; Dreyer, 2014), this volume contains on average 100 terminals, each of which was treated as an individual element. This is a reasonable approximation since each time a neuron fires, only a fraction of its terminals release transmitter. Estimates of the vesicular release probability of dopaminergic terminals are thus within the range of 6–20% (Dreyer, 2014; Pereira et al., 2016), and only ∼30% of the terminals may contain the molecular machinery for exocytosis (Liu et al., 2018). Based on this, we employed a deterministic mean-field model that approximates DA inside the ith subvolume as: dMidt=ΔνNi−VMNiMiKM+Mi + (D∇2Mi−ϵMi).

Here, M is the DA concentration, Δ is the amount of DA released by a terminal, ν is the neuronal firing frequency, V_M is the DA uptake per terminal, and N is the number of terminals within the subvolume. DA remains active in the extracellular space until it is removed by either transporters or degraded enzymatically (Fig. 1A), so we modeled transporter-mediated DA uptake after the Michaelis-Menten uptake equation. We also included a simple degradation term (ϵ), and a term to account for the diffusion between neighboring subvolumes (D). These are, however, so small that they can be neglected if a region has dopaminergic innervation and will only be considered for regions deprived of terminals, and hence they are placed in parentheses. As shown in previous studies, we included that each dopaminergic neuron can express one of three firing patterns: pauses, tonic, or phasic (Grace and Bunney, 1984a, b). Thus, in the model, the firing frequency could take one of three values: 0 Hz (pauses), 4 Hz (for tonic), or 20 Hz (for phasic). All employed parameter values (Extended Data Fig. 1-1) were adopted from previous theoretical work which used experimentally determined measurements to constrain their model (Dreyer et al., 2010; see their Table 1). From this, we obtained DA time courses that clearly reflected the underlying neuronal firing patterns, exhibiting periods of tonic and phasic DA signaling, and pauses where DA is cleared from the extracellular space (Fig. 1C).

We next characterized the dopaminergic innervation of the striatum at a mesoscale. To mimic the shape of the putamen in the human striatum, we modeled it as an ellipsoid. Dopaminergic innervation was constructed by filling the volume with axonal arbors from 10⁵ SNc neurons, based on estimates from human SNc (Hardman et al., 2002) and the fact that these neurons have wide-spread projection targets (Poulin et al., 2018). Each neuron contributed with a spherical arbor with a radius of 0.5 mm (which is our best estimate based on the existing data), wherein the density of terminals was constant (Doucet et al., 1986; Matsuda et al., 2009). To characterize the organization and spatial coverage of arbors within the striatum, we analyzed whether innervation was coherent within the striatum. This can be quantified by the number of what we termed contiguous arbor classes (CACs), inspired by the mathematical analysis of communication classes. We assumed that dopaminergic neurons belonged to the same CAC if their arbors considerably overlapped (their arbor centers <0.5 mm apart; Fig. 1D). Hence, if the number of CACs is low, it suggests a high degree of spatial coverage and cohesion, and vice versa. For classifying neurons into CACs, we used a Markov chain-inspired algorithm (see Extended Data 1). We found that all neurons in the healthy striatum belonged to the same CAC (Fig. 1E), suggesting a high degree of coverage (Fig. 1F). For each neuron, we also counted the number of overlapping arbors, and this metric followed a Poisson distribution (Fig. 1G).

From the equation above, it is evident that decreasing Ni does not affect the steady state DA concentration considerably unless this value is approximately zero. From calculations on the diffusion equation (see Extended Data 1), we determined that each point within the striatum with a distance larger than 0.1 mm to its nearest neighboring arbor was defined as isolated. We therefore searched for spatially isolated areas where the innervation was sparse, since such areas would be more susceptible to impairments in DA signaling during denervation. Using Monte Carlo simulations, we approximated the distribution of smallest distances and used Voronoi tessellation to find the most isolated points (Fig. 1H; see Extended Data 1). This showed that no isolated areas existed in the fully innervated striatum of our model. This result is dependent on the size of the arbors as well as their numbers in the healthy striatum, and we thus tested different values of both and found the system to be quite robust.

These results demonstrate that the modelled dopaminergic arbors comprise a network that densely covers the striatum, where no isolated areas exist.

Different denervation patterns break down the dopaminergic network with distinct evolutions

In biology, structure often informs function. We therefore probed the spatial landscape of dopaminergic arbors in the denervating model striatum. The molecular pathways involved in the loss of dopaminergic neurons are beyond the scope of this study. Instead, we sought to characterize the organization of the remaining innervation arising from distinct models of progressive neuron loss.

To describe denervation, we assumed that all neurons have a rate of dying. All models were simulated using the Gillespie algorithm, and neurons were removed according to their rate. We first modelled the denervation to be independent of the spatial position of a neuron, meaning that all neurons had the same rate of dying. This we termed random denervation (RD; Fig. 2A). The second model, prion-like denervation (PLD; Fig. 2B), one neuron was initially infected. All infected neurons then had a rate to infect others and dying, while all uninfected neurons had zero rate of dying. Based on the rates, one infected neuron was chosen to infect two neighboring neurons before being removed from the network. This algorithm was based on the proposed mechanisms where protein aggregates spread between neurons and cause their degeneration (Prusiner, 2012; Chu and Kordower, 2015; Surmeier et al., 2017). In the third model, stress-induced denervation (SID), each neuron has a rate of dying calculated as a sigmoidal function of its number of neighbors. Thus, neurons with few overlapping arbors have a higher risk of dying compared with neurons with many overlapping arbors. This algorithm was based on the proposed mechanism where the remaining neurons may upregulate their DA synthesis and firing activity to maintain DA signaling. However, these neurons may already be close to their maximal metabolic capacity (Bolam and Pissadaki, 2012), and increased activity could trigger stress-induced degeneration (Jenner, 2003; Sulzer, 2007). We want to emphasize, that we do not consider these three models as mutually exclusive nor the “ground truth” mechanisms for the process of dopaminergic denervation. However, they each represent a simple algorithm for studying the self-organization of these complex phenomena and has the potential to give important insight in the different denervation structures.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Different denervation patterns break down the dopaminergic network with distinct evolutions. A–C, Diagrams of network mechanism for RD, PLD, and SID. The color of each dopaminergic neuron (circle) corresponds to probability of death. In C, dotted lines denote overlap of arbors. D, Visualization of the dopaminergic axonal arbor network following RD, PLD, and SID. Colors correspond to separate CACs. E, Distributions of the number of overlapping arbors for each individual arbor. F, Distributions of the number of arbors in each CAC. G, Distributions of the smallest distance to the nearest neighboring arbor center for each arbor. Dotted line denotes threshold for classifying isolated areas. In E–G, denervation is 80%. H, Fraction of remaining arbors as a function of time. I, Fraction of arbors belonging to the largest CAC as a function of denervation. J, Fraction of striatal space with smallest distance to nearest arbor larger than 0.1 mm (isolated area) as a function of denervation. In H–J, full line is mean, and shading is SD. RD, random denervation; PLD, prion-like denervation; SID, stress-induced denervation; CACs, contiguous arbor classes; AU, arbitrary unit. See also Extended Data Figure 2-1.

Using simulations, we observed how the three models resulted in distinct spatial landscapes, each characterized by a unique dopaminergic network breakdown (Fig. 2D–G). For RD, the remaining arbors covered the entire striatal space but no longer belonged to the same CAC. In contrast, for PLD, large fractions of the striatum were deprived of arbors and instead dominated by one or two subregions with seemingly normal innervation. For SID, arbors were concentrated in small, isolated subregions, each forming its own CAC. We quantified these observations by the distribution of the number of overlapping arbors (Fig. 2E). For PLD, a notable fraction of arbors had very low numbers of overlapping arbors, while a larger fraction had numbers like those in the healthy striatum. In SID, only arbors with many overlapping neighbors remained. Importantly, a commonality of all models was that the dopaminergic network broke down into multiple CACs, but in distinct patterns (Fig. 2F): RD had only small classes remaining, PLD contained many small but also one dominating class, whereas SID contained many classes containing 100 or more arbors. We also assessed the emergence of isolated areas (Fig. 2G). For RD, no isolated areas existed. In contrast, for both PLD and SID, the striatum contained numerous isolated areas, deprived of arbors.

Next, we followed spatial characteristics as a function of denervation. First, we determined the percentage of remaining arbors as a function of time (Fig. 2H). For RD, this followed an exponential decay with a relatively slow temporal progression. Interestingly, for PLD, the curve followed a convex function, suggesting that neuron loss accelerated with time, whereas the curve for SID followed a concave function, indicating that denervation in this scheme started fast, but then slowed with time. These results have predictive strength and can be mathematically described by stretched exponentials of the form: N(t)∝e−btc , with b being the decay rate and c = 1 for RD, c > 1 for PLD and c < 1 for SID. We next characterized the breakdown of the spatial network, by calculating the fraction of arbors in the largest CAC (Fig. 2I). PLD kept one dominating class until the final stage of denervation, while RD and SID were characterized by a tipping point, at which the network dramatically transitioned from fully coherent to segregated into multiple classes. Notably, this transition occurred around 75% denervation, which often correlates with the onset of symptoms in PD patients (Bernheimer et al., 1973; Fearnley and Lees, 1991; Ma et al., 1997; Lee et al., 2000). Finally, we probed the emergence of isolated areas, practically devoid of DA signaling. At 75% denervation, isolated areas comprised ∼50% and 20% of the striatum in the PLD and SID models, respectively (Fig. 2J). At the same denervation level, no isolated areas existed for RD, but these emerged at severe denervation.

Overall, we found notable spatial and temporal differences between distinct models of dopaminergic denervation. These differences between models were remarkably robust to changes in key parameters, that is, axonal arbor volume and numbers of dopaminergic neurons in the healthy state (Extended Data Fig. 2-1).

Dopaminergic denervation affects cAMP production and the activity of striatal SPNs

The excitability of SPNs is, along with several other factors, strongly regulated by DA (Surmeier et al., 2007; Kreitzer and Malenka, 2008; Kreitzer, 2009; Lahiri and Bevan, 2020). Thus, we next asked how dopaminergic denervation affects the activity of individual SPNs. Previous work has shown that D1 and D2 receptors have low and high DA affinity, respectively (Richfield et al., 1989; Fig. 3A). DA regulation of SPN excitability is mediated by the signaling molecule cAMP. D1 and D2 receptor activation increases and decreases the production of cAMP, respectively, and cAMP in turn regulates SPN ion channels (Surmeier et al., 2007; Kreitzer and Malenka, 2008; Kreitzer, 2009). Inspired by previous work (Dreyer et al., 2010), we described the intracellular cAMP concentration in D1- and D2-SPNs as: dcAMPD1dt=α + λ1DAhDAh+κ1h−δ1cAMPD1 dcAMPD2dt=α + λ2κ2hDAh+κ2h−δ2cAMPD2.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Dopaminergic denervation affects cAMP signaling and excitability of striatal SPNs. A, Diagram of how DA stimulates and inhibits the production of cAMP in D1- and D2-SPNs, respectively. B, Traces showing cAMP in D1- and D2-SPNs as a function of DA signaling. C, Maximal cAMP concentration in D1- and D2-SPNs during dopaminergic phasic firing and firing pauses, respectively, as a function of the number of dopaminergic terminals. D–K, Membrane potential dynamics of D1- and D2-SPNs in response to synaptic barrages (D, H), cAMP stimulation (E, I), synaptic barrages in combination with cAMP stimulation in the healthy state (F, J), or synaptic barrages in combination with cAMP stimulation in the 75% denervated state (G, K). Raster plots show the firing rate across time for 100 simulated neurons in each condition. DA, dopamine; D1, D1-class dopamine receptor; D2, D2-class dopamine receptor; SPN, spiny projection neuron; SNc, substantia nigra pars compacta. See also Extended Data Figure 3-1.

Here, α is the steady state production of cAMP, and δ is its spontaneous decay. In addition, receptor-dependent cAMP production was implemented: cAMP in D1- and D2-SPNs increased and decreased with DA stimulation, respectively. With increasing denervation, cAMP production during phasic firing became progressively lower in D1-SPNs, while in D2-SPNs it became progressively lower during firing pauses (Fig. 3B,C).

We next asked how these impairments in cAMP signaling may manifest in the activity of SPNs. For this, we used a previously published Hodgkin–Huxley-inspired model (Tatsuki et al., 2016; Rasmussen et al., 2017) to simulate the membrane potential (V_m) of D1- and D2-SPNs. This model contains extrinsic and intrinsic ion channel conductances. The extrinsic conductances are NMDA, AMPA, and GABA_A ion channels. The intrinsic conductances are voltage-gated and persistent Na⁺ channels, voltage-gated Ca²⁺ channels (Ca_V), and voltage-gated, leak, fast A-type, inwardly rectifying, slowly inactivating (K_SI), and Ca²⁺-dependent (K_Ca) K⁺ channels (Tatsuki et al., 2016; Rasmussen et al., 2017; see Extended Data 1). Using this model, we could closely mimic the V_m dynamics of SPNs. These neurons are characterized by their transitions between downstates and upstates (Wickens and Wilson, 1998; Sippy et al., 2015). In our model, on increased levels of synaptic barrages (implemented by increasing the stochastic noise of the V_m), both D1- and D2-SPNs transitioned into a brief upstate in which multiple action potentials were fired (Fig. 3D,H).

To implement the effect of DA, via its regulation of cAMP, on the V_m dynamics of D1- and D2-SPNs, we targeted the high-threshold Ca_V (N-type and P-type), K_SI, and K_Ca channels, which are negatively influenced by DA and cAMP signaling (Surmeier and Kitai, 1993; Nisenbaum et al., 1994; Surmeier et al., 1995, 2007; Shen et al., 2004; Lahiri and Bevan, 2020). Thus, for increasing cAMP levels, the conductance of these channels decreases and vice versa. We note that other channels may also be subject to dopaminergic modulation, such as persistent Na⁺ or NMDA channels, but to limit the parameter space, we here focus on the above mentioned Ca²⁺ and K⁺ channels. For stimulating D1- and D2-SPNs, we used the cAMP concentrations observed during dopaminergic phasic firing and firing pauses, respectively. This was motivated by the result that, in the healthy striatum, the maximal cAMP production in D1- and D2-SPNs was observed during these two phases respectively (Fig. 3B). In itself, cAMP stimulation was very rarely sufficient to evoke a transition from the downstate to an upstate in types of SPNs (Fig. 3E,I), supporting the notion of DA as a “modulator” rather than a “driver” (Surmeier et al., 2007; Kreitzer, 2009; Lahiri and Bevan, 2020). However, if DA and cAMP stimulation coincided with increased levels of synaptic barrages, this triggered a robust upstate that lasted longer and elicited more action potentials than with synaptic barrages alone (Fig. 3F,J); this demonstrates that cAMP powerfully regulates the excitability of D1- and D2-SPNs. As a result, in the denervated state, the activity of SPNs was notably affected (Fig. 3G,K): the duration of the upstate and the firing rate during the upstate were strongly diminished in both D1- and D2-SPNs. These findings were replicated using the simpler Izhikevich model (Extended Data Fig. 3-1), suggesting model-invariance.

Together, these data demonstrate that DA signaling, via its downstream effector cAMP, can regulate the firing activity of SPNs, and this regulation is impaired in the denervated state.

Distinct denervation patterns differentially affect global SPN firing activity

Next, we sought to investigate how different denervation patterns affect the activity of D1- and D2-SPNs across striatal space. For this, we spatially mapped the maximal firing activity for both types of neurons; D1-SPN during phasic firing and D2-SPN during firing pauses using the Hodgkin–Huxley-inspired model (Fig. 4A,B). In RD, although almost all the subregions had relatively low DA levels compared with the healthy striatum, this was still sufficient to evoke intermediate D1-SPN firing rates across the extent of the striatum. In contrast, in both PLD and SID, D1-SPN firing was high only in the subregions with preserved DA innervation. Noticeably, PLD transformed the striatum into a strongly polarized activity map, whereas SID caused local heterogeneity. For D2-SPNs, the emergence of isolated areas, resulted in a very different outcome. Since the maximal DA concentration in isolated areas is zero (except for small diffusive fluctuations), D2-SPN firing rates were here very high, most profoundly expressed for PLD and SID. We note here that, under physiological conditions, D2-SPNs residing in regions deprived of DA signaling might adapt by downregulating their firing rates to maintain homeostasis. In this scenario, the results would likely be like those for D1-SPNs.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Distinct denervation patterns differentially affect global striatal SPN firing activity. A, B, Maximal firing activity of D1- and D2-SPNs across space in the healthy and 75% denervated striatum for the three denervation patterns. C, D, Spatial mean and SD of maximal firing activity in D1- and D2-SPNs as a function of denervation. RD, random denervation; PLD, prion-like denervation; SID, stress-induced denervation; D1, D1-class dopamine receptor; D2, D2-class dopamine receptor; SPN, spiny projection neuron.

Finally, we characterized SPN activity in the three denervation models as a function of denervation. The mean D1-SPN firing rates decreased linearly as a function of denervation in all models (Fig. 4C). We also noted that the SD of D1-SPN firing was smaller in RD compared with both PLD and SID, indicating spatial homogeneity of firing levels (Fig. 4D). The effect on mean D2-SPN firing rates was different: firing increased notably for PLD and slightly for SID, as a function of denervation (Fig. 4C). In RD, the firing rates decreased until it reached a minimum around 80% denervation, whereafter it rapidly increased. This observation is explained by the occurrence of isolated areas, deprived of DA signaling, resulting in a dramatic increase in cAMP production in D2-SPNs (Fig. 3C), in turn resulting in a profound increase in excitability. Comparing the three models, the early progression of denervation (up to ∼60%) resulted in increased SD of SPN firing rates for all denervation patterns (Fig. 4D). This increase in activity variance across neurons may thus be a fingerprint of the denervating striatum. All the described results were fully replicated with the Izhikevich model (Extended Data Fig. 3-1).

Overall, these results show that the global firing activities of D1- and D2-SPNs are strongly affected by the specific spatial pattern of dopaminergic denervation.

A dual presynaptic compensation strategy preserves DA signaling in the denervated striatum

Given that dopaminergic neurons loss may trigger compensatory mechanisms in the remaining neurons in an attempt to maintain normal DA signaling (Zigmond, 1997; Brotchie and Fitzer-Attas, 2009), we sought to probe the potency of such mechanisms, to attempt predicting ideal therapeutic strategies. We included three presynaptic compensatory mechanisms as perturbed parameter values (V_M and Δ) in our original model for DA concentration dynamics. First, remaining dopaminergic terminals may increase their DA release capacity (Zigmond et al., 1990; Zigmond, 1997; Greenbaum et al., 2013). We refer to this as enhanced release compensation (ERC; Fig. 5A): Δ+↦Δ01−δ.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

A dual presynaptic compensation strategy preserves DA signaling in the denervated striatum. A, Diagrams of mechanisms of the ERC, DUC, and DEC models. B, Traces showing cAMP in D1- and D2-SPNs as a function of DA signaling at 80% denervation in the compensation models. C, D, Spatial mean and SD of maximal firing activity in D1- and D2-SPNs as a function of denervation pattern and compensation model. E, Spatial mean of maximal firing activity in D1- and D2-SPNs as a function of denervation and compensation model in the randomly denervated striatum. ERC, enhanced release compensation; DUC, decreased uptake compensation; DEC, dual enhanced compensation; NC, no compensation; DA, dopamine; D1, D1-class dopamine receptor; D2, D2-class dopamine receptor; SPN, spiny projection neuron; RD, random denervation; PLD, prion-like denervation; SID, stress-induced denervation. See also Extended Data Figures 5-1, 5-2.

Here, we introduced the compensation parameter δ, a sigmoidal function going from zero to one as a function of the number of dopaminergic arbors covering a small volume. The parameter Δ0 refers to the DA release in healthy subregions, whereas Δ+ is the compensated release value. Second, DA transporters, expressed on terminals, may reduce their uptake capacity (Zigmond et al., 1990; Zigmond, 1997; Lee et al., 2000; Greenbaum et al., 2013). We refer to this as decreased uptake compensation (DUC; Fig. 5A): V−↦V0(1−δ).

As above, the parameter V0 refers to the uptake value in healthy subregions, whereas V− is the compensated uptake strength. Finally, we suggest a mechanism where neurons compensate by enhancing both DA release and uptake capacity in the terminals. This idea stems from the equations for DA concentration dynamics, from where we can show mathematically that this combination recovers the original equations. Such a compensatory mechanism has not previously been suggested, and we refer to this as dual enhanced compensation (DEC; Fig. 5A); this is included in the model through changes in both the uptake and release parameters: Δ+↦Δ01−δandV+=V01−δ.

Here, all parameters are defined as above.

With these implementations, we simulated DA signaling and the corresponding cAMP production in D1- and D2-SPNs with 80% denervation. As shown above, in the absence of compensation, DA release during tonic firing is unaffected, but notably affected during phasic firing and firing pauses (Fig. 3B). Here, during tonic firing, the DA concentration notably increased for ERC and DUC models; during phasic firing, DA was increased for ERC, and during firing pauses, DA removal was incomplete for DUC (Fig. 5B). Importantly, the DEC model preserved DA levels during both tonic and phasic firing at comparable levels to in the healthy state, while still allowing complete DA removal during firing pauses (Fig. 5B). We also tested several postsynaptic compensatory mechanisms previously proposed in the literature, including increased D2 receptor expression (Guttman and Seeman, 1985), enhanced D1 and D2 receptor sensitivity (Lee et al., 1978), and suppressed cAMP degradation in D1- and D2-SPNs (Niccolini et al., 2015). These mechanisms were inadequate to restore DA and cAMP signals in SPNs (Extended Data Fig. 5-1), and we therefore did not explore these further.

Next, we asked whether any of the compensation mechanisms were able to counteract the impairments in SPN firing activity in the denervated state (Figs. 3, 4). For this we calculated the spatial mean and SD of the maximal D1- and D2-SPN firing rates in the three denervation models and combined these with the presynaptic compensation mechanisms (Fig. 5C). Interestingly, in the scenario of RD, only the DEC model preserved the mean level of D1- and D2-SPN firing rates. In contrast, for PLD and SID, none of the three compensation models were able to counteract the decrease in D1-SPN firing with denervation, while all models performed relatively well for D2-SPN firing activity. For the SD of the D1-SPN firing rates, we note that in the RD scenario, all compensation models, as well as the noncompensation state, maintained this measure near the healthy level (Fig. 5D). In contrast, for PLD and SID, the SD was notably increased for all compensation models, and curiously, the noncompensated state was most similar to the healthy state. For the SD of the D2-SPN firing activity, none of the compensation models truly maintained this measure close to the healthy level, regardless of the denervation pattern (Fig. 5D). It is here worth noting that the DEC model, across all denervation patterns, maintained the SD of D2-SPN firing activity at a very low level. This is because, in regions with low dopaminergic coverage, DA signaling from remaining terminals in the DEC model can compensate, restoring coherent neuronal activity. The low SD in firing activity across neurons means that all striatal subregions are capable of generating a very similar firing response on dopaminergic stimulation. Overall, we conclude that the DEC model, in combination with the RD pattern, best preserved the SPN firing activity. In the final set of simulations, we thus explored this for all levels of denervation (Fig. 5E). For the DEC model, the firing rates of both D1- and D2-SPNs remained remarkably close to the healthy state, despite reaching severe denervation. In contrast, for the ERC and DUC models, even at intermediate denervation, SPN firing differed from the healthy state. Therefore, the ERC and DUC mechanisms do not seem ideal as therapeutic strategies. These findings were replicated with the Izhikevich model (Extended Data Fig. 5-2).

Taken together, these results show that an ideal strategy to maintain normal SPN activity is to locally introduce a dual compensation mechanism, increasing both DA release and uptake capacity, and to globally minimize the dopaminergic arbor density differences, or at least avoid the emergence of isolated areas.