Predicting complex spikes in striatal projection neurons of the direct pathway following neuromodulation by acetylcholine and dopamine

The input structure of the basal ganglia, striatum, receives dense neuromodulatory input in the form of dopamine and acetylcholine. The two systems are tightly connected, for example, synchronized activity of cholinergic interneurons, leading to increased acetylcholine release, has been shown to directly trigger dopamine release from dopaminergic terminals in striatum. Both signals are further needed for induction of locomotion. High dopamine concentration leads to increased excitability of the direct pathway striatal projection neurons. High cholinergic tone inhibits various potassium channels further increasing the excitability of striatal projection neurons. Here, we investigate the combined effect of concurrent high acetylcholine and dopamine using biophysically detailed models based on rodent data. The aim of the study is to investigate how neuromodulation affects dendritic integration. The result shows that neuromodulation paired with synaptic activation of dendrites can give rise to complex spiking patterns, resembling spike shapes seen in the hippocampus. In the hippocampus, these complex spikes are associated with behavioral time scale plasticity and place cell tuning. We further investigate the mechanisms behind the complex spikes and find that there are two components, one axo‐somatic and one dendritic in origin.


| INTRODUCTION
The striatum is the main input structure of the basal ganglia, and the striatal projection neurons (SPN) integrate glutamatergic input from a large part of cortex and thalamus. They also receive strong neuromodulatory input of which the two most well studied are the dopaminergic and the cholinergic. The dopaminergic input originates outside of striatum while the cholinergic primarily is produced by the intrastriatal cholinergic interneuron (ChIN).
Dopamine (DA) can act both as a learning signal and motivational cue while ChINs has been proposed as the switch between the two modes (Berke, 2018). At the same time, both signals are required for induction of locomotion (Howe et al., 2019). Transiently decreased activity of ChINs also gives decreased excitability of SPN in vivo (Zucca, Zucca, Nakano, Aoki, & Wickens, 2018).
Both dopamine and acetylcholine act through metabotropic receptors in the cell membrane, triggering intracellular signaling cascades, for example, leading to changed phosphorylation levels of ion channels. The SPNs can be further divided into subtypes based on expression of dopamine receptors. SPNs of the direct pathway (dSPN) express the dopamine receptor type 1 and project directly to the basal ganglia output targets. The second subtype projects indirectly to the output targets (iSPN) and expresses the type 2 dopamine receptor.
The main dopaminergic cascade in dSPN starts with receptor activation and leads to downstream activation of protein kinase A, etc. This pathway is inhibited by the activation of muscarinic receptor type 4 (M4R), activated by acetylcholine (Bruce et al., 2019). Another muscarinic pathway is, however, most commonly reported during cholinergic modulation. This pathway starts with the type 1 muscarinic receptor (M1R) and leads to the activation of protein kinase C (PKC). For a recent review of the cholinergic system in striatum, see Abudukeyoumu, Hernandez-Flores, Garcia-Munoz, and Arbuthnott (2019) In this study, we test the combined effect of the two neuromodulators acetylcholine and dopamine on dSPN under synaptic bombardment and obtain complex spike activity similar to what has been reported in hippocampus (Bittner et al., 2015). In the hippocampus, these complex spikes are observed together with spontaneous plateau potentials leading to behavioral time scale plasticity and place cell induction (Bittner et al., 2015;Bittner, Milstein, Grienberger, Romani, & Magee, 2017). These plateau potentials are N-methyl-D-aspartate (NMDA)-dependent events triggered in the dendrites by spatio-temporal clustered glutamatergic input. Similar dendritic non-linearities have been observed in SPNs (Du et al., 2017;Plotkin, Day, & Surmeier, 2011).

| Model
The set of models used here are based on a model previously used (Lindroos et al., 2018). Here, however, the distribution and maximal conductance of ion channels were instead randomly set based on computational techniques (Prinz, Bucher, & Marder, 2004;Taylor, Goaillard, & Marder, 2009). Additionally, the muscarinic K + current (I m ) was added to the axon initial segment (AIS; Petersen, Cotel, & Perrier, 2017) since it is expressed in SPNs, and is a major target for cholinergic neuromodulation (Shen, Hamilton, Nathanson, & Surmeier, 2005;Shen et al., 2007). The passive parameters and the morphology of the neuron model were not updated.

| Morphologies
The morphology used in the iSPNs was downloaded from neuromorpho.org: http://neuro morpho.org/neuron_info.jsp? neuron_name=WT-P270-09 (Martone et al., 2003). The following corrections were also made: The dendrites were corrected for shrinkage, the diameter of the dendrites was rescaled (following Lindroos et al., 2018), the multipoint soma was exchanged with a single point (with diameter 12.1 µm; Reinius et al., 2015), and an AIS was added (length 60 µm, diameter 1 µm). The morphology of the dSPNs was the same as was used in Lindroos et al., (2018).

