Striatal D1R medium spiny neuron, including a subcellular DA cascade (Lindroos et al 2018)

 Download zip file 
Help downloading and running models
We are investigating how dopaminergic modulation of single channels can be combined to make the D1R possitive MSN more excitable. We also connect multiple channels to substrates of a dopamine induced subcellular cascade to highlight that the classical pathway is too slow to explain DA induced kinetics in the subsecond range (Howe and Dombeck, 2016. doi: 10.1038/nature18942)
1 . Lindroos R, Dorst MC, Du K, Filipovic M, Keller D, Ketzef M, Kozlov AK, Kumar A, Lindahl M, Nair AG, Pérez-Fernández J, Grillner S, Silberberg G, Hellgren Kotaleski J (2018) Basal Ganglia Neuromodulation Over Multiple Temporal and Structural Scales-Simulations of Direct Pathway MSNs Investigate the Fast Onset of Dopaminergic Effects and Predict the Role of Kv4.2. Front Neural Circuits 12:3 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Axon; Channel/Receptor; Dendrite; Molecular Network; Synapse; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Basal ganglia; Striatum;
Cell Type(s): Neostriatum medium spiny direct pathway GABA cell; Neostriatum spiny neuron;
Channel(s): I A; I A, slow; I Calcium; I CAN; I K; I K,Ca; I K,leak; I Krp; I Na,t; I Potassium; I R; I T low threshold; Kir;
Gap Junctions:
Receptor(s): D1; Dopaminergic Receptor; AMPA; Gaba; NMDA;
Transmitter(s): Dopamine; Gaba; Glutamate;
Simulation Environment: NEURON; Python;
Model Concept(s): Action Potentials; Detailed Neuronal Models; Electrical-chemical; G-protein coupled; Membrane Properties; Neuromodulation; Multiscale; Synaptic noise;
Implementer(s): Lindroos, Robert [robert.lindroos at]; Du, Kai [kai.du at]; Keller, Daniel ; Kozlov, Alexander [akozlov at];
Search NeuronDB for information about:  Neostriatum medium spiny direct pathway GABA cell; D1; AMPA; NMDA; Gaba; Dopaminergic Receptor; I Na,t; I T low threshold; I A; I K; I K,leak; I K,Ca; I CAN; I Calcium; I Potassium; I A, slow; I Krp; I R; Kir; Dopamine; Gaba; Glutamate;
MSN model used in Lindroos et al., (2018). Frontiers

Robert Lindroos (RL) <robert.lindroos at>
The MSN class and most channels were implemented by 
Alexander Kozlov <akozlov at>
with updates by RL

Implemented in colaboration with Kai Du <kai.du at>

from __future__ import print_function, division
from neuron import h
from joblib import Parallel, delayed
import multiprocessing
import numpy                as np
import matplotlib.pyplot    as plt
import plot_functions       as fun
import MSN_builder          as build


def save_vector(t, v, outfile):
    with open(outfile, "w") as out:
        for time, y in zip(t, v):
            out.write("%g %g\n" % (time, y))                     

