Get_balanced_weights.m

[dieloz]

Dear Katharina,

I read your Neuroimage paper and I try to run 'sim_script.m' to play a bit with the model. However I can't find the 'Get_balanced_weights.m'. Could you please confirm me that this function corresponds the Deco et al 2014 JNeurosci (https://www.jneurosci.org/content/34/23/7886) section described below?

"Hence, in the large-scale model of interconnected brain areas, we aim to constraint in each brain area (i) the local feedback inhibitory weight Ji such that IiE − bE/aE = − 0.026 is fulfilled (with a tolerance of ±0.005, implying an excitatory firing rate between 2.63–3.55 Hz). To achieve this, we apply following procedure: we simulate during a period of 10 s the system of stochastic differential Equations 5–10 and compute the averaged level of the input to the local excitatory pool of each brain area, i.e., Ii(E). If IiE − bE/aE = − 0.026 then we upregulate the corresponding local feedback inhibition Ji = Ji + Δ; otherwise, we downregulate Ji = Ji − Δ. We recursively repeat this procedure until the constraint on the input to the local excitatory pool is fulfilled in all 66 brain areas."

https://github.com/Katharinski/tensor-pipeline/blob/63bbdba694a39c54a0befffff35d93a6dc1ac6f5/sim_script.m#L24

Thanks beforehand,

Diego

[Katharinski]

Hi Diego,

Thanks for getting in touch! I'm attaching the function - not sure why it's not online, I'll look into that. I think it's a newer version of the model from the paper you mentioned. I'm attaching the code that I used. Let me know if it works!

Best, Katharina

On Wed, Feb 12, 2020 at 5:11 PM Diego Lozano-Soldevilla < notifications@github.com> wrote:

Dear Katharina,

I read your Neuroimage paper and I try to run 'sim_script.m' to play a bit with the model. However I can't find the 'Get_balanced_weights.m'. Could you please confirm me that this function corresponds the Deco et al 2014 JNeurosci (https://www.jneurosci.org/content/34/23/7886) section described below?

"Hence, in the large-scale model of interconnected brain areas, we aim to constraint in each brain area (i) the local feedback inhibitory weight Ji such that IiE − bE/aE = − 0.026 is fulfilled (with a tolerance of ±0.005, implying an excitatory firing rate between 2.63–3.55 Hz). To achieve this, we apply following procedure: we simulate during a period of 10 s the system of stochastic differential Equations 5–10 and compute the averaged level of the input to the local excitatory pool of each brain area, i.e., Ii(E). If IiE − bE/aE = − 0.026 then we upregulate the corresponding local feedback inhibition Ji = Ji + Δ; otherwise, we downregulate Ji = Ji − Δ. We recursively repeat this procedure until the constraint on the input to the local excitatory pool is fulfilled in all 66 brain areas."

https://github.com/Katharinski/tensor-pipeline/blob/63bbdba694a39c54a0befffff35d93a6dc1ac6f5/sim_script.m#L24

Thanks beforehand,

Diego

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/Katharinski/tensor-pipeline/issues/2?email_source=notifications&email_token=AE5EISYRN5B7VJ36FBHSBBDRCQNUZA5CNFSM4KT4X3MKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4IM7SDGA, or unsubscribe https://github.com/notifications/unsubscribe-auth/AE5EIS2GX5BTEC7Q65FKRWTRCQNUZANCNFSM4KT4X3MA .

[dieloz]

Hi Katharina, Thank you very much for your prompt reply. Where exactly did you attached the code? I can't find it neither here nor in your repository https://github.com/Katharinski/tensor-pipeline. Sorry beforehand in case I miss it Best Diego

[Katharinski]

Ah, sorry, it doesn't allow me to attach "that file type". Just copying it below. Copyright by Adrian Ponce-Alvarez.

Hope it helps!

Best, Katharina

%
%
% Calculates the steady states of the balanced model vs. global couplings
% Saves the rates, the feedback inhibition (J).
%
% Adrian Ponce
%--------------------------------------------------------------------------