| Optimization
During model optimization, multiple model candidates were created by drawing parameters randomly from predefined ranges of different conductance distributions, see Table 1. Each candidate was validated against experimental data regarding dendritic calcium (Ca) response following a backpropagating action potential (dCa-bAP, Figure 1a4) and somatic current-frequency response (FI, Figure 1a3). For the dCa-bAP protocol, one action potential was triggered in the soma using a short-duration, high-amplitude current injection (2 ms, 2 nA), while monitoring the Ca concentration in the | 2119 LINDROOS aND HELLGREN KOTaLESKI dendrites. The peaks of the following transients were extracted and normalized to the response in the proximal dendrites (at about 40 µm distance). Finally, an exponential function was fitted to the measurements and compared to the experimentally measured counterpart (Day, Wokosin, Plotkin, Tian, & Surmeier, 2008). The comparison was done by calculating the difference ("error") between the two curves at every 10 µm from 40 to 200 µm somatic distance. If the sum of the squared error was large (>0.2), the candidate was rejected, if not, the FI protocol was run. The FI protocol consisted of a longer duration current (900 ms) with amplitude of typical rheobase currents in SPN (Planert, Berger, & Silberberg, 2013). Here, the lower and upper boundaries for the rheobase were set to 240 and 480 pA, respectively, for dSPN and 140 and 380 for iSPN. If the model candidate spiked for the low bound injection, or did not spike using the high bound injection, it was rejected. A candidate was also rejected if it spiked more than 70 spikes for a single injection, had an afterhyperpolarization crossing −65 mV, had oscillations in the membrane potential leading up to spikes, or if the slope of the FI curve was not within 3 standard deviations from the experimentally measured mean. A model candidate had to pass all tests to be saved to the model library and used during simulations. The response of an example model to suprathreshold current injections of increasing magnitude is shown in Figure 1a2.
The random parameter set and the resulting ion channel distributions of the optimized models are shown in Figure 1a5, top and bottom, respectively (for dSPNs). Equations and parameter ranges of the dendritic ion channel distributions are given in Table 1.

| Ion channels
Some of the ion channels from the parent model included "pointers" (Lindroos et al., 2018), linking them to an intracel-lular cascade. These were all removed, and direct modulation was instead included (see Neuromodulation). Additionally, the I m channel was added. Since no striatal I m channel model existed, the choice fell on a model used in an influential pyramidal neuron model (Hay, Hill, Schürmann, Markram, & Segev, 2011), based on bullfrog data (Adams, Brown, & Constanti, 1982).
The distribution of ion channels in the dendrites followed the parent model (Lindroos et al., 2018), with low expression of the sodium channel (Naf) beyond the proximal dendrites and transient potassium (Kaf, Kv4.2) increasing with somatic distance. Together, these two distributions generated a valid dCa-bAP-response. The slow A-type potassium channel (Kas, Kv1.2) and the N-type Ca channel (CaN) also decreased with somatic distance while the T-type Ca channels increased (see Figure 1c, bottom). These non-uniform distributions were modeled using sigmoid functions: distribution = p 1a + p 1b / (1 + exp{(d-p 2 )/p 3 }), where d is the somatic distance along the dendrite and p 1-3 are randomly drawn parameters. The parameters p 1a and p 1b set the magnitude of the steady and transient part of the curve, respectively. The parameter p 2 sets the half transition distance, and p 3 sets the steepness of the transition. Additionally, the whole curve, as well as the uniformly distributed inward rectifying potassium channel (Kir) and the small Ca activated potassium channel (Sk), were scaled by another random factor, p 0 . The ranges for p 0-3 are given in Table 1, together with the equations for each channel and the resulting distribution curves are plotted in Figure 1c.

| Calcium pools
As in the parent model, two separate Ca pools were used. One connected to the N-P-Q-R-type Ca channels and the other to synaptic channels (α-amino-3-hydroxy-5-methyl-4-isoxazoleprop- Not<del author="Jeanette Hellgren Kotaleski" command="Delete" timestamp="1594793505392" title="Deleted by Jeanette Hellgren Kotaleski on 7/15/2020, 8:11:45 AM" class="reU3">e</del>: gbar is the maximal conductance from the parent model (Lindroos et al., 2018), and d is the somatic distance. See Figure 1c for resulting distributions and illustrations on how the ranges are spanned. The channels are the fast sodium channel (Naf), the fast A-type potassium channel (Kaf), the slow A-type potassium channel (Kas), the inward rectifying potassium channel (Kir), the small conductance potassium channel (calcium dependent, Sk), the N-type calcium channel (CaN), and the two T-type calcium channels (subtypes 2 and 3, CaT3.2 and CaT3.3, respectively).
ionic acid, AMPA and NMDA) as well as L-type and T-type Ca channels. The minimum Ca concentration was set to 10 nM for both pools.

| Synaptic input
Synaptic inputs came in two forms. One type of input was spread out synaptic bombardment, where each section of the morphology was given one randomly activated synapse of each type (glutamate and gamma-aminobutyric acid, GABA). The primary purpose of this input was to depolarize the cell, mimicking ongoing activity seen under in vivo like conditions, but was also used to drive the cell to spike by gradually increasing the frequency of the input (ramping). The second type of input was a specific pattern of activation, where 16 synaptic activations were given sequentially with an interspike interval (ISI) of 1 ms. These 16 activations were given in 1 to 16 discrete locations, and the purpose of this input was to trigger non-linear processes in the dendrites (Du et al., 2017;Plotkin et al., 2011). Gradual increasing spatial clustering was implemented by starting with 16 randomly F I G U R E 1 Validation of models and neuromodulation. In a 1, the morphology of the models is shown. In a 2 , example traces from one model version are shown following a set of suprathreshold current injections (rheobase to rheobase + 90 pA, in steps of 15 pA). The same protocol is used in a 3 where the FI curve of all models (N = 54) is compared to experimental data (Planert et al., 2013). In a 4 , dendritic excitability is validated as the relative calcium change following a backpropagating action potential and plotted as a function of somatic distance. The protocol used here is a short duration, high amplitude current injection (2 ms, 2 nA), exemplified in the inlay. Gray dots are individual calcium peak values at a dendritic location (normalized to the proximal values); black lines are the best fit to these measurements for a single model version, and the red line is experimental data (Day et al., 2008). In a 5 , the distribution of random factors over their respective intervals (normalized from 0 to 1) is shown, together with the resulting ion channel distribution over the dendrites (top and bottom rows, respectively). In b 1 and b 2, the dopaminergic modulation is illustrated and validated, respectively (Planert et al., 2013). In b 3, the different modulation protocols are investigated. The number of trials, N, in each protocol is 10,800 (54 models with 200 trials each). Uniform modulation over cellular compartments with low (± 5%) modulation of the fast potassium channel (Kaf) makes the cell less excitable (0.1% of the trials are more or equal excitable). Increasing the modulation of the Kaf channel increases the percentage of the trials that are more or equal excitable to 28.6% (b 3 , left panel). If additionally the axon initial segment is not modulated, most trials are more or equal excitable (91.1%, b 3 , right panel). Non-uniform modulation without modulation of the kaf channel results in 1.55% trials with increased or equal excitability. Color and size of the dots are given by the size of the modulation factors for the inward rectifying potassium channel (Kir) and slow A-type potassium channel (Kas), respectively. In c, the cholinergic modulation is validated against experimental data. c 1 illustrates the setup. c 2-4 shows the resulting effect on excitability (Pineda et al., 1995), spike induced afterhyperpolarization (Pineda et al., 1995), and inward rectification (Figueroa et al., 2002), respectively. The protocol used in c 2-3 is a constant current injection of rheobase + 15 pA amplitude. In c 4, subthreshold injections of −600 to + 300 pA amplitudes are used. The red arrows indicate the direction described in the literature (see Table 2)

