Prosthetic electrostimulation for information flow repair in a neocortical simulation (Kerr 2012)

 Download zip file   Auto-launch 
Help downloading and running models
This model is an extension of a model (<a href="">138379</a>) recently published in Frontiers in Computational Neuroscience. This model consists of 4700 event-driven, rule-based neurons, wired according to anatomical data, and driven by both white-noise synaptic inputs and a sensory signal recorded from a rat thalamus. Its purpose is to explore the effects of cortical damage, along with the repair of this damage via a neuroprosthesis.
1 . Kerr CC, Neymotin SA, Chadderdon GL, Fietkiewicz CT, Francis JT, Lytton WW (2012) Electrostimulation as a prosthesis for repair of information flow in a computer model of neocortex IEEE Transactions on Neural Systems & Rehabilitation Engineering 20(2):153-60 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 pyramidal corticothalamic L6 cell; Neocortex V1 pyramidal intratelencephalic L2-5 cell; Neocortex V1 interneuron basket PV cell; Neocortex fast spiking (FS) interneuron; Neocortex spiny stellate cell;
Channel(s): I Chloride; I Sodium; I Potassium;
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA; Gaba;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Deep brain stimulation; Information transfer; Brain Rhythms;
Implementer(s): Lytton, William [billl at]; Neymotin, Sam [samn at]; Kerr, Cliff [cliffk at];
Search NeuronDB for information about:  Neocortex V1 pyramidal corticothalamic L6 cell; Neocortex V1 pyramidal intratelencephalic L2-5 cell; Neocortex V1 interneuron basket PV cell; GabaA; AMPA; NMDA; Gaba; I Chloride; I Sodium; I Potassium; Gaba; Glutamate;
infot.mod *
intf6_.mod *
intfsw.mod *
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
decmat.hoc *
decnqs.hoc *
default.hoc *
drline.hoc *
infot.hoc *
local.hoc *
misc.h *
ratlfp.dat *
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
xgetargs.hoc *
This file contains all the function definitions necessary for running spectral
Granger causality. It is based on Mingzhou Ding's Matlab code package BSMART,
available from

Typical usage is as follows:
from bsmart import pwcausalr

	F is the frequency vector for the remaining quantities
	pp is the spectral power
	cohe is the coherence
	Fx2y is the causality of channel X to channel Y
	Fy2x is the causality of channel Y to channel X
	Fxy is the "instantaneous" causality (cohe-Fx2y-Fy2x I think)
	x is the data for at least two channels, e.g. a 2x8000 array consisting of two LFP time series
	ntrls is the number of trials (whatever that means -- just leave it at 1)
	npts is the number of points in the data (in this example, 8000)
	p is the order of the polynomial fit (e.g. 10 for a smooth fit, 20 for a less smooth fit)
	fs is the sampling rate (e.g. 200 Hz)
	freq is the maximum frequency to calculate (e.g. fs/2=100, which will return 0:100 Hz)

The other two functions (armorf and spectrum_AR) can also be called directly, but
more typically they are used by pwcausalr in intermediate calculations. Note that the 
sampling rate of the returned quantities is calculated as fs/2.

To calculate the power spectrum powspec of a single time series x over the frequency range 0:freq, 
use the following (NB: now accessible via "from spectrum import ar")
from bsmart import armorf, spectrum_AR
[A,Z,tmp]=armorf(x,ntrls,npts,p) # Calculate autoregressive fit
for i in range(freq+1): # Loop over frequencies
    [S,H]=spectrum_AR(A,Z,p,i,fs) # Calculate spectrum
    powspec[i]=abs(S**2) # Calculate and store power

In either case (pwcausalr or spectrum_AR), the smoothness of the spectra is determined by the
polynomial order p. Larger values of p give less-smooth spectra.

