The microcircuits of striatum in silico

Significance Our aim is to reconstruct a full-scale mouse striatal cellular level model to provide a framework to integrate and interpret striatal data. We represent the main striatal neuronal subtypes, the two types of projection neurons (dSPNs and iSPNs) giving rise to the direct and indirect pathways, the fast-spiking interneurons, the low threshold spiking interneurons, and the cholinergic interneurons as detailed compartmental models, with properties close to their biological counterparts. Both intrastriatal and afferent synaptic inputs (cortex, thalamus, dopamine system) are optimized against existing data, including short-term plasticity. This model platform will be used to generate new hypotheses on striatal function or network dynamic phenomena.

injected and 10 min after the injection the pipette was slowly retracted from the brain. Postoperative analgesics were given (Buprenorphine, 0.08 mg/kg, i.p.).

Ex vivo experiments of single neurons
P30 -P60 days old mice were deeply anaesthetized with isoflurane and decapitated. The brain was removed and immediately immersed in ice-cold cutting solution containing 205 mM sucrose, 10 mM glucose, 25 mM NaHCO 3 , 2.5 mM KCl, 1.25 mM NaH 2 PO 4 , 0.5 mM CaCl 2 and 7.5 mM MgCl 2 . Parasagittal brain slices (thickness 250 µm) were cut with a Leica VT 1000S or VT1200S vibratome and incubated for 30-60 min at 35 °C in a submerged chamber filled with artificial cerebrospinal fluid (ACSF) saturated with 95% oxygen and 5% carbon dioxide. The ACSF contained 125 mM NaCl, 25 mM glucose, 25 mM NaHCO 3 , 2.5 mM KCl, 2 mM CaCl 2 , 1.25 mM NaH 2 PO 4 , 1 mM MgCl 2 . Subsequently, slices were kept for at least 30 min at room temperature before recording.
Whole-cell patch clamp recordings were obtained in oxygenated ACSF at 35°C. Neurons were visualized using infrared differential interference contrast (IR-DIC) microscopy [Zeiss FS Axioskop, Oberkochen, Germany; BX51WI (Olympus, Japan)]. tdTomato-expressing neurons were identified by switching from infrared to epifluorescence mode using a mercury lamp (X-Cite, 120Q, Lumen Dynamics) or epi-fluorescence LED source (Mightex bioLED, USA). Patch pipettes (borosilicate, Hilgenberg) were pulled with a Flaming / Brown micropipette puller P-1000 (Sutter Instruments) and had resistances of 5-9 MOhm. Recordings were performed in current clamp mode with an intracellular solution composed of either 130 mM K-gluconate, 5 mM KCl, 10 mM HEPES, 4 mM Mg-ATP, 0.3 mM GTP, 10 mM Na 2 -phosphocreatine or 105 mM K-gluconate, 30 mM KCl, 10 mM HEPES, 4 mM Mg-ATP, 0.3 mM GTP, 10 mM Na 2phosphocreatine (pH 7.25, 285 mOsm). In some experiments 0.3% Neurobiotin was added to the intracellular solution for subsequent reconstructions of the morphology of the recorded neurons (Vector laboratories, CA). Signals were amplified using a MultiClamp 700B amplifier (Molecular Devices, CA, USA), filtered at 5 kHz, digitized at 10 -20 kHz using ITC-18 (HEKA Elektronik, InstruTECH, NY, USA), and acquired using custom-made routines running on Igor Pro (Wavemetrics, OR, USA). Liquid junction potential was not corrected for (estimated at 9.5 mV, it is subtracted during feature extraction and model fitting). Throughout all recordings pipette capacitance and access resistance were compensated for and data was discarded when access resistance increased beyond 30 MOhm. The intrinsic properties of the different neuron types were determined by a series of hyperpolarizing and depolarizing current steps and ramps, enabling the extraction of sub-and suprathreshold properties. Based on electrophysiological features such as resting membrane potential, input resistance, membrane time constant (τ), sag, and firing frequency, SPNs, LTSs, FSs, and ChINs were in some cases also identified in the absence of fluorescent marker proteins.