| 2121
LINDROOS aND HELLGREN KOTaLESKI distributed synapses and then on each iteration removing half of them. One of the remaining synapses inherited spike times from one removed synapse (with the synapse model used here, this is equivalent to moving the synapse to the same isopotential dendritic subcompartment). The procedure was repeated until only one synaptic site remained (clustered input) with all the 16 incoming spikes. Simulations were repeated on each iteration. Technically, the synaptic input was given based on the Exp2Syn object in Neuron. For the glutamatergic input, a new mechanism was implemented, including both the AMPA current and the magnesium-blocked NMDA current. For this new mechanism, time constants of AMPA and NMDA followed (Hjorth et al., 2020). The maximal conductance used in the specific pattern protocol was 1,500 pS and 500 pS for NMDA and AMPA, respectively, while both NMDA and AMPA were set to 300 pS in the synaptic bombardment protocol. The maximal conductance of GABA was set to 900 pS in the synaptic bombardment and 1,500 pS when specifically investigating inhibition (see Section 3). The included synaptic model did not include short-term dynamics and was therefore simpler than the one used in (Hjorth et al., 2020).

| Neuromodulation
Neuromodulation was implemented as a change of ion channel conductances. Channels were either updated at the beginning of the simulation (static), mimicking bath application of neuromodulators, or transiently during the simulation. The static protocol was used to assess the combined effect of neuromodulation on cellular behavior. Here, sets that modulated the behavior of the cell according to the literature were selected T A B L E 2 Cholinergic modulation of striatal projection neuron (SPN), summary of the literature review  Pineda et al. (1995) and saved for the transient protocol. In the transient protocol, the conductance then varied between the non-modulated control and the maximal value obtained in the static protocol. The transient protocol also included dopaminergic modulation of synaptic channels. Simultaneous dopaminergic and cholinergic modulation was set in an additive fashion, using the formula: modulation = 1 + (T DA -1)*l DA + (T ACh -1)*l ACh , where T x was the target value for the respective neuromodulator and l x gave the instant value (by scaling it from 0 to 1). The combined modulation calculated this way, was then multiplied onto the conductance of the channel. A modulation of 1 represents no modulation in this setup, while a smaller or larger value represents a down-or upregulation, respectively. For example, if dopamine reduces a channel by 20% and acetylcholine also reduces the same channel by 30% the individual target values would be 0.8 and 0.7, giving a modulation factor of 0.5 when fully activated (I DA = I ACh = 1). In order to avoid negative conductance, the minimal value of the combined modulation was set to 0. For cholinergic modulation of single channels, see Table 2, and for dopaminergic modulation, see older studies (Hjorth et al., 2020;Lindroos et al., 2018). Two sets of each DA and ACh were used during simulation while synaptic channels were modulated with single representative values in simulations investigating mechanisms of complex spikes (reducing run time and hard drive storage). Simulations on the direct involvement of specific factors used ranges.
Each set was randomly drawn from ranges (see Table 3), and the resulting model behavior was validated against experimental data using the FI protocol with current amplitudes of rheobase + 0-12 pA (Figure 1b2, Planert et al., 2013). To increase the excitability of the models, the Kaf channel was reduced (15%-25%) and the AIS was not modulated (Lindroos et al., 2018). This protocol was validated against protocols with uniform modulation (including modulation of the AIS) and without modulation of the Kaf channel (allowing a ± 5% change), and was found to be the only protocol to reliably increase the excitability of the cells (Figure 1b3, Lindroos et al., 2018;Planert et al., 2013;Surmeier, Carrillo-Reid, & Bargas, 2011).

| Cholinergic modulation
Cholinergic modulation (acetylcholine, ACh) was implemented using the same technique as for the DAergic modulation, except for the Kaf channel where the activation and inactivation curves were shifted into more negative potentials (rather than scaling the conductance, Akins, Surmeier, & Kitai, 1990). Random modulation factors were drawn from ranges reported in the literature (see Tables 2 and 3), applied to the models in the library, and the resulting behavior was validated against experimentally reported changes. These changes were increased discharge (Figure 1c2, Pineda, Bargas, Flores-Hernández, & Galarraga, 1995), decreased spike afterhyperpolarization (Figure 1c3, Pineda et al., 1995) and decreased inward rectification (Figure 1c4, Figeroa, Galarraga, & Bargas, 2002). Sets that pushed the models in the right direction were collected in a library and used in simulations.

