Striatal GABAergic microcircuit, dopamine-modulated cell assemblies (Humphries et al. 2009)

 Download zip file 
Help downloading and running models
To begin identifying potential dynamically-defined computational elements within the striatum, we constructed a new three-dimensional model of the striatal microcircuit's connectivity, and instantiated this with our dopamine-modulated neuron models of the MSNs and FSIs. A new model of gap junctions between the FSIs was introduced and tuned to experimental data. We introduced a novel multiple spike-train analysis method, and apply this to the outputs of the model to find groups of synchronised neurons at multiple time-scales. We found that, with realistic in vivo background input, small assemblies of synchronised MSNs spontaneously appeared, consistent with experimental observations, and that the number of assemblies and the time-scale of synchronisation was strongly dependent on the simulated concentration of dopamine. We also showed that feed-forward inhibition from the FSIs counter-intuitively increases the firing rate of the MSNs.
1 . Humphries MD, Wood R, Gurney K (2009) Dopamine-modulated dynamic cell assemblies generated by the GABAergic striatal microcircuit. Neural Netw 22:1174-88 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neostriatum fast spiking interneuron;
Gap Junctions: Gap junctions;
Receptor(s): D1; D2; GabaA; AMPA; NMDA; Dopaminergic Receptor;
Transmitter(s): Dopamine; Gaba; Glutamate;
Simulation Environment: MATLAB;
Model Concept(s): Activity Patterns; Temporal Pattern Generation; Synchronization; Spatio-temporal Activity Patterns; Parkinson's; Action Selection/Decision Making; Connectivity matrix;
Implementer(s): Humphries, Mark D [m.d.humphries at]; Wood, Ric [ric.wood at];
Search NeuronDB for information about:  D1; D2; GabaA; AMPA; NMDA; Dopaminergic Receptor; Dopamine; Gaba; Glutamate;
function [grps,maxEV] = multileadevsplit(A,varargin)

% MULTILEADEVSPLIT partition graph using leading eigenvector of modularity matrix
%   [S,U] = MULTILEADEVSPLIT(A) splits the vertices of the graph in adjacency matrix
%   A into n groups, defined in the vector S (each element taking value
%   {1,...,n} to indicate group membership of the vertex at that index). It
%   also returns U, the leading eigenvector of the modularity matrix.
%   MULTILEADEVSPLIT(...,FLAG) where FLAG is a string featuring any
%   combination of:
%       'r': adds the refinement step after each sub-graph,
%               so that each split maximises Q.
%       'a': keeps all splits found, even if Q < 0
%   Notes: 
%   (1) this repeats the leading eigenvector split on each sub-graph
%   until the change in modularity is not positive...
%   (2) or could use singular value decomposition to find eigenvalues,
%   as standard MatLab routine EIG() does not converge for some sample
%   directed graphs (e.g. mac95) even though they are made symmetric.
%   (Note though that using SVD does not split the dolphin graph into 3)
%   References: 
%   (1) Newman, M. E. J. (2006a) "Finding community structure in
%       networks using the eigenvectors of matrices". Phys Rev E, in press.
%   (2) Newman, M. E. J. (2006b) "Modularity and community structure in
%       networks". PNAS, in press.
%   (3) Reichardt & Bornholdt (2006) "Statistical mechanics of community detection" 
%       Phys Rev E. 74, 016110
%   Mark Humphries 27/1/2009

blnRefine = 0;
blnAll = 0;
if nargin >= 2 
    if strfind(varargin{1},'r') blnRefine = 1; end
    if strfind(varargin{1},'a') blnAll = 1; end

% apply correction to A for directed graphs. If graph undirected then this
% has no effect
Abar = (A + A') / 2;    

[n c] = size(A);
grps = zeros(n,1);
ks = sum(Abar);        
m = sum(ks)/2;         

[S,maxEV,B] = splitsubgraph(A,ks,m);

Q = S' * B * S;

% check modularity has increased
if ~blnAll & Q <= 1e-5    % very close or less than zero (deals with rounding errors)
    % if not increased then return nothing
    S = zeros(n,1);

% make subgraphs
n1 = sum(S>0);
n2 = sum(S<=0);
subG1 = A(S>0,S>0);
subG2 = A(S<=0,S<=0);
grps = real(S>0) + 1;     % group 1 and group 2

subqueue = {subG1, subG2};
grpvertices = {find(S>0), find(S<=0)};

while ~isempty(subqueue)
    newA = subqueue{1};
    [S,maxEV,B] = splitsubgraph(newA,ks(grpvertices{1}),m); 
    if blnRefine
        S = refinesplit(S,B);
    % compute Q
    % keyboard
        Q = S' * B * S; 
    % split this group if:
    % (a) all splits are kept, and there is a split OR
    % (b) not all splits are kept, but Q > 0 (Newman 2006a suggestion)  
    if (blnAll & any(S>0) & any(S<0)) | Q > 0
        subG1 = newA(S>0,S>0);
        subG2 = newA(S<=0,S<=0);
        grp1 = grpvertices{1}(S>0);     % the ids of the original vertices
        grp2 = grpvertices{1}(S<=0);    % the ids of the original vertices
        grps(grp1) = max(grps)+1;       % re-number first group, second group retains numbers
        %grps(grp2) = max(grps)+1;
        % delete split graph from stack
        subqueue(1) = [];
        grpvertices(1) = [];
        % add new ones
        grpvertices = [{grp1,grp2} grpvertices];        
        subqueue = [{subG1,subG2} subqueue];
        % don't split any further, delete top of stack
        subqueue(1) = [];
        grpvertices(1) = [];

    % keyboard

function [S,maxEV,B] = splitsubgraph(A,ks,m) 
    % Note: if graph is directed, the correction to A should be
    % applied
    Abar = (A + A') / 2;   
    S = zeros(length(A),1);
    % generate expected graph
    P = expectedA(A);
    % correct P
    Pbar = (P + P') / 2;
    kg = sum(Abar);        % in-degrees of sub-graph
    dg = sum(kg);       % total edges of sub-graph
    % this correction should apply without problem to a directed graph
    % because by this stage the graph will be undirected. The key is pass
    % the correct value of ks to this function.
    correction = kg - ks.*dg/(2*m);     % definitely 2*m here: is divided by 2 at start :-)
    % create modularity matrix - correction applied only to self-self
    % connections
    B = Abar - Pbar - diag(correction);   
    % find eigenvalues (D) and eigenvectors (V) of B
    [V,D] = eig(B);         % assumes symmetry
    eg = diag(D);
    if ~isreal(eg)
        warning('Complex eigenvalues have occurred')
    % keyboard
     % largest positive eigenvalue is last element
    maxeig = max(eg); %% note - can have the same eigenvalue if nodes have no links in subgraph!
    ix = find(eg == maxeig(1));
    maxEV = V(:,ix(1)); % if more than one with the same max eigenvalue, just pick first...
    %%%% use singular value decomposition %%%%%%%%%%%%
%     [u,s,v] = svd(B);
%     sv = diag(s);
%     maxEV = v(:,find(sv == max(sv)));
    S(maxEV > 0) = 1;       % group 1
    S(maxEV <= 0) = -1;     % group 2