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 ( ) 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 Trans Neural Syst Rehabil Eng 20:153-60 [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: Neocortex;
Cell Type(s): Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA 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 [bill.lytton at]; Neymotin, Sam [Samuel.Neymotin at]; Kerr, Cliff [cliffk at];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA 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,:]