| Spike sorting
Action potentials were extracted from the voltage trace using the same function as in the parent model (Lindroos et al., 2018) with a threshold of 0 mV. Complex spikes were detected in a low pass filtered version of the membrane potential, by checking if a sliding average crossed a threshold. A second order low pass butterworth filter was used, implemented using the scipy.signal libraries butter and filtfilt and the threshold was set to −37 mV. F I G U R E 2 Neuromodulation induces complex spikes in dSPN. The models are activated using synaptic patterns of different spatial clustering. An example of these patterns is shown in a. The arrows in the middle panel illustrate how activated synapses are gradually moved from being distributed over many sites (16 × 1, left panel) to all being activated in a single site (1 × 16, right panel). The color (red to black) also indicates low to high degree of clustering. If activated with no other stimuli, the most spatially clustered pattern gives the highest amplitude of the depolarization (b, left panel). The amplitude of the clustered input is further dependent on the somatic distance of the location (right panel), and different model versions have peak values at different somatic distances (small yellow dots). In c, the models are bombarded with additional synaptic input (glutamate and GABA), depolarizing the somatic membrane to about −70 ± 10 mV. Under these conditions, the most clustered input is still most efficient in driving the cell (right panel). In d, additional neuromodulation is added, which increases the spike probability of the cell. Dopamine alone (green trace) or combined with acetylcholine (purple trace) is increasing the spike probability the most. In e, the ramping activation protocol is illustrated (ramping synaptic bombardment, without additional input). In f, the results from this protocol are shown. The right column shows the spike probability and time to first spike under different combinations of neuromodulators. The left column shows example traces including an example of the complex spike pattern resembling the complex spikes seen in the hippocampus (Bittner et al., 2015) 2.11 | Experimental data Experimental data were extracted from published articles using an online extraction tool, the WebPlotDigitizer https:// autom eris.io/WebPl otDig itize r/.

| Synaptically driven dSPNs can trigger complex spikes during neuromodulation
In an earlier study, we showed that spatio-temporally activated glutamatergic input (from here on referred to as clustered input) can trigger dendritic NMDA spikes in the distal dendrites, but not in the proximal (Du et al., 2017).
Here, we further developed that framework by investigating the response to clustered input over the whole morphology of the models. As control, we used the same number of synapses, but randomly distributed, and with the same activation strengths. For a description of the activation patterns used, see Section 2, and for an illustration, see Figure 2a.
We first activated patterns of 16 glutamatergic synapses alone and quantified the amplitude and duration of the resulting depolarization. The duration was quantified as the full width of the depolarization at half amplitude. The results showed that clustered input gave a slightly larger amplitude than distributed (around 9 mV on average for both), but with a large difference in duration (about 20 ms for distributed and 35 ms for clustered), see Figure 2b, left row. Sorting the clustered response on somatic distance revealed that the duration increased with distance while the amplitude had a maximum at around 70-80 µm distance and then decreased with distance ( Figure 2b, right row). The decreased amplitude was not necessarily a monotonous function of distance from the soma, and not all models had peak values for the same distance, see yellow dots. This shows that both the morphology and ion channel distribution along the dendrites dictate the response to a given signal.
Similarly, clustered input also gave a higher spike probability if the patterns were activated from more depolarized potentials. The additional depolarization was achieved by adding synaptic background noise, GABAergic and glutamatergic synaptic input, spread over the morphology (synaptic bombardment, see the Section 2). Only a few models spiked using this stimulation protocol but the general trend with a peak around 70-80 µm from the soma still persisted (Figure 2c).
If additionally neuromodulation was added to the protocol, spike probability increased, but both relationships still held. Clustered input was more efficient, specifically if triggered at 70-80 µm somatic distance. The neuromodulation was here activated using a step like sigmoidal function activated simultaneously as the first synapse in the pattern (during the synaptic activation period, SAP, see Figure 1a). The neuromodulator DA (green traces alone, or purple traces combined with ACh) increased the response most compared to control (black traces) or ACh alone (orange traces, Figure 1d).
If we instead gradually increased the frequency of the synaptic bombardment (illustrated in Figure 1e), DA + ACh still gave the highest spike probability and the fastest time to first spike ( Figure 1f).
Furthermore, during neuromodulation the spike shape sometimes changed into a shape resembling the complex spikes seen in pyramidal neurons of CA1 (left panel of Figure 1f, purple trace, Bittner et al., 2015). These complex spikes were at rare occasions also seen with the activation of the synaptic pattern, illustrated in later figures.
To summarize this section, both morphology and ion channel distribution along the dendrites dictates the somatic response of dendritic excitation. If the excitation is large, it can trigger complex spikes in the soma.

| Dendritic NMDA spikes are triggered at the onset of complex spikes
Complex spikes (CS) are associated with increased NMDA current in the dendrites (Figure 3). Specifically, NMDA spikes are triggered in a few dendritic branches (light gray traces of Figure 3a1-2 and c1-2), leading to increased Ca influx (Figure 3a3 and c3). The increased Ca concentration is much larger than the influx caused by the burst of action potentials leading up to the CS (AP, marked with arrows in Figure 3a3 and c3).
Further, the CS initiation zone (the minimum voltage of the last afterhyperpolarization before the CS, Figure 3a1) is close to the NMDA-spike threshold (the peak acceleration point of the voltage-dependent magnesium block, d 2 Mg/dt 2 , Figure 3a1 and b1). The average dendritic NMDA current is also highly correlated with the low pass filtered membrane potential of the soma, both during and before the ramp (Figure 3b2). The AMPA current is instead decreased during CS and not correlated with the somatic membrane potential (Figure 3a2, b2, c2).
The duration of the CS is generally longer following ramping input compared to the ones triggered using the activation pattern (Figure 3b3). The reason for the longer duration is sequentially triggered NMDA spikes in the dendrites using the ramping protocol, exemplified in the left panel of Figure 3b3.
To summarize this section, NMDA spikes are triggered during complex spikes, leading to large calcium influx in the dendrites.

