Spinal Dorsal Horn Network Model (Medlock et al 2022)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:267056
To explore spinal dorsal horn (SDH) network function, we developed a computational model of the circuit that is tightly constrained by experimental data. Our model comprises conductance-based model neurons that reproduce the characteristic firing patterns of excitatory and inhibitory spinal neurons. Excitatory spinal neuron subtypes defined by calretinin, somatostatin, delta-opioid receptor, protein kinase C gamma, or vesicular glutamate transporter 3 expression or by transient/central spiking/morphology and inhibitory neuron subtypes defined by parvalbumin or dynorphin expression or by islet morphology were synaptically connected according to available qualitative data. Synaptic weights were adjusted to produce firing in projection neurons, defined by neurokinin-1 expression, matching experimentally measured responses to a range of mechanical stimulus intensities. Input to the circuit was provided by three types of afferents (Aß, Ad, and C-fibres) whose firing rates were also matched to experimental data.
Reference:
1 . Medlock L, Sekiguchi K, Hong S, Dura-Bernal S, Lytton WW, Prescott SA (2022) Multiscale computer model of the spinal dorsal horn reveals changes in network processing associated with chronic pain The Journal of Neuroscience [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Spinal cord lamina I neuron; Spinal cord lamina I-III interneuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NetPyNE; NEURON;
Model Concept(s): Sensory processing; Pain processing;
Implementer(s): Medlock, Laura [laura.medlock at mail.utoronto.ca]; Sekiguchi, Kazutaka [kazutaka.sekiguchi at shionogi.co.jp];
'''
this is the code for generating the spike times for 10seconds of A and C fibers
created by K. Sekiguchi (18th July 20)
'''

import numpy as np
from cfg_mechanical import cfg

def rate_SAI():
    t = [ x for x in np.arange(0, 10001, 1) ]

    rate = []

    for i, key in enumerate(t):
        if i < len(t):
            freq_SAI = (-1.45433609113392e-10 * key ** 3 + 1.340603396708e-6 * key ** 2 - 0.00378224210238498 * key + 4.52737468545426) * 8
            rate.append(freq_SAI)

    return rate, t

def rate_SAII():
    t = [ x for x in np.arange(0, 10001, 1) ]

    rate = []

    for i, key in enumerate(t):
        if i < len(t):
            freq_SAII = (-1.45433609113392e-10 * key ** 3 + 1.340603396708e-6 * key ** 2 - 0.00378224210238498 * key + 4.52737468545426) * 8
            rate.append(freq_SAII)

    return rate, t

def rate_Ad():
    t = [ x for x in np.arange(0, 10001, 1) ]

    rate = []

    for i, key in enumerate(t):
        if i < len(t):
            freq_Ad = ( -7.265297e-15 * key ** 4 + 1.573831e-10 * key ** 3 - 8.201529e-7 * key ** 2 - 0.002539930 * key + 27.18412 ) * cfg.stim_ratios
            rate.append(freq_Ad)

    return rate, t

def rate_C():
    t = [ x for x in np.arange(0, 10001, 1) ]

    rate  = []

    for i, key in enumerate(t):
        if i < len(t):
            freq_C = (  9.630390e-15 * key ** 4 - 2.844577e-10 * key ** 3 + 3.012887e-6 * key ** 2 - 0.01400381 * key + 31.86546 ) * cfg.stim_ratios
            rate.append(freq_C)

    return rate, t

def poisson_generator(rate, t_start=0.0, t_stop=1000.0, seed=None):
    """
    Returns a SpikeTrain whose spikes are a realization of a Poisson process
    with the given rate (Hz) and stopping time t_stop (milliseconds).
    Note: t_start is always 0.0, thus all realizations are as if 
    they spiked at t=0.0, though this spike is not included in the SpikeList.
    Inputs:
    -------
        rate    - the rate of the discharge (in Hz)
        t_start - the beginning of the SpikeTrain (in ms)
        t_stop  - the end of the SpikeTrain (in ms)
        array   - if True, a np array of sorted spikes is returned,
                    rather than a SpikeTrain object.
    Examples:
    --------
        >> gen.poisson_generator(50, 0, 1000)
        >> gen.poisson_generator(20, 5000, 10000, array=True)
    See also:
    --------
        inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
    """
    rng = np.random.RandomState(seed)
    #number = int((t_stop-t_start)/1000.0*2.0*rate)
    # less wasteful than double length method above
    n = (t_stop-t_start)/1000.0*rate
    number = np.ceil(n+3*np.sqrt(n))
    if number<100:
        number = min(5+np.ceil(2*n),100)
    if number > 0:
        isi = rng.exponential(1.0/rate, int(number))*1000.0
        if number > 1:
            spikes = np.add.accumulate(isi)
        else:
            spikes = isi
    else:
        spikes = np.array([])
    spikes+=t_start
    i = np.searchsorted(spikes, t_stop)
    extra_spikes = []
    if i==len(spikes):
        # ISI buf overrun
        
        t_last = spikes[-1] + rng.exponential(1.0/rate, 1)[0]*1000.0
        while (t_last<t_stop):
            extra_spikes.append(t_last)
            t_last += rng.exponential(1.0/rate, 1)[0]*1000.0
        
        spikes = np.concatenate((spikes,extra_spikes))
    else:
        spikes = np.resize(spikes,(i,))
        
    return spikes

def inh_poisson_generator(rate, t, t_stop, seed=None):
    """
    Returns a SpikeTrain whose spikes are a realization of an inhomogeneous 
    poisson process (dynamic rate). The implementation uses the thinning 
    method, as presented in the references.
    Inputs:
    -------
        rate   - an array of the rates (Hz) where rate[i] is active on interval 
                    [t[i],t[i+1]]
        t      - an array specifying the time bins (in milliseconds) at which to 
                    specify the rate
        t_stop - length of time to simulate process (in ms)
        array  - if True, a np array of sorted spikes is returned,
                    rather than a SpikeList object.
    Note:
    -----
        t_start=t[0]
    References:
    -----------
    Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier 
    Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories
    Neural Comput. 2007 19: 2958-3010.
    Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
    Examples:
    --------
        >> time = arange(0,1000)
        >> stgen.inh_poisson_generator(time,sin(time), 1000)
    See also:
    --------
        poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
    """
    rng = np.random.RandomState(seed)
    if np.shape(t)!=np.shape(rate):
        raise ValueError('shape mismatch: t,rate must be of the same shape')
    # get max rate and generate poisson process to be thinned
    rmax = np.max(rate)
    ps = poisson_generator(rate=rmax, t_start=t[0], t_stop=t_stop, seed=None)
    # return empty if no spikes
    if len(ps) == 0:
        np.array([])
        
    # gen uniform rand on 0,1 for each spike
    rn = np.array(rng.uniform(0, 1, len(ps)))
    # instantaneous rate for each spike
    idx = np.searchsorted(t, ps) - 1
    #spike_rate = rate[idx]
    spike_rate = np.array([rate[i] for i in idx])
    # thin and return spikes
    spike_train = ps[rn<spike_rate/rmax]
    return list(spike_train)