Version: 2012dec10 by Cliff Kerr (

# ARMORF -- AR parameter estimation via LWR method modified by Morf. 
#   X is a matrix whose every row is one variable's time series 
#   ntrls is the number of realizations, npts is the length of every realization 
#   If the time series are stationary long, just let ntrls=1, npts=length(x) 
#   A = ARMORF(X,NR,NL,ORDER) returns the polynomial coefficients A corresponding to 
#   the AR model estimate of matrix X using Morf's method. 
#   ORDER is the order of the AR model. 
#   [A,E] = ARMORF(...) returns the final prediction error E (the variance 
#   estimate of the white noise input to the AR model). 
#   [A,E,K] = ARMORF(...) returns the vector K of reflection coefficients (parcor coefficients). 
#   Ref: M. Morf, etal, Recursive Multichannel Maximum Entropy Spectral Estimation, 
#              IEEE trans. GeoSci. Elec., 1978, Vol.GE-16, No.2, pp85-94. 
#        S. Haykin, Nonlinear Methods of Spectral Analysis, 2nd Ed. 
#              Springer-Verlag, 1983, Chapter 2 

def timefreq(x,fs=200):
    This function takes the time series and the sampling rate and calculates the
    total number of points, the maximum frequency, the minimum (or change in)
    frequency, and the vector of frequency points F.
    Version: 2011may04
    from numpy import size, arange, append
    maxfreq=float(fs)/2.0 # Maximum frequency
    minfreq=float(fs)/float(size(x,0)) # Minimum and delta frequency -- simply the inverse of the length of the recording in seconds
    F=arange(minfreq,maxfreq+minfreq,minfreq) # Create frequencies evenly spaced from 0:minfreq:maxfreq
    F=append(0,F) # Add zero-frequency component
    return F

def ckchol(M):
    This function computes the Cholesky decomposition of the matrix if it's
    positive-definite; else it returns the identity matrix. It was written
    to handle the "matrix must be positive definite" error in linalg.cholesky.
    Version: 2011may03
    from numpy import linalg, matrix, eye, size

    try: # First, try the Cholesky decomposition
    except: # If not, just return garbage
        print 'WARNING: Cholesky failed, so returning (invalid) identity matrix!'    
    return output

def armorf(x,ntrls,npts,p):
    from scipy import shape, array, matrix, zeros, concatenate, eye, dstack
    from numpy import linalg # for inverse and Cholesky factorization;
    inv=linalg.inv; # Make name consistent with Matlab
    # Initialization 
    [L,N]=shape(x);      # L is the number of channels, N is the npts*ntrls 
    pf=pb=pfb=ap=bp=En=matrix(zeros((L,L,1)));    # covariance matrix at 0, 
    # calculate the covariance matrix? 
    for i in range(ntrls):
    ap = inv((ckchol(ap/ntrls*(npts-1)).T).H); 
    bp = inv((ckchol(bp/ntrls*(npts-1)).T).H); 
    for i in range(ntrls): 
       efp = ap*x[:,i*npts+1:(i+1)*npts]; 
       ebp = bp*x[:,i*npts:(i+1)*npts-1]; 
       pf = pf + efp*efp.H; 
       pb = pb + ebp*ebp.H; 
       pfb = pfb + efp*ebp.H; 
    En = (ckchol(En/N).T).H;       # Covariance of the noise 
    # Initial output variables
    for i in range(L): tmp.append([]) # In Matlab, coeff=[], and anything can be appended to that.
    coeff = matrix(tmp);#  Coefficient matrices of the AR model 
    kr = matrix(tmp);  # reflection coefficients 
    aparr=array(ap) # Convert AP matrix to an array, so it can be dstacked
    for m in range(p): 
      # Calculate the next order reflection (parcor) coefficient 
      ck = inv((ckchol(pf).T).H)*pfb*inv(ckchol(pb).T);  
      # Update the forward and backward prediction errors 
      ef = eye(L)- ck*ck.H; 
      eb = eye(L)- ck.H*ck; 
      # Update the prediction error 
      En = En*(ckchol(ef).T).H;

      # Update the coefficients of the forward and backward prediction errors 
      Z=zeros((L,L)) # Make it easier to define this
      pf = pb = pfb = Z
      # Do some variable juggling to handle Python's array/matrix limitations

      for i in range(m+2):  
          tmpap1=matrix(aparr[:,:,i]) # Need to convert back to matrix to perform operations
          tmpa = inv((ckchol(ef).T).H)*(tmpap1-ck*tmpbp2); 
          tmpb = inv((ckchol(eb).T).H)*(tmpbp1-ck.H*tmpap2); 

      for k in range(ntrls):
          efp = zeros((L,npts-m-2)); 
          ebp = zeros((L,npts-m-2)); 
          for i in range(m+2): 
              efp = efp+matrix(a[:,:,i])*matrix(x[:,k1:k2]); 
              ebp = ebp+matrix(b[:,:,m+1-i])*matrix(x[:,k1-1:k2-1]); 
          pf = pf + efp*efp.H; 
          pb = pb + ebp*ebp.H; 
          pfb = pfb + efp*ebp.H; 

      aparr = a; 
      bparr = b; 
    for j in range(p):
       coeff = concatenate((coeff,inv(matrix(a[:,:,0]))*matrix(a[:,:,j+1])),1); 

    return coeff, En*En.H, kr

#Port of spectrum_AR.m
# Version: 2010jan18

def spectrum_AR(A,Z,M,f,fs): # Get the spectrum in one specific frequency-f
    from scipy import eye, size, exp, pi
    from numpy import linalg; inv=linalg.inv
    N = size(Z,0); H = eye(N,N); # identity matrix
    for m in range(M):
        H = H + A[:,m*N:(m+1)*N]*exp(-1j*(m+1)*2*pi*f/fs);   # Multiply f in the exponent by sampling interval (=1/fs). See Richard Shiavi
    H = inv(H);
    S = H*Z*H.H/fs;
    return S,H

# Using Geweke's method to compute the causality between any two channels 
#   x is a two dimentional matrix whose each row is one variable's time series 
#   Nr is the number of realizations, 
#   Nl is the length of every realization 
#      If the time series have one ralization and are stationary long, just let Nr=1, Nl=length(x) 
#   porder is the order of AR model 
#   fs is sampling frequency 
#   freq is a vector of frequencies of interest, usually freq=0:fs/2 
#       CK: WRONG!! freq must be a scalar, else the for loop doesn't work.
#   Fx2y is the causality measure from x to y 
#   Fy2x is causality from y to x 
#   Fxy is instantaneous causality between x and y 
#        the order of Fx2y/Fy2x is 1 to 2:L, 2 to 3:L,....,L-1 to L.  That is, 
#        1st column: 1&2; 2nd: 1&3; ...; (L-1)th: 1&L; ...; (L(L-1))th: (L-1)&L. 

# revised Jan. 2006 by Yonghong Chen 
# Note: remove the ensemble mean before using this code 

def pwcausalr(x,Nr,Nl,porder,fs,freq=0): # Note: freq determines whether the frequency points are calculated or chosen
    from pylab import size, shape, real, log, conj, zeros, arange, array
    from numpy import linalg; det=linalg.det
    import numpy as np # Just for "sum"; can't remember what's wrong with pylab's sum
    [L,N] = shape(x); #L is the number of channels, N is the total points in every channel 
    if freq==0: F=timefreq(x[0,:],fs) # Define the frequency points
    else: F=array(range(0,freq+1)) # Or just pick them
    # Initialize arrays
    # Had these all defined on one line, and stupidly they STAY linked!!
    index = 0;

    for i in range(1,L):
        for j in range(i+1,L+1):
            y=zeros((2,N)) # Initialize y
            index = index + 1; 
            y[0,:] = x[i-1,:]; 
            y[1,:] = x[j-1,:];   
            A2,Z2,tmp = armorf(y,Nr,Nl,porder); #fitting a model on every possible pair 
            eyx = Z2[1,1] - Z2[0,1]**2/Z2[0,0]; #corrected covariance 
            exy = Z2[0,0] - Z2[1,0]**2/Z2[1,1]; 
            f_ind = 0; 
            for f in F:
                f_ind = f_ind + 1; 
                S2,H2 = spectrum_AR(A2,Z2,porder,f,fs); 
                pp[i-1,f_ind-1] = abs(S2[0,0]*2);      # revised 
                if (i==L-1) & (j==L):
                    pp[j-1,f_ind-1] = abs(S2[1,1]*2);  # revised 
                cohe[index-1,f_ind-1] = real(abs(S2[0,1])**2 / S2[0,0]/S2[1,1]);   
                Fy2x[index-1,f_ind-1] = log(abs(S2[0,0])/abs(S2[0,0]-(H2[0,1]*eyx*conj(H2[0,1]))/fs)); #Geweke's original measure 
                Fx2y[index-1,f_ind-1] = log(abs(S2[1,1])/abs(S2[1,1]-(H2[1,0]*exy*conj(H2[1,0]))/fs)); 
                Fxy[index-1,f_ind-1] = log(abs(S2[0,0]-(H2[0,1]*eyx*conj(H2[0,1]))/fs)*abs(S2[1,1]-(H2[1,0]*exy*conj(H2[1,0]))/fs)/abs(det(S2))); 
    return F,pp,cohe,Fx2y,Fy2x,Fxy

def granger(vec1,vec2,order=10,rate=200,maxfreq=0):
    Provide a simple way of calculating the key quantities.
        F is a 1xN vector of frequencies
        pp is a 2xN array of power spectra
        cohe is the coherence between vec1 and vec2
        Fx2y is the causality from vec1->vec2
        Fy2x is the causality from vec2->vec1
        Fxy is non-directional causality (cohe-Fx2y-Fy2x)
        vec1 is a time series of length N
        vec2 is another time series of length N
        rate is the sampling rate, in Hz
        maxfreq is the maximum frequency to be returned, in Hz
    Version: 2011jul18
    from bsmart import timefreq, pwcausalr
    from scipy import array, size
    if maxfreq==0: F=timefreq(vec1,rate) # Define the frequency points
    else: F=array(range(0,maxfreq+1)) # Or just pick them
    return F,pp[0,:],cohe[0,:],Fx2y[0,:],Fy2x[0,:],Fxy[0,:]

Kerr CC, Neymotin SA, Chadderdon GL, Fietkiewicz CT, Francis JT, Lytton WW (2012) Electrostimulation as a prosthesis for repair of information flow in a computer model of neocortex IEEE Transactions on Neural Systems & Rehabilitation Engineering 20(2):153-60[PubMed]

References and models cited by this paper

References and models that cite this paper

Adesnik H, Scanziani M (2010) Lateral competition for cortical space by layer-specific horizontal circuits. Nature 464:1155-60 [PubMed]

Binzegger T, Douglas RJ, Martin KA (2004) A quantitative map of the circuit of cat primary visual cortex. J Neurosci 24:8441-53 [PubMed]

Carnevale NT, Hines ML (2006) The NEURON Book

Cui J, Xu L, Bressler SL, Ding M, Liang H (2008) BSMART: a Matlab-C toolbox for analysis of multichannel neural time series. Neural Netw 21:1094-104 [PubMed]

Francis JT, Xu S, Chapin JK (2008) Proprioceptive and cutaneous representations in the rat ventral posterolateral thalamus. J Neurophysiol 99:2291-304 [PubMed]

Gisiger T, Boukadoum M (2011) Mechanisms Gating the Flow of Information in the Cortex: What They Might Look Like and What Their Uses may be. Front Comput Neurosci 5:1-304 [PubMed]

Hines ML, Carnevale NT (2001) NEURON: a tool for neuroscientists. Neuroscientist 7:123-35 [Journal] [PubMed]

   Spatial gridding and temporal accuracy in NEURON (Hines and Carnevale 2001) [Model]

Kamiński M, Ding M, Truccolo WA, Bressler SL (2001) Evaluating causal relations in neural systems: granger causality, directed transfer function and statistical assessment of significance. Biol Cybern 85:145-57 [PubMed]

Lefort S, Tomm C, Floyd Sarria JC, Petersen CC (2009) The excitatory neuronal network of the C2 barrel column in mouse primary somatosensory cortex. Neuron 61:301-16 [PubMed]

Lizier JT, Heinzle J, Horstmann A, Haynes JD, Prokopenko M (2011) Multivariate information-theoretic measures reveal directed information structure and task relevant changes in fMRI connectivity. J Comput Neurosci 30:85-107

Lloyd KG, Davidson L, Hornykiewicz O (1975) The neurochemistry of Parkinson's disease: effect of L-dopa therapy. J Pharmacol Exp Ther 195:453-64 [PubMed]

Lytton WW, Neymotin SA, Hines ML (2008) The virtual slice setup. J Neurosci Methods 171:309-15 [Journal] [PubMed]

   The virtual slice setup (Lytton et al. 2008) [Model]

Lytton WW, Omurtag A (2007) Tonic-clonic transitions in computer simulation. J Clin Neurophysiol 24:175-81 [PubMed]

   Tonic-clonic transitions in a seizure simulation (Lytton and Omurtag 2007) [Model]

Lytton WW, Omurtag A, Neymotin SA, Hines ML (2008) Just in time connectivity for large spiking networks Neural Comput 20(11):2745-56 [Journal] [PubMed]

   JitCon: Just in time connectivity for large spiking networks (Lytton et al. 2008) [Model]

Lytton WW, Stewart M (2005) A rule-based firing model for neural networks Int J Bioelectromagn 7:47-50

Lytton WW, Stewart M (2006) Rule-based firing for network simulations. Neurocomputing 69:1160-1164

Meyer JS, Obara K, Muramatsu K (1993) Diaschisis. Neurol Res 15:362-6 [PubMed]

Neymotin SA, Jacobs KM, Fenton AA, Lytton WW (2011) Synaptic information transfer in computer models of neocortical columns. J Comput Neurosci. 30(1):69-84 [Journal] [PubMed]

   Synaptic information transfer in computer models of neocortical columns (Neymotin et al. 2010) [Model]

Quilodran R, Gariel MA, Markov NT, Falchier A, Vezoli J, Sallet J, Anderson JC, Dehay C, Doug (2008) Strong loops in the neocortex Society for Neuroscience Abstracts 853.4

Rasche D, Rinaldi PC, Young RF, Tronnier VM (2006) Deep brain stimulation for the treatment of various chronic pain syndromes. Neurosurg Focus 21:E8

Richman JS, Moorman JR (2000) Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol Heart Circ Physiol 278:H2039-49 [PubMed]

Schiff ND, Giacino JT, Kalmar K, Victor JD, Baker K, Gerber M, Fritz B, Eisenberg B, Biondi T (2007) Behavioural improvements with thalamic stimulation after severe traumatic brain injury. Nature 448:600-3 [PubMed]

Schroeder CE, Mehta AD, Foxe JJ (2001) Determinants and mechanisms of attentional modulation of neural processing. Front Biosci 6:D672-84

Shipp S (2005) The importance of being agranular: a comparative account of visual and motor cortex. Philos Trans R Soc Lond B Biol Sci 360:797-814 [PubMed]

Stefani A, Lozano AM, Peppe A, Stanzione P, Galati S, Tropepi D, Pierantozzi M, Brusa L, Scar (2007) Bilateral deep brain stimulation of the pedunculopontine and subthalamic nuclei in severe Parkinson's disease. Brain 130:1596-607 [PubMed]

Stoerig P, Cowey A (1997) Blindsight in man and monkey. Brain 120 ( Pt 3):535-59 [PubMed]

Traub RD, Contreras D, Cunningham MO, Murray H, Lebeau FE, Roopun A, Bibbig A, et al (2005) A single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles and epileptogenic bursts J Neurophysiol 93(4):2194-232 [Journal] [PubMed]

   A single column thalamocortical network model (Traub et al 2005) [Model]
   Collection of simulated data from a thalamocortical network model (Glabska, Chintaluri, Wojcik 2017) [Model]

Van Essen DC, Anderson CH, Felleman DJ (1992) Information processing in the primate visual system: an integrated systems perspective. Science 255:419-23 [PubMed]

Von_monakow C (1914) Die Lokalisation im Grosshirn und der Abbau der Funktion durch kortikale Herde

Chadderdon GL, Mohan A, Suter BA, Neymotin SA, Kerr CC, Francis JT, Shepherd GM, Lytton WW (2014) Motor cortex microcircuit simulation based on brain activity mapping. Neural Comput 26:1239-62 [Journal] [PubMed]

   Motor cortex microcircuit simulation based on brain activity mapping (Chadderdon et al. 2014) [Model]

Dura-Bernal S, Li K, Neymotin SA, Francis JT, Principe JC, Lytton WW (2016) Restoring behavior via inverse neurocontroller in a lesioned cortical spiking model driving a virtual arm. Front. Neurosci. Neuroprosthetics 10:28 [Journal]

   Cortical model with reinforcement learning drives realistic virtual arm (Dura-Bernal et al 2015) [Model]

Dura-Bernal S, Neymotin SA, Kerr CC, Sivagnanam S, Majumdar A, Francis JT, Lytton WW (2017) Evolutionary algorithm optimization of biological learning parameters in a biomimetic neuroprosthesis. IBM Journal of Research and Development (Computational Neuroscience special issue) 61(2/3):6:1-6:14 [Journal]

   Motor system model with reinforcement learning drives virtual arm (Dura-Bernal et al 2017) [Model]

Dura-Bernal S, Zhou X, Neymotin SA, Przekwas A, Francis JT, Lytton WW (2015) Cortical Spiking Network Interfaced with Virtual Musculoskeletal Arm and Robotic Arm. Front Neurorobot 9:13 [Journal] [PubMed]

   Cortical model with reinforcement learning drives realistic virtual arm (Dura-Bernal et al 2015) [Model]

Kerr CC, Van Albada SJ, Neymotin SA, Chadderdon GL, Robinson PA, Lytton WW (2013) Cortical information flow in Parkinson's disease: a composite network-field model. Front Comput Neurosci 7:39:1-14 [Journal] [PubMed]

   Composite spiking network/neural field model of Parkinsons (Kerr et al 2013) [Model]

Neymotin SA, Chadderdon GL, Kerr CC, Francis JT, Lytton WW (2013) Reinforcement learning of 2-joint virtual arm reaching in a computer model of sensorimotor cortex Neural Computation 25(12):3263-93 [Journal] [PubMed]

   Sensorimotor cortex reinforcement learning of 2-joint virtual arm reaching (Neymotin et al. 2013) [Model]

Neymotin SA, Dura-Bernal S, Lakatos P, Sanger TD, Lytton WW (2016) Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex. Front Pharmacol 7:157 [Journal] [PubMed]

   Multitarget pharmacology for Dystonia in M1 (Neymotin et al 2016) [Model]

Neymotin SA, Lee H, Park E, Fenton AA, Lytton WW (2011) Emergence of physiological oscillation frequencies in a computer model of neocortex. Front Comput Neurosci 5:19-75 [Journal] [PubMed]

   Emergence of physiological oscillation frequencies in neocortex simulations (Neymotin et al. 2011) [Model]

Rowan MS, Neymotin SA, Lytton WW (2014) Electrostimulation to reduce synaptic scaling driven progression of Alzheimer's disease. Front Comput Neurosci 8:39 [Journal] [PubMed]

   Electrostimulation to reduce synaptic scaling driven progression of Alzheimers (Rowan et al. 2014) [Model]

(38 refs)