| Persistent activation of the sodium window current
The duration of the complex spike is not only set by dendritic processes, there is also an axo-somatic component. During complex spikes, the local membrane in the axon is depolarized beyond the maximal value of the sodium window current (see Figure 4a). This causes a situation where a hyperpolarizing drive will cause an increase in the persistent sodium, which stabilizes and prolongs the depolarized state. This is exemplified in Figure 4b. See specifically the p1 line in the right panel, showing a zoomed in version of the areas marked with dotted boxes in the left panel of Figure 4a,b. As the drive to hyperpolarize increases, the membrane gradually crosses the peak and the current decreases (see p2-p4 lines).

F I G U R E 3
Complex spikes are associated with multiple dendritic NMDA spikes. In a, the ramping protocol (see Section 2 and Figure 2e) is used to trigger a complex spike (CS). a 1 shows the somatic and dendritic voltage during the ramping phase together with an illustration of the setup. The last afterhyperpolarization before the CS is also marked with an arrow (cs initiation zone). The right panel shows the voltage dependence of the magnesium block (Mg-block, black) together with its first and second derivative (blue and gray, respectively). The maximum of the second derivative (−42 mV) is also shown as a gray line through the voltage plot. In a 2, the glutamatergic current driving the cell is shown. The NMDA current (gray) is increasing throughout the ramp, but spikes during the complex spike. The AMPA current (green), instead decreases during the complex spike (due to reduced driving force). In A3, the resulting intracellular calcium concentration is shown for the two pools in the models (see Section 2). The timing of action potentials (AP) is also marked with arrows. In b 1 , the voltage of the CS initiation zone is shown for all complex spikes together with the second derivative of the Mg-block. All values are close to the peak, where the Mg-block changes the fastest. In b 2, the Pearson correlation coefficient is shown for current versus membrane potential for the full trace (−500 to 500 ms), the flat part before the ramp (−500 to 0 ms) and the ramping (0-500 ms). The NMDA current is highly correlated with both phases (flat and ramping). In b 3, the duration of the CS (full width from the CS initiation voltage, illustrated in the inlay) is plotted for the two stimulation protocols (ramping and clustered pattern). The duration is generally longer for the ramping protocol due to sequentially triggered NMDA spikes, left panel. In c 1-3, the same figures as in a 1-3 are shown for CS triggered using a spatio-temporal clustered pattern (see Section 2 and Figure 2a). The dendrites where the synapses are activated are highlighted using arrows 2126 |

LINDROOS aND HELLGREN KOTaLESKI
The hyperpolarizing drive originates in the dendrites and/ or soma since the difference in membrane potential between the soma and the axon increases following the peak of the complex spike (see the bottom right panel of Figure 4a, p1-p4 lines).
To summarize this section, persistent sodium current in the AIS prolongs the duration of the complex spike by counteracting the hyperpolarizing drive, originating in the dendrites, with increased current.

| Neuromodulatory induced changes in spike machinery leads to complex spikes
The probability of getting complex spikes in our data is generally low. Using the ramping protocol, around 3% (36/1080) of the traces were found to contain complex spikes under combined neuromodulation (ACh + DA, Figure 5a). The probability of regular spikes was on the contrary high, at 97% (1046/1080). We also found complex spikes with a low probability under only dopaminergic modulation (DA, 4/540). No complex spikes were obtained under control conditions (ctrl) or acetylcholine (ACh). The regular spike probability was also lower under these three conditions (DA 456/540, ACh 71/540, and control 19/270, Figure 5a). Under the clustered activation protocol, both types of spike probabilities were lower (n = 32,400; regular, ACh + DA 63202/4n, DA 24535/2n, ACh 625/2n and control 159/n; complex, ACh + DA 475/4n, DA 35/2n and zero for ACh and control).
Further, only a few of the models in the library triggered complex spikes under combined neuromodulation (ACh + DA) under the standard dopaminergic modulation protocol (with no dopaminergic modulation of the axon initial segment, Figure 5c1-2, top row). If we instead applied the uniform dopaminergic protocol, the regular spike probability went down, but many more models triggered complex spikes (49 versus 8 with the standard dopaminergic protocol), and with a higher percentage (Figure 5c1-2, bottom  row). Since the AIS contains a high expression of sodium channels, we performed simulations where the modulation factor for the sodium channel was allowed to vary (between 0.7 and 1.0), while the other factors stayed at their respective values. This showed that a low value of the factor increased the proportion of complex spikes under uniform modulation (Naf, Figure 5c3, bottom panel). Under the standard protocol, there was no clear such a relationship (Figure 5c3, top panel).

F I G U R E 4
Complex spikes are associated with persistent activation of the sodium window current. In a 1, the somatic (black) and axonal initial segment (AIS, purple and red) voltage traces are shown for an example complex spike (CS) together with the sodium channel window current (right panel). The maximum of the window current (−25 mV) is shown as a gray line through the voltage plot. In a 2, a magnification of the area in the box of a 1 is shown together with the difference in voltage between the three compartments (distal axon minus proximal axon, and proximal axon minus soma; bottom panel). The difference increases during the second half of the CS (from p1), indicating that the dendrosomatic region is hyperpolarizing first. In b 1, the sodium current in the three compartments plotted in A are shown. In b 2, the area in the box of b 1 is magnified (bottom) together with the area in a 1 (top). As the membrane is gradually hyperpolarized, the sodium current is pushed over the peak of the window current and reduces (p1-4). When the distal part of the AIS passes the peak, the CS rapidly ends (p3-p4) Repeating the same type of simulation for the fast potassium channel (Kaf) also gave an increased proportion of complex spikes when the modulation factor was low (for both modulation protocols, Figure 5c).
Finally, we investigated if complex spikes could be triggered using current injection under concurrent dopaminergic and cholinergic modulation (ACh + DA). Under the standard protocol, there were no complex spikes following current injection, but if strongly stimulated, the cells instead went into depolarization block. Average current needed to push a model into depolarization block ranged from about 40-400 pA above the unmodulated rheobase current (Figure 5d, top row). However, under the uniform dopaminergic protocol complex spikes were triggered in about ⅔ of the models (37/54). Specifically models with low peak amplitude of action potentials, indicating a low sodium current, were prone to trigger complex spikes (Figure 5d, bottom row).
Together, the results above indicate that increased excitability alone does not lead to increased complex spike probability. On the contrary, decreased sodium current seems to be directly involved in complex spike formation.

