A neural mass model for critical assessment of brain connectivity (Ursino et al 2020)

 Download zip file 
Help downloading and running models
Accession:263637
We use a neural mass model of interconnected regions of interest to simulate reliable neuroelectrical signals in the cortex. In particular, signals simulating mean field potentials were generated assuming two, three or four ROIs, connected via excitatory or by-synaptic inhibitory links. Then we investigated whether bivariate Transfer Entropy (TE) can be used to detect a statistically significant connection from data (as in binary 0/1 networks), and even if connection strength can be quantified (i.e., the occurrence of a linear relationship between TE and connection strength). Results suggest that TE can reliably estimate the strength of connectivity if neural populations work in their linear regions. However, nonlinear phenomena dramatically affect the assessment of connectivity, since they may significantly reduce TE estimation. Software included here allows the simulation of neural mass models with a variable number of ROIs and connections, the estimation of TE using the free package Trentool, and the realization of figures to compare true connectivity with estimated values.
Reference:
1 . Ursino M, Ricci G, Magosso E (2020) Transfer Entropy as a Measure of Brain Connectivity: A Critical Analysis With the Help of Neural Mass Models Front Comput Neurosci . [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neural mass; Connectionist Network; Synapse;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex L5/6 pyramidal GLU cell; Neocortex layer 5 interneuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s): Glutamate; Gaba;
Simulation Environment: MATLAB (web link to model); MATLAB; Trentool;
Model Concept(s): Brain Rhythms; Connectivity matrix; Delay;
Implementer(s): Ursino, Mauro [mauro.ursino at unibo.it]; Ricci, Giulia [Giulia.Ricci at unibo.it]; Magosso, Elisa [elisa.magosso at unibo.it];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Gaba; Glutamate; 
clear all
close all 
clc



Npop=3; % Number of ROIs

% time definition
dt=0.0001;
f_eulero=1/dt;
tend=57; % 56 sec, because one second will be excluded for transitory effects
t=(0:dt:tend);
N=length(t);

rng(11)  % noise seed

%% parameters definition
% Connectivity constants 
C(:,1) = 40.*ones(1,Npop); %Cep
C(:,2) = 40.*ones(1,Npop); %Cpe
C(:,3) = 40.*ones(1,Npop); %Csp
C(:,4) = 50.*ones(1,Npop); %Cps   
C(:,5) = 20.*ones(1,Npop); %Cfs
C(:,6) = 40.*ones(1,Npop); %Cfp
C(:,7) = 60.*ones(1,Npop); %Cpf
C(:,8) = 20.*ones(1,Npop); %Cff
 

% definition of excitatory and inhibitory synapses

Wp=zeros(Npop);

%excitatorty synapses
Wp(1,2)=0;   % Wp_12 = 0, 20, 40, 60, 80 : 5 points (a,b,c,d,e) for each panel
Wp(2,1)=40;
Wp(1,3)=0;  % 0, 20, 40, 60 : 4 panels 
Wp(2,3)=0;  % 0, 20, 40, 60 : 4 panels 

% inhibitory synapses
Wf=zeros(Npop); 


e0 = 2.5; % Saturation value of the sigmoid
r = 0.56; % Slope of the sigmoid

D=0.0166*ones(1,Npop); % Delay between regions
                
a=[75 30 300 ]; % Reciprocal of synaptic time constants (\omega)             

G=[5.17 4.45 57.1]; % Synaptic gains
                   
 %% Simulation
for trial = 1: 10

disp(trial)

sigma = sqrt(9/dt); % Standard deviation of the input noise
np = randn(Npop,N)*sigma; % input noise to excitatory neurons
nf = randn(Npop,N)*sigma; % input noise to inhibitory neurons

% defining equations of a single ROI
yp=zeros(Npop,N);
xp=zeros(Npop,N);
vp=zeros(Npop,1);
zp=zeros(Npop,N);
ye=zeros(Npop,N);
xe=zeros(Npop,N);
ve=zeros(Npop,1);
ze=zeros(Npop,N);
ys=zeros(Npop,N);
xs=zeros(Npop,N);
vs=zeros(Npop,1);
zs=zeros(Npop,N);
yf=zeros(Npop,N);
xf=zeros(Npop,N);
zf=zeros(Npop,N);
vf=zeros(Npop,1);
xl=zeros(Npop,N);
yl=zeros(Npop,N);