def main(   par="./params_dMSN.json",        \
                            sim='vm',       \
                            amp=0.265,      \
                            run=None,       \
                            simDur=1000,    \
                            stimDur=900     ): 
    # initiate cell
    cell    =   build.MSN(  params=par,                             \
                            morphology='latest_WT-P270-20-14ak.swc' )
    # set cascade--not activated in this script, 
    # but used for setting pointers needed in the channel mechnisms
    casc    =   h.D1_reduced_cascade2_0(0.5, sec=cell.soma) 
    # set pointer target in cascade
    pointer =   casc._ref_Target1p    
    # set edge of soma as reference for dendritic distance 
    h.distance(1, sec=h.soma[0])
    # set current injection
    stim        =   h.IClamp(0.5, sec=cell.soma)
    stim.amp    =   amp  
    stim.delay  =   100
    stim.dur    =   stimDur    
    # record vectors
    tm  = h.Vector()
    vm  = h.Vector()
    tstop       = simDur 
    # dt = default value; 0.025 ms (25 us)
    # set pointers; need since same mechanisms are used for dynamic modulation of channels.
    # Modulation of channels is not used in this script
    for sec in h.allsec():
        for seg in sec:
            # naf and kas is in all sections
            h.setpointer(pointer, 'pka', seg.kas )
            h.setpointer(pointer, 'pka', seg.naf )
            if'axon') < 0:    
                # these channels are not in the axon section
                h.setpointer(pointer, 'pka', seg.kaf )
                h.setpointer(pointer, 'pka', seg.cal12 )
                h.setpointer(pointer, 'pka', seg.cal13 )
                h.setpointer(pointer, 'pka', seg.kir )
                if'soma') >= 0:
                    # N-type Ca (can) is only distributed to the soma section
                    h.setpointer(pointer, 'pka', seg.can )
    # configure simulation to record from both calcium pools.
    # the concentration is here summed, instead of averaged. 
    # This doesn't matter for the validation fig, since relative concentration is reported.
    # For Fig 5B, where concentration is reported, this is fixed when plotting.
    # -> see the plot_Ca_updated function in plot_functions.
    if sim == 'ca':
        print('configure', sim, 'simulation')
        for i,sec in enumerate(h.allsec()):
            if'axon') < 0: # don't record in axon
                for j,seg in enumerate(sec):
                    sName ='[')[0]
                    # N, P/Q, R Ca pool
                    cmd = 'ca_%s%s_%s = h.Vector()' % (sName, str(i), str(j))
                    cmd = 'ca_%s%s_%s.record(seg._ref_cai)' % (sName, str(i), str(j))
                    # the L-type Ca
                    cmd = 'cal_%s%s_%s = h.Vector()' % (sName, str(i), str(j))
                    cmd = 'cal_%s%s_%s.record(seg._ref_cali)' % (sName, str(i), str(j))
                    # uncomment here if testing kaf blocking effect on bAP
                    #block_fraction = 0.2
                    #gbar           = seg.kaf.gbar
                    #seg.kaf.gbar   = (1 - block_fraction) * gbar
    # solver------------------------------------------------------------------------------            
    cvode = h.CVode()
    # run simulation
    while h.t < tstop:
    # save output ------------------------------------------------------------------------
    if sim == 'ca':
        print('saving', sim, 'simulation')
        # vm
        save_vector(tm, vm, ''.join(['Results/Ca/vm_', sim, '_', str(int(amp*1e3)), '.out']) )        
        # ca
        for i,sec in enumerate(h.allsec()):
            if'axon') < 0:
                for j,seg in enumerate(sec):
                    sName       ='[')[0]
                    vName       =   'ca_%s%s_%s'  %  ( sName, str(i), str(j)  )
                    v2Name      =   'cal_%s%s_%s' %  ( sName, str(i), str(j)  )
                    fName       =   'Results/Ca/ca_%s_%s.out'  %  ( str(int(np.round(h.distance(seg.x)))), vName )
                    cmd     = 'save_vector(tm, np.add(%s, %s), %s)' % (vName, v2Name, 'fName' ) # this is were concentrations are summed (see above)
    elif sim == 'vm':
        print('saving', sim, 'simulation', str(int(amp*1e3)))
        # vm
        save_vector(tm, vm, ''.join(['Results/FI/vm_', sim, '_', str(int(amp*1e3)), '.out']) )

# Start the simulation.
# Function needed for HBP compability  ===================================================
if __name__ == "__main__":
    print('starting sim')
    # dendritic validation: change in [Ca] following a bAP (validated against Day et al., 2008)
    current = 2000
    main( par="./params_dMSN.json",          \
                amp=current*1e-3,           \
                simDur=200,                 \
                stimDur=2,                  \
                sim='ca'                    )
    print('starting somatic excitability simulation')                                               
    # somatic excitability (validated against FI curves in Planert et al., 2013)  
    currents    = np.arange(-100,445,40)
    num_cores   = multiprocessing.cpu_count()
    Parallel(n_jobs=num_cores)(delayed(main)(   par="./params_dMSN.json",   \
                                                amp=current*1e-3,           \
                                                run=1,                      \
                                                simDur=1000,                \
                                                stimDur=900                 \
                        ) for current in currents)
    currents    = np.arange(320,445,40)
    Parallel(n_jobs=num_cores)(delayed(main)(   par="./params_dMSN.json",   \
                                                amp=current*1e-3,           \
                                                run=1,                      \
                                                simDur=1000,                \
                                                stimDur=900                 \
                        ) for current in currents)
    print('all simulations done! Now plotting')