| Complex spikes in iSPN are also dependent on reduced sodium current
In contrast to the effect on dSPN, the excitability of the iSPN is reduced by dopamine (Hernandez-Lopez et al., 2000;Planert et al., 2013). Modulation of Na and NMDA are both inconsistently reported in the literature for this SPN subtype (see supplementary table 8 in Hjorth et al., 2020). Here, we allowed a small modulation in either direction (see Table 3), simulated in models with iSPN characteristics obtained using the same technique as for the dSPN (n models = 34, see the Section 2 section and Figure 6a). Under these conditions, F I G U R E 5 Decreasing the sodium conductance increases the probability of complex spikes. In a, right panel, the spike probability is shown for the ramping protocol. Complex spikes (CS) are triggered in 3% of the trials under concurrent neuromodulation (acetylcholine, ACh and dopamine, DA). The left panel illustrates the setup, and the middle panel shows two spiking example traces (regular, gray and complex, red). In b, the same is shown for the spatio-temporally clustered pattern (without example traces). Here, the CS probability is only 0.4%, and regular spike probability is also lower compared to the ramping protocol (gray bars). In c, the spike probability is shown over individual models under concurrent ACh + DA modulation. The top row shows the standard modulation protocol, without DAergic modulation of the axon, and the bottom row uniform DAergic modulation. The regular spike probability is lower with uniform modulation but the CS probability is higher. In C3, the factors for kaf and naf are redrawn from the interval of 0.7 to 1, and CS probability is plotted over factors. Using the standard protocol, there is a small increase as the fast A-type potassium channel is reduced (Kaf), while for the sodium channel (Naf), the trend is not clear (top). Under uniform modulation, however, the CS proportion increases as the Naf channel is decreased. To a lesser extent, this is also true for the Kaf channel. In D, current injection is used to trigger CS under uniform modulation. Also, here the models with low sodium current are more likely to trigger CS as indicated by the low mean peak voltage of the traces with CS than without (bottom-right). Under the standard modulation protocol, CS cannot be triggered in the models, but if the amplitude of the injection is large, it can push the cell into spike depolarization block. The current needed to do this varies with model and modulation set, but is typically around 200 pA above the unmodulated rheobase (40-400 pA, upper-right) DA decreased while ACh increased the spike probability following ramping input (Figure 6b1-2), but none of these simulations triggered complex spikes (n trials >6,000 for each condition). However, if we reversed the cholinergic modulation of the sodium channel, either alone or together with increased NMDA we could trigger complex spikes (Figure 6c1-2). Plotting the modulation factors of Na and NMDA in these traces with complex spikes showed that complex spikes were most common if Na was maximally reduced (30%) and NMDA was maximally increased (60%, Figure 6c3). In summary, with the current models of the two types of SPNs, we thus predict that dSPNs are more likely to show complex spikes under neuromodulation than iSPNs.

| Inhibition with slow dynamics reduces the complex spike probability the most
In a previous study, we showed that on-site inhibition was the most efficient inhibition to prevent or terminate dendritic NMDA spikes (Du et al., 2017). Here, we extend that investigation and show that when using ramping input for eliciting complex spikes, the most efficient inhibition strategy is rather perisomatic synaptic locations and slow synaptic time constants (rise time constant, τ 1 = 10 ms and decay time constant, τ 2 = 80 ms). Proximal GABAergic inputs with fast time constants (τ 1 = 0.25 ms, τ 2 = 3.75 ms) can also block complex spikes, but to a lower extent (Figure 7a1-4). With fast inhibition, the probability is high that a complex spike develops shortly after the first spike was stopped (illustrated in Figure 7a2). This is similar to what has also been observed using fast synaptic inhibition to terminate dendritic NMDA spikes (Doron, Chindemi, Muller, Markram, & Segev, 2017;Du et al., 2017). If instead the inhibition is placed in the dendrites, either in locations with large (red traces) or small (dark traces) NMDA current, the efficacy goes down (Figure 7b1-2). This is specifically true if the location is distal (158 µm somatic distance) with a small NMDA transient (black trace). Close to the peak inhibition (at about −30 ms), there is also a trend with less complex spikes the closer the inhibition is to the soma (most visible with slow inhibition).
When using the clustered input and either on-site or somatic inhibition, the on-site strategy is slightly more effective because it can kill the dendritic plateau as shown earlier (Figure 7c). Comparing the reduction of complex spikes over the two activation protocols, it is clear that ramping input is harder to inhibit (Figure 7a-c).
The results above show that inhibition can stop complex spikes, directly if triggered shortly before the initiation of the complex spike, but also indirectly if triggered earlier. This indicates that small changes induced by inhibition propages over long time scales, changing the dynamics of the system. F I G U R E 6 Complex spikes in iSPN are also dependent on reduced sodium current. In A, the iSPN models are compared to the dSPN models and experimental data (Day et al., 2008;Planert et al., 2013). The iSPNs are on average more excitable and dSPNs both somatically (as shown in a 1 ) but also dendritically, as indicated by the ability of backpropagating action potentials (bAP) to trigger calcium (Ca) transients (a 2 ). In B, synaptic bombardment is used to investigate the ability of the model to trigger regular and complex spikes. b 1 shows example responses under the unmodulated control, and under cholinergic (ACh) and dopaminergic (DA) modulation. ACh increases while DA reduces spike probability (b 2 ). In b 3-5, the ACh modulation protocol is used, but the sodium channel is reduced (naf) and the NMDA channel is increased (NMDA). Under these conditions, complex spikes can be triggered in a small proportion of trials (b 4 ). The sodium and NMDA factors (f-naf and f-NMDA) in the trials triggering complex spikes are largely reduced and increased, respectively (b 5 ). Gray areas mark the factor ranges normally used