step_red = 100;  % step reduction from 10000 to 100 Hz
fs = f_eulero/step_red;
eeg=zeros(Npop,(N-1-10000)/step_red);  % exclusion of the first second due to a possible transitory

m = zeros(Npop,1); % mean value of the input noise

kmax=round(max(D)/dt);

for k=1:N-1
   up=np(:,k)+m; % input of exogenous contributions to excitatory neurons
   uf=nf(:,k);  % input of exogenous contributions to inhibitory neurons
    
    if(k>kmax)
        for i=1:Npop
            up(i)=up(i)+Wp(i,:)*zp(:,round(k-D(i)/dt));
            uf(i)=uf(i)+Wf(i,:)*zp(:,round(k-D(i)/dt));
        end
    end
   
    % post-synaptic membrane potentials
    vp(:)=C(:,2).*ye(:,k)-C(:,4).*ys(:,k)-C(:,7).*yf(:,k);
    ve(:)=C(:,1).*yp(:,k);
    vs(:)=C(:,3).*yp(:,k);
    vf(:)=C(:,6).*yp(:,k)-C(:,5).*ys(:,k)-C(:,8).*yf(:,k)+yl(:,k);
    
    % average spike density
    zp(:,k)=2*e0./(1+exp(-r*(vp(:))))-e0;
    ze(:,k)=2*e0./(1+exp(-r*(ve(:))))-e0;
    zs(:,k)=2*e0./(1+exp(-r*(vs(:))))-e0;
    zf(:,k)=2*e0./(1+exp(-r*(vf(:))))-e0;
    
    % post synaptic potential for pyramidal neurons
    xp(:,k+1)=xp(:,k)+(G(1)*a(1)*zp(:,k)-2*a(1)*xp(:,k)-a(1)*a(1)*yp(:,k))*dt;   
    yp(:,k+1)=yp(:,k)+xp(:,k)*dt; 
    
    % post synaptic potential for excitatory interneurons
    xe(:,k+1)=xe(:,k)+(G(1)*a(1)*(ze(:,k)+up(:)./C(:,2))-2*a(1)*xe(:,k)-a(1)*a(1)*ye(:,k))*dt;  
    ye(:,k+1)=ye(:,k)+xe(:,k)*dt; 
    
    % post synaptic potential for slow inhibitory interneurons
    xs(:,k+1)=xs(:,k)+(G(2)*a(2)*zs(:,k)-2*a(2)*xs(:,k)-a(2)*a(2)*ys(:,k))*dt;   
    ys(:,k+1)=ys(:,k)+xs(:,k)*dt; 
    
    % post synaptic potential for fast inhibitory interneurons
    xl(:,k+1)=xl(:,k)+(G(1)*a(1)*uf(:)-2*a(1)*xl(:,k)-a(1)*a(1)*yl(:,k))*dt;  
    yl(:,k+1)=yl(:,k)+xl(:,k)*dt; 
    xf(:,k+1)=xf(:,k)+(G(3)*a(3)*zf(:,k)-2*a(3)*xf(:,k)-a(3)*a(3)*yf(:,k))*dt;   
    yf(:,k+1)=yf(:,k)+xf(:,k)*dt; 


end

% 3 ROIs data generation
start = 10000; % exclusion of the first second due to a possible transitory
eeg=diag(C(:,2))*ye(:,start:step_red:end)-diag(C(:,4))*ys(:,start:step_red:end)-diag(C(:,7))*yf(:,start:step_red:end);

switch trial
    case 1
        eeg1 = eeg;
    case 2
        eeg2 = eeg;
    case 3
        eeg3 = eeg;
    case 4
        eeg4 = eeg;
    case 5
        eeg5 = eeg;
    case 6
        eeg6 = eeg;
    case 7
        eeg7 = eeg;
    case 8
        eeg8 = eeg;
    case 9
        eeg9 = eeg;
    case 10
        eeg10 = eeg;
end
end

tt=t(start:step_red:end); % time vector

% save data for TRENTOOL
save sim_data_Fig5_0a eeg1 eeg2 eeg3 eeg4 eeg5 eeg6 eeg7 eeg8 eeg9 eeg10 tt Wf Wp