Data on cortical and thalamic inputs
Optogenetic activation of striatal inputs was done after at least 5 min of bath application of 10 μM gabazine (SR-95531, Sigma-Aldrich). The stimulation was generated by a 1-Watt blue LED (wavelength 465 nm) and delivered through the 64x objective. Duration and intensity of the stimulation was controlled by an LED driver (Mightex Systems) connected to the ITC-18 acquisition board. The photostimulation diameter was approximately 400 μm with an illumination intensity of 0.4 mW / mm 2 . Light pulses of 2 ms duration were applied for activating cortical and thalamic inputs and stimuli were repeated for at least 8 sweeps at 0.1 Hz. Light intensity was reduced if necessary to ensure subthreshold EPSPs of at most 15 mV. In these recordings, all neurons were held at -75 + 2 mV to ensure that the driving force is equal for all cells when studying the synaptic strength and short-term plasticity.
Morphological reconstruction and optimization of electrophysiological models for the neurons Morphology All morphological reconstructions were digitally processed to ensure uniform quality and appearance using a set of morphology processing and repair scripts (https://github.com/a1eko/morphon). Single-point specification for somata was used, with center at the soma center and radius of the effective sphere. Correction for tissue shrinkage in zdirection was applied by expanding reconstruction progressively with the slice depth so that the lowest points remain unchanged and the points on the slice surface are shifted up the most (25%, if not reported otherwise). Axons were scaled up uniformly in (x, y)-plane by 10% and the dendrites were stretched along the principal axis preserving the length, which decreased their contraction and increased the reach by 10-15% [similar to unraveling in (2)]. Axons and dendrites cut at the borders of the slice were "repaired", i.e. extended, using intact branches of the same topological order from the inner part of the slice.
SPN. Digital reconstructions of SPNs from the dorsal striatum of the mice were obtained from a web-based inventory NeuroMorpho.Org, NMO (http://neuromorpho.org), see (3). Reconstructions were assigned to dSPN (n=4) and iSPN (n=4) types based on their morphometric features (see Tables S2, S3). Since no local axon collaterals were available for the selected SPN cells, axons of cortical GABAergic FS cells were instead used as in Humphries et al. (4). Selected synthesized morphologies for SPNs are available via DOI:10.25493/YP8N-HD2.
FS. Morphological reconstruction of the striatal FS cells (n=3) were selected from a large pool of the cortical PV-expressing FS interneurons at NMO and validated against striatal FS morphometry (see Table S4; selected morphologies DOI:10.25493/3FNQ-5KG).
LTS. The cells were filled with Neurobiotin during patch clamp recordings. Following experiments, slices were transferred to 4% paraformaldehyde solution containing 14% picric acid in 0.01 M PBS for 12h at 4°C for fixation. Cells filled with 0.3% Neurobiotin were visualized by rinsing the slices in PBS followed by incubation in Cy5-conjugated streptavidin antibody (1:1000, Jackson ImmunoResearch Laboratories) diluted in 0.01 M PBS containing 0.3% Triton-X100 and 1% BSA for at least 6h. The slices were imaged using Confocal (ZEISS LSM 800) and z-stacks were retrieved. The confocal z-stacks were used in a semi-manual reconstruction using Neutube (5) and custom code. The most complete reconstruction (n=1) was selected for optimization (DOI:10.25493/DVPH-RDE).
ChIN. The slices were fixated as described above and processed with the Vectastain EliteABC kit (Vector Laboratories). Slices were thoroughly washed in PBS and incubated in 0.6% hydrogen peroxidase in methanol for 20 min, rinsed in PBS and transferred to the ABC solution for 3 hr. After rinsing in PBS, slices were incubated in diaminobenzidine (DAB; ImmPACT DAB, Vector Laboratories) for 3-5 min, rinsed, and dehydrated in alcohol prior to mounting in Entellan (Merck). The neurons were manually reconstructed (n=10; DOI:10.25493/3EV4-TDG, DOI:10.25493/ADRK-VJP) using Neurolucida (MBF Bioscience) coupled to Zeiss Axio Imager.A1. Also, here, the most complete reconstruction was selected for optimization.

Optimization of multi-compartmental neurons
Biophysically detailed single-cell models were optimized using BluePyOpt (6). Electrophysiological features for the optimizer were extracted from the experimental traces in IGOR Pro binary formats using the Electrophys Feature Extraction Library, eFEL (https://github.com/BlueBrain/eFEL). The optimization procedure was divided into several stages, in order to fit parameters of the passive membrane, subthreshold and higher-voltage activated ion channels.
Optimizer aims to match multiple features of electrophysiological traces simultaneously (e.g. Figures 2B and 3B) and returns several parameter sets corresponding to acceptable models, where the features were less than 2 standard deviations from the population average.

SPN.
The models of the SPN are based on Lindroos et al. (7) and previous studies (8)(9)(10)(11). They use detailed reconstructed neuronal morphology and active ion channels distributed over the soma, dendrites and axon initial segment. Choice of the ion channels reflects expression of the channel markers (RNA-seq) shown in Fig. S1. Inward-rectifying potassium channel Kir 2.3 ensures low resting membrane potential as well as current rectification characteristic to SPNs. Slowly inactivating A-type channel Kv 1.2 determines delay to the first spike, while the fast inactivating sodium Nav 1 and potassium Kv 4.2 channels shape the action potential. Currently only Nav 1 and Kv 1 are placed on the axon initial segment although other options are also possible (12). Multiple calcium channels of types L, N, P/Q, R and T (see Table S1) are implemented.
Calcium-dependent potassium channels KCa 1.1 (BK) and KCa 2.2 (SK) affect AHP and control the firing rate. Distribution of the ion channels in the SPN model is shown in Table S6.
Optimization of the electrophysiological model is performed in two stages. Initially, only passive membrane parameters (membrane conductance, axial membrane resistance and membrane leak reversal potential) are estimated together with conductances of the potassium channels for the subthreshold stimuli. Obtained optimal values of the axial membrane resistance and the leak reversal potential are fixed and transferred to the next optimization stage unchanged. Values of the membrane conductance and conductance of Kir channel are used as reference and varied within two orders of magnitude at the final stage. Each model is then validated with respect to the subthreshold response, discharge rates, plateau potential and dendritic calcium elevation to select the best individuals (as in Fig. 2C-E and Fig. S3C-E). Resulting models based on the morphology of four dSPNs and four iSPNs, respectively, are fitted to electrophysiological recordings identified via DOI:10.25493/MZE0-BH5.
The calcium response to a backpropagating action potential was tested by triggering one action potential in the soma, using a short duration, high amplitude current injection (2 ms, 2.5 nA), and measuring the calcium response in the dendrites. The signals were normalized to the response in the proximal part (30-40 μm for dSPN and 40-50 μm for iSPN) and plotted in relation to somatic distance. More proximal measurements were discarded. The measurements were fitted to an exponential curve and compared to experiments (13), see Fig. 2E and Fig. S3E.
SPNs can trigger dendritic NMDA elicited plateau potentials (11,14). In the models these events can be induced using sequential activation (inter activation interval 1 ms) of 16 glutamatergic synapses in a distal part of the dendrite (at 90-120 μm somatic distance). To validate dendritic plateaus in the used SPNs, maximal conductance of single NMDA components was set to 1880 pS and AMPA/NMDA ratio to 1/5 for dSPN (11). For iSPN the maximal NMDA conductance was reduced to 1400 pS to avoid spiking. Following Du et al. (11), a constant somatic current injection was given throughout the simulation to depolarize the cell to about -75 mV.
FS. Electrophysiology of the fast-spiking interneurons is characterized by the short AP and deep AHP attributed to the fast delayed-rectifying Kv 3 channels [ (15,16]; Fig. S1). Inactivating potassium channels Kv 1.2 , Kv 4.2 and inwardly-rectifying Kir 2.3 are also included in the model to account for occasional delayed first spike and a degree of subthreshold rectification displayed by the recorded cells. P/Q-type calcium channel Cav 2.1 and BK-type calcium-dependent potassium channel KCa 1.1 are included as well to match the channel expression pattern (see Fig. S1). Uniform distribution of the ion channels in the soma and dendrites is assumed. The FS models were optimized in a single pass; see DOI:10.25493/E883-NFA.
LTS. Due to the significant spherical asymmetry of LTS cells and lack of the morphometric statistics, repair of cut dendrites was not attempted and the best available reconstruction with its electrophysiology was used for single-cell model optimizations (DOI:10.25493/5GE0-6MF). We restricted the model to have ion channels of the subtypes most expressed. Among these are the T-type calcium channels, slow delayed rectifier potassium channels and M-type channels.
ChIN. The ion channels were selected based on RNA-seq data (Fig. S1) and Maurice et al. (17), and their distribution followed previous publications (18,19). The dendrites contain models for Kir, HCN, Kv 4 , BK, SK channels and L-type and P-type calcium channels. The ion channels in the soma include Kv 2 , KCNQ and sodium channels while the HCN and Kir channels were excluded. The final models for the selected morphology were validated as previously described and verified to reproduce a spontaneous activity with minimal stimulation, rebound response and pause response (DOI:10.25493/3NTS-Q0B).

Dopamine modulation of ion channels and receptors
Dopamine modulation was implemented as a change of maximal conductances as per the current knowledge in the field [see Table S7-S10 and Lindroos et al. (7)]. Since the experimental conditions, such as model species, drugs applied and experimental protocols differ between studies, we randomly drew modulation factors from within reported ranges following Lindroos et al. (7). The resulting set of factors was then applied to ion channels during simulation of selected ex vivo protocols for each cell type. The sets that reproduce behavior in accordance with the literature (Table S7-S10) were accepted and stored for network simulation, the rest were rejected. For SPN and FS the protocols used for optimization were also used to assess correct modulation. Here quantified as a change of excitability. The iSPN excitability is decreased by dopamine while for dSPN and FS it is increased (20)(21)(22). For ChINs we quantified the modulation as the relative increase in time to spike following a burst-pause protocol (23). The protocol here consisted of a stimulation for 800 ms using a depolarizing current injection (50 pA). For LTS the effects of dopamine are not well studied. Here we randomly modulated the sodium and HCN channels within a range of +/-20%, with the goal to induce spiking from a hyperpolarized state (24). The factors that gave this behavior (both positive) were opposite to the ones giving correct results for ChINs. This is in line with the different dopamine receptor subtypes responsible for the modulation in LTS (D1R-type) and ChINs (D2R-type), see Tables S9-S10. The protocol used for LTS consisted of a constant negative current (-16.5 pA), hyperpolarizing the membrane to about -60 mV.
Uniform modulations over cellular compartments were used in all cells, except for the dSPN where the axon initial segment was not modulated. Non-uniform modulation of dSPN has been proposed to give increased excitability (7).
During simulation individual cells were randomly paired with one set of factors. The modulation itself was implemented in two ways, static and transient. In the transient protocol channels were gradually updated from zero modulation up to the modulated maximum. This was used for network simulation, although the maximal effects were set to 50% in Fig. 10. The static protocol used the modulated maximum throughout, and was used to assess the effect in single cells (except for LTS, where only transient activation was used). Changes were set to have immediate effect in both protocols.

Placement of the neurons and synapse prediction in silico
The Snudda software is used for the steps described below.

Initialization of the network
The network initialization defines all subcomponents of the network in a configuration file. Specifically, the network volume, cell models and synaptic dynamics are specified. The user specifies multiple morphological reconstructions for each neuron type, and a family of parameter sets for each reconstruction generated by the BluePyOpt optimiser as described above. In the case of the striatal network the initialization procedure has been streamlined, and the user only needs to specify the number of neurons and volume to use, i.e. the full striatum, a cube or a slice preparation. The configuration file is then defined with the requested number of neurons with the correct proportion of the different neuron types. For other brain structures the user has to generate their own configuration file.

Placement of neurons
Striatal neurons are placed within a mesh downloaded from the Allen Brain Atlas with a volume of 21.5 mm 3 [see also (25)]. The neurons are placed sequentially by drawing coordinates uniformly within the striatal volume. If there are previously placed neurons within an exclusion zone with a radius of 15 μm, then the position is rejected. Only positions within the striatal volume count towards the total number of neurons, however, positions outside are tracked to avoid artificial inflation of the density at the border. The placed neurons are then rotated according to the specification in the configuration file, the default is random rotation. The average neuron density within the volume is uniform, but future versions can be extended to handle density gradients.

Connecting the microcircuit
The connectivity of the microcircuit is derived from the morphology of the neurons placed in the volume. Only axons and dendrites in close proximity to each other can form putative synapses, analogously dendrites need to be in close apposition to form gap junctions. The touch detection algorithm (Fig. 7) divides the space into 3x3x3 μm voxels, and if two neurites occupy the same voxel a connection between them is possible. To be able to distribute the work efficiently 100x100x100 voxels are grouped into a hypervoxel and are processed together. The touch detection uses ipyparallel to distribute the workload between the processors.
The algorithm has three steps: 1) Identify which neurons are present in each hypervoxel. 2) Placing putative synapses in the hypervoxels by first marking the voxels that axons and dendrites of different neurons occupy, and then detecting when they occupy the same voxel, i.e. in close apposition to each other. 3) Pruning the putative synapses based on a set of rules. Below each of the four steps is discussed in more detail.
Step 1: Setting up the touch detection. For each neuron, the algorithm calculates which hypervoxels it occupies, by checking every point in the reconstructed neuron. Each hypervoxel thus has a list of which neurons it contains. The hypervoxels are then sorted in descending order based on the number of neurons it has, so that the hypervoxel with most neurons is processed first. This creates an even distribution of the workload, making sure all the processors finish at around the same time. Typically, a neuron will span several hypervoxels, and the workload increases proportionally to the number of hypervoxels the neurons occupy. Increasing the size of the hypervoxels will reduce the total amount of processing needed to be done, but would increase the memory requirements.
Step 2: Detecting putative synapses. A line drawing algorithm traverses each part of the neuron in voxel sized steps x v = x 0 + k*dx. Here x 0 , x 1 are the endpoints of the compartment, dx = (x 1x 0 )/n is the step size along the direction of the compartment, and n is the number of voxel steps along x, y or z dimension (whichever is largest). The algorithm also handles marking the voxels belonging to the somas. Putative synapses are then detected by finding which voxels are occupied by axons and dendrites of different neurons. If connections between the neurons are allowed by the rules defined in the configuration file, then a putative synapse is added. For each pair of neurons only one synapse is allowed per voxel, however multiple neuron pairs can have synapses in the same voxel. The conductances of the synapses are drawn from a normal distribution.
Step 3: Pruning of synapses to match experimental data. There are many more putative locations for synapses than real synapses. The pruning consists of randomly 1) removing a fraction of all synapses, 2) (optional) distance-dependent filter favoring either synapses proximally or distally, 3) removal of synapses between pairs with excessive number of synapses 4) removal of neuron pairs with too few synapses, and 5) removing all synapses between a fraction of coupled pairs. Table S11 gives the parameters used for the pruning of the striatal network. These rules are defined to allow the researcher to shape the number of synapses connecting a pair of neurons and the connection probability between the pairs. First a set of randomly selected synapses are removed, leaving a fraction f 1 . This reduces the average number of synapses between coupled pairs, and also disconnects some neurons that were originally only connected by a few putative synapses. Some neuron types are known to synapse more proximally on their target neuron while others prefer to form synapses more distally. To account for this, distance dependent pruning has been implemented between LTS-SPN and FS-SPN pairs (equations given in the caption of Table  S11), and later also added for synapses originating from SPNs (Fig. 8, S4) to match experimental data (see Results). Neurons located close to each other can have an excessive number of putative synapses, and removing a part of those synapses by using f 1 would disconnect too many distal neurons. To combat this, we increase the probability of removing synapses if n > softMax, where softMax is a tuned parameter. Here the probability of keeping each synapse is P sm (n) =(2*softMax) / ((1+exp(-(n-softMax)/5))*n). This means that for each synapse past the softMax threshold the likelihood of keeping that synapse is reduced. We also implement a rule that if there are too few synapses between a pair, then they are decoupled with probability Pμ = 1.0/(1.0 + exp(-8.0/μ 2 *(n -μ 2 )), similarly to Markram et al. (2). All these rules shape the distribution of how many synapses connect coupled pairs. The last pruning rule randomly removes all synapses between coupled pairs, leaving a fraction a 3 of the pairs. This maintains the distribution of synapses between pairs of coupled neurons, while modifying the connection probability between pairs. Splitting the pruning into these different steps simplifies matching experimental data. The pruning parameters of the network were manually tuned (Table S11). We also use the same algorithm for gap junctions, with the exception that they connect dendrites together. In the striatum gap junctions are only present between dendrites, and synapses exist only between axons and dendrites or somas. Currently axon-axonal connections are not supported, but that functionality could easily be added using the same algorithm. Our approach is similar to the rules employed for the cortical microcircuit (2). SONATA (26) export functionality is being implemented for compatibility with other tools supplied by the Human Brain Project, Blue Brain Project and Allen Institute for Brain Sciences.

Density estimation when axon is missing
The touch detection algorithm uses the morphology to place putative synapses in the circuit, however the axonal arborization is not always available. In those cases the software allows the user to instead use a density estimation of the axonal arborization which is specified as a function of either radius from the center of the soma, or of (x, y, z)-coordinates. Using this scheme, a radial symmetry density was defined for ChIN axons and for LTS axons a sum of three coordinate-based functions was used. The magnitude of the density was set so that there would be enough synapses to support the connectivity probability seen in experiments. To verify the LTS axonal distribution a set of points was placed using the density, and compared to existing axonal Sholl analysis (27).

Analyzing the network connectivity
One of the strengths of a model is that we have direct access to the entire state of the system. To compare it to biological experiments we employ a set of virtual experiments in the model. For the microcircuit connectivity experiments a 1 mm 3 cube with 80 500 neurons were created. The pairwise connection probability as a function of distance ( Fig. 8A and Figs. S4A, S5A, S6A, S7A) is calculated, for each presynaptic neuron, by selecting a set of potential postsynaptic neurons, and checking if they were connected. The data is then binned based on the distance between preand postsynaptic neurons. Here we have restricted the calculations to about 10 million pairs with the presynaptic neuron in the center. The distribution of the number of synapses between connected pairs (Fig. 8B and Figs. S4B, S5B, S6B, S7B) and the number of connected neighbors ( Fig. 8C and Figs. S4C, S5C, S6C, S7C) are based on postsynaptic neurons within a cube of side 200 μm placed in the center of the 1 mm 3 modelled volume. This includes the majority of the neurons that can connect to each postsynaptic neuron analyzed. The cumulative distribution of synapse locations along the dendrites employs a different counting scheme where only the neurons in the corners are included, and each synapse is counted eight-fold to compensate for only counting part of the surrounding volume ( Fig. 8D and Figs. S4D, S5D, S6D, S7D). In addition, we also employ simulations with activation of either single neurons, or a set of neurons, to verify connection strength (described further below).
Constraining the synaptic properties and generating input for the model

Synaptic dynamics
The synaptic dynamics of the cortical and thalamic inputs are derived from experimental data (28). Using optogenetic techniques, axonal terminals from ipsilateral M1, S1, thalamus and contralateral M1 are selectively activated, while simultaneously recording SPN, FS, LTS and ChIN in the striatum (See Materials and Methods -data on cortical and thalamic input). Additionally, traces were recorded with a brief 2 ms hyperpolarizing delta pulse, from which the membrane time constant and input resistance were extracted. The data traces are available in the Knowledge Graph in EBRAINS, and the resulting families of synapse models are included in the GitHub repository. Custom written code based on NEURON and PySwarm was used to fit the glutamate receptor dynamics to the amplitudes of postsynaptic potentials (PSP) using a minimal neuron model having a soma with membrane time constant and input resistance similar to each particular experiment. The amplitude of each PSP was based on an extrapolation of the decay following the previous PSP, and the error was the absolute difference in amplitudes between the experimental data and the model trace. Here the errors of the first and last PSPs were weighted three times higher in the error function. A term was also added for the error in the decay of the recovery PSP, to better capture the AMPA/NMDA ratio. Example synaptic dynamics are shown in Fig. 9A and Fig. S8. For the intrastriatal synapses we created a surrogate dataset from Table 1 in Planert et al. (20). Here, the mean and standard deviation of the parameters u, tauF and tauD of a Tsodyks-Markram model (29) generating the peak amplitudes that match experiments were given. To fit parameters for a dynamic Tsodyks-Markram model we generated a set of surrogate data sampled from the parameters of the original model. The table in Planert et al. (20) also specified the pair pulse ratio (PPR, second to first peak), as well as the recovery test ratio (RTR, last to first peak). Using the newly generated parameter set we calculated the sequence of peaks from the original model, and iteratively excluded traces so that the population of parameters would give the same mean and standard deviation for the PPR and RTR listed in the paper. The reduced parameter set that passed this validation was then used to generate a surrogate data set used to fit the dynamic Tsodyks-Markam synaptic model for the GABAergic synapses similar to those recorded in Planert et al. (20).
Synaptic conductances for SPN-SPN synapses were drawn from a normal distribution with mean 240 pS and standard deviation 100 pS, where we assume single channel conductance to be 24 pS (30) and each synapse having 10 channels. With an average of four synapses between SPN-SPN we match the total conductance measured between neurons in Taverna et al. (30). A simulation of a virtual slice preparation was then set up with chloride reversal potential -40 mV and voltage clamp at -80mV to match Planert et al. (20), and failure rates from Taverna et al. (30). Pairwise neurons were checked and it was found that in order to match the depolarization in experiments we needed to remove proximal synapses between SPN-SPN (See Results, Fig. 8D i-iii and Fig.  S4D i-iii). To validate the SPN-SPN connection strength we then replicated the experiment of Chuhma et al. (31) where 10% of SPN population are activated. The response in non-activated SPN, ChIN and FS were measured. To match those experiments a chloride reversal potential of -2mV and voltage clamp at -60mV was used. The SPN-ChIN data from that experiment was not used for validation as it had been used for tuning.
Szydlowski et al. (32) reported that 2/12 FS-LTS pairs were connected, and showed synaptic potentials of 0.5 and 0.8mV. While replicating this as a virtual experiment with calculated chloride reversal potential of -39mV and a holding voltage of -76mV it became apparent that the FS-LTS synapses had to have more pruning than FS-SPN synapses, and that we also had to weaken the conductance for FS-LTS synapses by a factor of 10 as compared to FS-SPN synapses to match the recorded voltage deflection.
In Straub et al. (33) FS and LTS populations were independently stimulated optogenetically, and the response was recorded in dSPN, iSPN and ChIN. Virtual experiments with estimated chloride reversal potential of 0 mV from the chloride concentrations in the recording pipette and holding voltage of -70mV were performed. We had previously estimated the LTS axon shape and density from Ibáñez-Sandoval et al. (27) and these new virtual experiments led us to an increase of the connection probability between LTS-SPN and a scale-up of synapse conductances to match the current response seen in SPNs. The FS to SPN response matched without any additional tuning. Distance dependent synaptic delay was introduced into the model by assuming an axonal propagation speed of 0.8 m/s (34,35) and a synaptic transmission delay of 1 ms.

Optimizing synapse dynamics, strength and activation frequency
To generate correlated input spikes for the microcircuit a shared set of input spikes is first generated with frequency f. Spikes from this shared set propagate into the individual spike trains with a probability p. Additional independent spikes are then added with a frequency (1-p)*f Hz to generate the final correlated spike trains with a total average frequency of f Hz (36). Additionally, the Snudda framework allows the spikes to be jittered.
The synaptic input is specified in a JSON file. For each input type and neuron type the user specifies the parameters in the section above, as well as any additional jitter. The JSON file also specifies which NEURON MOD file is used for the synapse, as well as any file specifying additional synapse parameters. Optionally the user can specify a CSV file where the explicit spike times could be loaded from instead of pre-generated by Snudda. In addition, it is possible to specify if the synapses should be uniformly distributed on the dendrites, or if there should be a distance dependence. Here cortical inputs are set based on the density function: P c = 0.05/(1+exp(-(d-30)/5)) while thalamic inputs are set using: P t = 0.05*exp(-d/200), where d is the arc distance from the soma (in μm).
It is also possible to include reconstructed long-distance axons from cortex and thalamus, and activate these virtual neurons with specific spike trains to drive the network. Here, the location of synapses onto the dendrites in the striatum is then determined using the touch detection and pruning rules.

Simulating the microcircuit
The microcircuit can be simulated with Parallel Neuron through the Snudda interface. Using a supercomputer like the Cray XC40 it is possible to simulate large networks. For the dopamine network we simulated 10,000 detailed neurons. Here the neurons are distributed in a round robin fashion between the available processors. The code then sets up the connectivity based on the synapses specified in the configuration file, external synaptic input from cortex and thalamus are read from the input file and connected to the neurons. Extrapolating from in vivo recordings in mice, we estimated the baseline spontaneous activity of cortical neurons to be in the range of 1-5 Hz (37). The baseline thalamic activity was estimated to be 1-2 Hz (38). When simulated with these values for the cortical and thalamic inputs, the average activity of SPN, FS and ChIN in our striatal microcircuit simulation was similar to frequencies recorded in vivo (39). By default, the simulation will save the spike times for all neurons, optionally the membrane potential of some or all of the neurons can also be saved. Scripts are included on GitHub for running both on a local machine and on supercomputers using SLURM.
How to use the software The neuron population and the volume they occupy, as well as rules for which neurons can connect to each other, and pruning rules are specified using a network JSON file.
To generate an example network with 5000 striatal neurons: snudda init networks/helloworld --size 5000 snudda detect networks/helloworld snudda prune networks/helloworld snudda input networks/helloworld --input config/input-config.json snudda simulate networks/helloworld The first line creates the JSON network config file that is used to define the network. The next line runs the touch detection to locate putative synapses, followed by the command that does the pruning of the synapses. After that the network synaptic input is generated, and finally the simulation is run. All files are stored in the networks/helloworld directory. By default, only the spike times are saved from the simulation, but by using the --voltOut myVoltfile.csv option somatic voltages can be saved. Fig. S1. Expression of ion channels in the main neuron phenotypes from the dorsal striatum of the mouse derived from dataset A in (40), Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo), accession number GSE97478. Color density represents molecular counts from low expression (light) to high expression (dark) in logarithmic scale. The database contains records for 694 neurons, identified as cells expressing Fox3 (NeuN), Eno2 (NSE) or Map2. Neuron types are selected using expressions of Drd1a for dSPN, Adora2a for iSPN, Pvalb for FS, Npy/Sst for LTS, and Chat for ChIN. Molecular markers used for selection of the ion channels are listed in Table S1. The channel expression pattern is consistent between multiple studies available for mouse striatum (40)(41)(42).  supra-threshold responses to current injections for a model neuron (black) and the corresponding experimental data (red). An example model fit to experimental data, with the current protocol used. (C) Population behavior for models and experiment: voltage-current and frequency-current relations shown for 4 iSPN models optimized to corresponding data. (D) Somatic potential response to spatio-temporal clustered synaptic input, demonstrating the model's ability to trigger NMDA dependent plateau potentials. Experimental data (red) digitized from Du et al. (11). (E) Normalized change in calcium concentration in response to a backpropagating action potential (triggered with a short duration 2 ms, high amplitude 2.0 nA current injection). Experimental data digitized from Day et al. (13). Model data in black, experimental data in red. Here we show statistics for neurons in the center of the volume to avoid edge effects. Note that the multimodal distribution seen here is a consequence of only using a limited number of reconstructions for dSPN and iSPN. Preliminary modelling shows that adding a larger number of reconstructions creates a unimodal distribution, however, currently we only have optimized models for the morphology of four iSPNs. Future versions will include more reconstructions. (C) (i,ii) the black line shows the distribution obtained for a bigger set of reconstructions (n=100,000) using a jitter to promote morphological variability (each reconstruction scaled randomly in 3D by a factor (0.8-1.2); and in addition, 20 randomly selected dendritic branches are rotated to a random angle around its principal axis).       (20) and GABAergic Tsodyks-Markram model (black) fitted using a single compartment model (described in Materials and Methods). Protocol includes eight pulses at 20 Hz followed by a recovery pulse to probe the synaptic dynamics. Table S1. List of genes and related ion channels. Adapted from the Mouse Genome Informatics resource (http://www.informatics.jax.org).