function [JI,Se_init,Si_init] = Get_balanced_weights(SC,G)
% since code was written for Hagmann SC, matrix has to be rescaled
sumHagSC = 15.3014;
% CHag = load('Hagmann_matrix');
% sumHagSC = sum(CHag(:));
fac = sumHagSC/sum(SC(:));
SC = SC*fac;

N = size(SC,1);

% Fixed parameters:
%----------------------------

dt=0.1;
tmax=10000;
tspan=0:dt:tmax;

taon=100;
taog=10;
gamma=0.641;
sigma=0.01;
JN=0.15;
J=ones(N,1);
I0=0.382; %%397;
Jexte=1.;
Jexti=0.7;
w=1.4;

% transfer functions:
% transfert function: excitatory
%--------------------------------------------------------------------------
ae=310;
be=125;
de=0.16;
He=@(x) (ae*x-be)./(1-exp(-de*(ae*x-be)));

% transfert function: inhibitory
%--------------------------------------------------------------------------
ai=615;
bi=177;
di=0.087;
Hi=@(x) (ai*x-bi)./(1-exp(-di*(ai*x-bi)));

% all tested global couplings:
numG = length(G);

% initialization:
%-------------------------

curr=zeros(tmax,N);
delta=0.02*ones(N,1);

Se_init = zeros(N,numG);
Si_init = zeros(N,numG);

JI=zeros(N,numG);

kk=0;
for g=G

    display(g)
    kk=kk+1;

    %%% Balance (greedy algorithm)
    % note that we used stochastic equations to estimate the JIs
    % Doing that gives more stable solutions as the JIs for each node will be
    % a function of the variance.

    for k=1:5000
        sn=0.001*ones(N,1);
        sg=0.001*ones(N,1);
        nn=1;
        j=0;
        for i=2:1:length(tspan)
            xn=I0*Jexte+w*JN*sn+g*JN*SC*sn-J.*sg;
            xg=I0*Jexti+JN*sn-sg;
            rn=feval(He,xn);
            rg=feval(Hi,xg);
            sn=sn+dt*(-sn/taon+(1-sn)*gamma.*rn./1000.)+sqrt(dt)*sigma*randn(N,1);
            sn(sn>1) = 1;
            sn(sn<0) = 0;
            sg=sg+dt*(-sg/taog+rg./1000.)+sqrt(dt)*sigma*randn(N,1);
            sg(sg>1) = 1;
            sg(sg<0) = 0;
            j=j+1;
            if j==10
                curr(nn,:)=xn'-125/310;
                nn=nn+1;
                j=0;
            end
        end

        currm=mean(curr,1);
        flag=0;
        for n=1:1:N
            if abs(currm(n)+0.026)>0.005
                if currm(n)<-0.026
                    J(n)=J(n)-delta(n);
                    delta(n)=delta(n)-0.001;
                    if delta(n)<0.001
                        delta(n)=0.001;
                    end
                else
                    J(n)=J(n)+delta(n);
                end
            else
                flag=flag+1;
            end
        end
        display(flag)
        if flag==N
            break;
        end
    end

    Se_init(:,kk)=sn; % store steady end states
    Si_init(:,kk)=sg;
    JI(:,kk)=J; % store feedback inhibition values
end
end

[dieloz]

Thank you very much Katharina!

I just finished the "sim_script.m" run and all looks good.

I want to let you know few minor things:

1.- I was able to run Get_balanced_weights.m using Hagmann's matrix taken from https://doi.org/10.1371/journal.pcbi.1004445.s009.

2.- Your repository still misses some functions (BOLD.m, phie.m and phii.m; DMF_sim.m uses them). Luckily I found them from here: https://github.com/anuragrpatil/resting_state_brain_dynamics_asd/blob/master/DMF_original_code/phii.m https://github.com/anuragrpatil/resting_state_brain_dynamics_asd/blob/master/DMF_original_code/phie.m https://github.com/anuragrpatil/resting_state_brain_dynamics_asd/blob/master/DMF_original_code/BOLD.m Best Diego