| DISCUSSION
In the hippocampus, complex spikes are observed with dendritic plateaus, that in turn are thought to trigger behavioral time scale plasticity and place cell formation (Bittner et al., 2015(Bittner et al., , 2017. The simulations run in this study suggest that similar complex spikes can also be triggered in SPNs under certain conditions. Two major components of these complex spikes were further identified: (a) NMDA spikes originating in the dendrites and (b) persistent activation of the sodium channel window current in the axo-somatic region. Neuromodulation was involved in both mechanisms as it increases the dendritic response to an incoming signal and decreases the transient part of the sodium current. The resulting signal was further reliably propagated into the dendrites in the form of elevated Ca levels, in contrast to the more transient backpropagating action potentials (Day et al., 2008). Due to the long duration of the elevated Ca level, it is likely that such a signal would be involved in plasticity (Evans & Blackwell, 2015;Jędrzejewska-Szmek, Damodaran, Dorman, & Blackwell, 2017). The duration of the complex spike could also be prolonged if sequential NMDA spikes were triggered in the dendrites, for example, assuming ramping synaptic activation. Perhaps, this is also how the plateau potentials in the hippocampus are formed; complex spikes triggered by individual NMDA spikes lead to plasticity so that next time similar stimuli are presented additional spikes are triggered and eventually enough spikes to cause a full blown distributed plateau potential, ensuring plasticity globally.
Complex spikes have not been reported experimentally in SPNs to date. One reason for this might be that both distributed synaptic activation and reduced sodium current were needed. To test for complex spikes experimentally, one should use techniques to reduce the sodium current, preferably in an awake in vivo setup. Reduced sodium has been reported following neuromodulation (Chen et al., 2008), but should also result from periods of sustained activation (Payandeh, Gamal El-Din, Scheuer, Zheng, & Catterall, 2012;Rudy, 1978). Under ex vivo conditions, one could perhaps try to alter the extracellular concentration of sodium and potassium, or to use blockers of them, to decrease the amplitude of the AP and the AHP.

| Neuromodulation
Acetylcholine pushes many ion channels in the direction of increased probability for complex spikes (pCS), see Table 2. F I G U R E 7 Inhibition of complex spikes. In A, the soma of the cell is inhibited with two slow or 20 fast GABAergic input during the ramping protocol (illustrated in a 1 ). Both types of input can stop complex spikes, but the short duration of the fast input increases the probability for complex spikes (pCS) at a later stage (illustrated in a 2-3 ). In a 3 , the timing of the input is shifted relative to the onset of the complex spike and the percentage of traces that still include complex spikes are calculated. In b, the inhibition is placed in the dendrites, specifically in dendrites with either low (dark traces) or high (red traces) NMDA current (b 1 ). In b 2 , the pCS is shown for both slow GABA and fast GABA (illustrated in a 2-3 ), together with an illustration of the setup illustrating the somatic distance of the inhibited dendrites. Only one of the five backgrounds is used in b 2 . In C, the same is repeated for the clustered activation pattern. Here, the dendritic inhibition is placed in the same dendrite as the clustered stimuli, illustrated in the left panel. The right panel shows the pCS following the fast and slow GABA for somatic (black) and dendritic (gray) inhibition Potassium channels, like Kir, KCNQ, and Kaf, are reduced, which make the cell more excitable and reduce the afterhyperpolarization (AHP) of the action potential. The decreased AHP will likely reduce the proportion of sodium channels that are ready to fire when the threshold is reached again, and vice versa, reducing the spike amplitude and leading to increased pCS. Similarly, blocking of potassium channels in hippocampal CA1 neurons leads to spike patterns resembling the ones we see (Chen, Benninger, & Yaari, 2014).
Also, various calcium channels are modulated by ACh. The N-type calcium channel is reduced by ACh modulation (Howe & Surmeier, 1995;Perez-Rosello et al., 2005) which should further reduce the potassium current, since it is connected to the small conductance calcium-activated potassium channel (Sk). In agreement with a study in a non-neuronal cell (Pemberton & Jones, 1995), the L-type calcium channel has also been reported to be upregulated via an M4 dependent mechanism in dSPN (Hernández-Flores et al., 2015). This is perhaps controversial since DA is also reported to increase these channels in the same cell type (Galarraga, Hernández-López, Reyes, Barral, & Bargas, 1997;Surmeier, Bargas, Hemmings, Nairn, & Greengard, 1995). The intracellular effects of M4R and D1R activation are antagonistic in dSPNs (increased versus decreased cAMP-PKA production; see, e.g., Bruce et al., 2019;Nair, Gutierrez-Arenas, Eriksson, Vincent, & Hellgren Kotaleski, 2015). In this study, we simulated cholinergic modulation of the L-type channel with a decrease, since this has been reported following M1 activation (see Tables 2 and 3). It should, however, be mentioned that the general influence of the L-type calcium channels on excitability in the current models is low. L-type channel influence was limited not to affect NMDA spikes induction and shape; instead, T-type and R-type channels were upregulated in accordance with (Plotkin et al., 2011). We have further found no evidence currently in the literature for modulation of R-type and T-type channels by DA or ACh in dSPN (see table 2 and Lindroos et al., 2018).
Further, the persistent sodium current (Nap) is increased by ACh (Carrillo-Reid et al., 2009;Figueroa et al., 2002;Galarraga et al., 1999), while activation of the D1 type receptor decreases the peak of the transient sodium current in striatum (Schiffmann, Lledo, & Vincent, 1995), as well as in hippocampus (Cantrell, Smith, Goldin, Scheuer, & Catterall, 1997). Both these changes should further increase the pCS according to our simulations in dSPNs. Nap has also been connected to place cell function in the hippocampus (Hsu, Zhao, Milstein, & Spruston, 2018).
The sodium model used in this study is a compound model, including currents from multiple subcomponents and including a significant window current (Ogata & Tatebayashi, 1990). The major components of the sodium current should be the sodium channels Nav1.1, Nav1.2, Nav1.3, and Nav1.6 (see Table S1 in Hjorth et al., 2020;Muñoz-Manchado et al., 2018). Of these, the Nav1.6 channel carries a larger extent of the persistent sodium current (Y. Chen et al., 2008;Rush, Dib-Hajj, & Waxman, 2005). Additionally, the Nav1.6 channel expression is high in the AIS and less prone to neuromodulation (Chen et al., 2008). This hence gives support to our standard modulation protocol for DA (not modulating the axon initial segment).
The reduced sodium current following dopaminergic modulation of the dSPN has further been proposed to act as a safety mechanism, guarding against over-excitability . We here found that decreased transient sodium leads to increased pCS, so perhaps this reduction is also facilitating plasticity?
Synaptic currents are also modulated by DA and ACh. Here, we did not include the large increase of the NMDA current reported following cholinergic modulation (Calabresi, Centonze, Gubellini, Pisani, & Bernardi, 1998), since such an increase would cause overexcitation by itself (under our stimulation protocol) and a more recent study report no modulation (Higley, Soler-Llavina, & Sabatini, 2009). A change in this direction would, however, further increase the pCS since NMDA spikes are one component of complex spikes. Dopaminergic modulation of synaptic currents was included and contributed to increased pCS. It is also possible that direct and indirect modulation of GABAergic input could affect the pCS since inhibition was able to decrease pCS in our simulations.

| Crosstalk between dopamine and acetylcholine
Acetylcholine acts on targets in dSPNs via primarily muscarinic type 1 and 4 receptors (M1R and M4R), where the M4R is antagonistically coupled to the dopamine cascade. In this study, we did not investigate possible crosstalk between the two pathways (at pre-or postsynaptic sites). It is hence possible that some of the effects described here under concurrent cholinergic and dopaminergic modulation would not happen during physiological conditions. Most of the muscarinic effects reported here, however, belong to the M1R pathway and should not interfere directly with the dopaminergic effects, at least not postsynaptically, see Table 2. Time scale differences were neither investigated, but we tested also the effects of DA and ACh alone and in general then the effects were smaller.
Another question is if DA and ACh would be simultaneously high and how simultaneous activation would affect intracellular modulation pathways controlling synaptic plasticity? Recent studies suggest that coordinated ChIN activity can trigger DA release directly from DAergic axons (Threlfell et al., 2012). Further, using fast measuring techniques (Patriarchi et al., 2018), it was recently shown that | 2131 LINDROOS aND HELLGREN KOTaLESKI dopamine concentration is ramping up on a subsecond time scale in expectancy of a reward and that this increase is independent of the firing rate of the dopaminergic cell bodies (Mohebi et al., 2019). The decoupling of concentration and cell body firing also indicates a local regulation of DA concentration. The above stated DA release following coordinated ChIN activity makes them the most likely candidate (Threlfell et al., 2012). Concurrent DA and ACh signals in the dorsal striatum are also associated with locomotion induction (Howe et al., 2019).
Despite this recent evidence for concurrent DA and ACh transients, it is still unclear what the function of these concurrent transients is and how they influence cellular behavior of the two SPN subtypes. Sequential activation of ACh and DA is also indicated by LTP protocols where a ChIN burst-pause activation may be required in conjunction with DA release in dSPNs (Bruce et al., 2019;Nair et al., 2015). It is therefore clear that more studies are needed using minimal invasive patching techniques, preserving the intracellular signaling pathways (Lahiri & Bevan, 2020).

| Striatal up-states, dendritic plateau potentials and complex spikes
Under in vivo conditions, preferentially during sleep or anesthesia, the membrane potentials of SPNs generally fluctuate between two states, a silent, hyperpolarized state referred to as the down-state and an active depolarized state referred to as the up-state (see e.g., Wilson & Kawaguchi, 1996). These two states represent periods with low and high glutamatergic drive, respectively. Similarly to the up-states, our complex spikes also require high glutamatergic drive to form, but in contrast to the up-state, an additional decrease of the transient sodium current was also required. It is therefore possible that these two phenomena represent different operational modes.
In relation to striatal up-states, it has also been proposed that spatio-temporal clustered activation of synapses, which increase the probability of dendritic plateau potentials, could reduce the number of inputs needed to trigger up-states (Plotkin et al., 2011). In contrast to this, our simulations predict a higher proportion of complex spikes using a ramping synaptic input than using a brief clustered input. This may seem contradictory but can be explained by the difference in total activation duration. While the patterned activation (including spatio-temporal clustered activation) represented a brief increase in global input frequency of 1,000 Hz for 16 ms, the ramping input gave a more continuous increase of about 700 Hz when fully activated. The ramping input was also able to continuously trigger small NMDA-spike-like events in the dendrites, sometimes large enough to trigger complex spikes.

| Technical aspects
We here reported a crosstalk between the dendrites and the soma in the sense that excitation in a relatively small part of the dendrites is getting amplified by an axo-somatic component leading to cellwide depolarization and calcium influx. The widespread depolarization is caused by the release of potential locked up in the NMDA receptor, by reducing the magnesium block. Since the normal direction of information flow in neurons is from the dendrites to the soma, this mechanism is a more robust form of backpropagation than the more transient backpropagation following action potentials. Perhaps, this type of mechanism could sometimes function as the hypothesized neural correlate of the backprop algorithm used in artificial neural networks (Rumelhart, Hinton, & Williams, 1986)?