function [link_n,link_wAMPA,link_wNMDA,link_wGABAa,link_t,n_link,delays,bounds,n_paths,max_prox,max_soma,M_hat] = ...
BG_net_shunt_AMPA_NMDA2(scalar_constants,ae_AMPA,ae_NMDA,ai_GABAa,as,fast_weights,slow_weights,delay_times,proportions, ...
tau_AMPA,tau_NMDA,tau_GABAa,tau_m,R,stnda,gpeda)
% BG_NET_SHUNT_AMPA_NMDA2 - creates structure of model, with both AMPA and
% NMDA synapes, and allows separate weighting of them...
% constants
p_connect = scalar_constants(1);
dt = scalar_constants(2);
dop1 = scalar_constants(3);
dop2 = scalar_constants(4);
scale = scalar_constants(5);
AMPA_PSP_size = scalar_constants(6);
NMDA_PSP_size = scalar_constants(7);
GABAa_PSP_size = scalar_constants(8);
n_channels = scalar_constants(9);
neurons_per_channel = scalar_constants(10);
neurons_per_nucleus = scalar_constants(11);
n_neurons = scalar_constants(12);
n_sources = scalar_constants(13);
max_scale = scalar_constants(14);
%specify nuclei indices
SD1 = 1:neurons_per_nucleus; %Striatal D1 neurons
SD2 = neurons_per_nucleus+1:2*neurons_per_nucleus; %Striatal D2 neurons
STN = 2*neurons_per_nucleus+1:3*neurons_per_nucleus; %Subthalamic neurons
GPe = 3*neurons_per_nucleus+1:4*neurons_per_nucleus; %Globus Pallidus internus
GPi = 4*neurons_per_nucleus+1:5*neurons_per_nucleus; %Globus Pallidus externus
EXT = 5*neurons_per_nucleus+1:6*neurons_per_nucleus; %Extrinsic input
%Memory efficient cell array connectivity structure
link_n = cell(n_sources,1); %link (target) neurons
link_wAMPA = cell(n_sources,1); %link weights
link_wNMDA = cell(n_sources,1); %link weights
link_wGABAa = cell(n_sources,1); %link weights
link_t = cell(n_sources,1); %link types
n_link = repmat(uint32(0),n_sources,1); %number of links
n_link(SD1) = neurons_per_channel;
n_link(SD2) = neurons_per_channel;
n_link(STN) = 2*neurons_per_nucleus;
n_link(GPe) = 2*neurons_per_channel + neurons_per_nucleus;
n_link(GPi) = neurons_per_nucleus;
n_link(EXT) = 3*neurons_per_channel;
%%%%%%%%% SET THESE TO ZERO TO LESION STRUCTURE %%%%%%%%%%%%%%%%%%%
n_paths = repmat(uint32(0),n_sources,1);
n_paths(SD1) = 1; %GPi
n_paths(SD2) = 1; %GPe
n_paths(STN) = 2; %GPe and GPi
n_paths(GPe) = 3; %STN, GPi, and itself
n_paths(GPi) = 1;
n_paths(EXT) = 3; %SD1, SD2 and STN
%introduce delays and bounds as more efficient alternative to paths
delays = repmat(uint32(0),n_sources,double(max(n_paths)));
bounds = repmat(uint32(0),n_sources,double(max(n_paths))+1); %Slightly bigger for n_link add-on
delays(SD1,1) = delay_times(1) / dt;
delays(SD2,1) = delay_times(2) / dt;
delays(STN,1) = delay_times(3) / dt;
delays(STN,2) = delay_times(4) / dt;
delays(GPe,1) = delay_times(5) / dt;
delays(GPe,2) = delay_times(6) / dt;
delays(GPe,3) = delay_times(7) / dt;
delays(GPi,1) = delay_times(8) / dt;
delays(EXT,1) = delay_times(9) / dt;
delays(EXT,2) = delay_times(10) / dt;
delays(EXT,3) = delay_times(11) / dt;
bounds(:,1) = 0; %all start from zero in mex - all index first entry of a given pathway
bounds(SD1,2) = neurons_per_channel;
bounds(SD2,2) = neurons_per_channel;
bounds(STN,2) = neurons_per_nucleus;
bounds(STN,3) = 2*neurons_per_nucleus;
bounds(GPe,2) = neurons_per_channel;
bounds(GPe,3) = 2*neurons_per_channel;
bounds(GPe,4) = 2*neurons_per_channel+neurons_per_nucleus;
bounds(GPi,2) = neurons_per_nucleus;
bounds(EXT,2) = neurons_per_channel;
bounds(EXT,3) = 2*neurons_per_channel;
bounds(EXT,4) = 3*neurons_per_channel;
%%%% determine weights %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate weights -> [relative strength * PSP-generating current (given mean R, Tm, Ts) * scalar (for no. of synapses)]
SD1_w_GABAa = fast_weights(1) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_GABAa,tau_m(5),R(5),GABAa_PSP_size,'step') *scale;
SD2_w_GABAa = fast_weights(2) * ones(1,neurons_per_nucleus) * PSPtoPSC(tau_GABAa,tau_m(4),R(4),GABAa_PSP_size,'step') *scale;
STN_GPe_w_AMPA = fast_weights(3) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_AMPA,tau_m(4),R(4),AMPA_PSP_size,'step') * scale;
STN_GPe_w_NMDA = slow_weights(1) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_NMDA,tau_m(4),R(4),NMDA_PSP_size,'step') * scale;
STN_GPi_w_AMPA = fast_weights(4) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_AMPA,tau_m(5),R(5),AMPA_PSP_size,'step') * scale;
STN_GPi_w_NMDA = slow_weights(2) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_NMDA,tau_m(5),R(5),NMDA_PSP_size,'step') * scale;
% discrete STN output
%STN_GPew = weights(3) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_e1,tau_m(4),R(4),PSP_size,'step')) * scale;
%STN_GPiw = weights(4) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_e1,tau_m(5),R(5),PSP_size,'step')) * scale;
GPe_STN_w_GABAa = fast_weights(5) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_GABAa,tau_m(3),R(3),GABAa_PSP_size,'step') *scale;
GPe_GPi_w_GABAa = fast_weights(6) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_GABAa,tau_m(5),R(5),GABAa_PSP_size,'step')*scale;
GPe_GPe_w_GABAa = fast_weights(7) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_GABAa,tau_m(4),R(4),GABAa_PSP_size,'step')*scale;
GPi_GPi_w_GABAa = fast_weights(8) * ones(1,neurons_per_nucleus)* PSPtoPSC(tau_GABAa,tau_m(5),R(5),GABAa_PSP_size,'step')*scale;
EXT_w_AMPA = fast_weights(9) * ones(1,3*neurons_per_nucleus)*scale;
EXT_w_AMPA(1:neurons_per_nucleus) = EXT_w_AMPA(1:neurons_per_nucleus)*PSPtoPSC(tau_AMPA,tau_m(1),R(1),AMPA_PSP_size,'step');
EXT_w_AMPA(1+neurons_per_nucleus:2*neurons_per_nucleus) = EXT_w_AMPA(1+neurons_per_nucleus:2*neurons_per_nucleus)*PSPtoPSC(tau_AMPA,tau_m(2),R(2),AMPA_PSP_size,'step');
EXT_w_AMPA(1+2*neurons_per_nucleus:3*neurons_per_nucleus) = EXT_w_AMPA(1+2*neurons_per_nucleus:3*neurons_per_nucleus)*PSPtoPSC(tau_AMPA,tau_m(3),R(3),AMPA_PSP_size,'step');
EXT_w_NMDA = slow_weights(3) * ones(1,3*neurons_per_nucleus)*scale;
EXT_w_NMDA(1:neurons_per_nucleus) = EXT_w_NMDA(1:neurons_per_nucleus)*PSPtoPSC(tau_NMDA,tau_m(1),R(1),NMDA_PSP_size,'step');
EXT_w_NMDA(1+neurons_per_nucleus:2*neurons_per_nucleus) = EXT_w_NMDA(1+neurons_per_nucleus:2*neurons_per_nucleus)*PSPtoPSC(tau_NMDA,tau_m(2),R(2),NMDA_PSP_size,'step');
EXT_w_NMDA(1+2*neurons_per_nucleus:3*neurons_per_nucleus) = EXT_w_NMDA(1+2*neurons_per_nucleus:3*neurons_per_nucleus)*PSPtoPSC(tau_NMDA,tau_m(3),R(3),NMDA_PSP_size,'step');
%dopamine modulation via weight adjustment
% 1. Striatum
EXT_w_AMPA(1:neurons_per_nucleus) = EXT_w_AMPA(1:neurons_per_nucleus)*(1+dop1); %D1
EXT_w_AMPA(1+neurons_per_nucleus:2*neurons_per_nucleus) = EXT_w_AMPA(1+neurons_per_nucleus:2*neurons_per_nucleus)*(1-dop2); %D2
EXT_w_NMDA(1:neurons_per_nucleus) = EXT_w_NMDA(1:neurons_per_nucleus)*(1+dop1); %D1
EXT_w_NMDA(1+neurons_per_nucleus:2*neurons_per_nucleus) = EXT_w_NMDA(1+neurons_per_nucleus:2*neurons_per_nucleus)*(1-dop2); %D2
% 2. STN (includes cortico-subthalamic weight ratio with cortico-striatal
% strength) %%% note that fast_weights(10) is the STN weight scaling
% ratio...
EXT_w_AMPA(1+2*neurons_per_nucleus:3*neurons_per_nucleus) = fast_weights(10) * EXT_w_AMPA(1+2*neurons_per_nucleus:3*neurons_per_nucleus)*(1-stnda(1)*dop2);
EXT_w_NMDA(1+2*neurons_per_nucleus:3*neurons_per_nucleus) = slow_weights(4) * EXT_w_NMDA(1+2*neurons_per_nucleus:3*neurons_per_nucleus)*(1-stnda(1)*dop2);
GPe_STN_w_GABAa = GPe_STN_w_GABAa .* (1-stnda(2)*dop2);
% 3. GPe
SD2_w_GABAa = SD2_w_GABAa .* (1 - gpeda(1)*dop2);
STN_GPe_w_AMPA = STN_GPe_w_AMPA .* (1 - gpeda(2)*dop2);
STN_GPe_w_NMDA = STN_GPe_w_NMDA .* (1 - gpeda(2)*dop2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Connections:
%1. connect, on average, to proportion p_connect neurons in the corresponding channel
max_prox = zeros(n_neurons,1);
max_soma = zeros(n_neurons,1);
for c = 1:n_channels
% input indexes
in = (c-1)*neurons_per_channel+(1:neurons_per_channel);
for n = 1:neurons_per_channel
out = (c-1)*neurons_per_channel + n; % current out neuron
% channel-wise input
ext_probs = rand(3,neurons_per_channel);
blnExt = ext_probs <= p_connect;
EXT_STN_in = uint32((STN(in)-1) .* blnExt(1,:));
EXT_SD1_in = uint32((SD1(in)-1) .* blnExt(2,:));
EXT_SD2_in = uint32((SD2(in)-1) .* blnExt(3,:));
EXT_out = [EXT_SD1_in,EXT_SD2_in,EXT_STN_in];
% channel-wise output
probs = rand(neurons_per_channel,4); % 4 channel-wise connections to sparsify
blnIn = probs <= p_connect; % logical array of selected unit
D1_GPi_in = uint32((GPi(in)-1).*blnIn(:,1)'); %shifted back to zero-base for mex file
D2_GPe_in = uint32((GPe(in)-1).*blnIn(:,2)');
GPe_STN_in = uint32((STN(in)-1).*blnIn(:,3)');
GPe_GPi_in = uint32((GPi(in)-1).*blnIn(:,4)');
GPe_out = [GPe_STN_in GPe_GPi_in];
% diffuse connections
diffuse_probs = rand(neurons_per_nucleus,4); % 3 diffuse connections to sparsify
blnDiffuseIn = diffuse_probs <= (p_connect / n_channels); % scale probability by number of channels
% diffuse STN output
STN_out = [uint32((GPe-1).*blnDiffuseIn(:,1)') uint32((GPi-1).*blnDiffuseIn(:,2)')];
% diffuse GPe-GPe
blnDiffuseIn(out,3) = 0; % do not connect to self
GPe_out = [GPe_out uint32((GPe-1).*blnDiffuseIn(:,3)')];
% diffuse GPi-GPi
blnDiffuseIn(out,4) = 0; % do not connect to self
GPi_out = uint32((GPi-1).*blnDiffuseIn(:,4)');
% discrete STN output
% stn_probs = rand(neurons_per_channel,2); % 2 discrete connections
% blnSTN_in = stn_probs <= p_connect;
% STN_out = [uint32((GPe(in)-1).*blnSTN_in(:,1)') uint32((GPi(in)-1).*blnSTN_in(:,2)')];
% assign connections
syn = rand(neurons_per_channel,4) .* blnIn; % all of blnIn as it defines all inhibitory connections
syn_diffuse = rand(neurons_per_nucleus,2) .* blnDiffuseIn(:,3:4);
link_n{SD1(out)} = D1_GPi_in;
link_wGABAa{SD1(out)} = SD1_w_GABAa(in) .* ai_GABAa(GPi(in))' .* blnIn(:,1)';
link_wAMPA{SD1(out)} = 0;
link_wNMDA{SD1(out)} = 0;
type1 = zeros(neurons_per_channel,1);
type1(syn(:,1) > 0 & syn(:,1) <= proportions(1,1)) = 1; % distal
type1(syn(:,1) > proportions(1,1) & syn(:,1) <= (proportions(1,1) + proportions(1,2))) = 2; % proximal
type1(syn(:,1) > (proportions(1,1) + proportions(1,2))) = 3; % somatic
link_t{SD1(out)} = uint32(type1);
link_n{SD2(out)} = D2_GPe_in;
link_wGABAa{SD2(out)} = SD2_w_GABAa(in) .*ai_GABAa(GPe(in))' .* blnIn(:,2)';
link_wAMPA{SD2(out)} = 0;
link_wNMDA{SD2(out)} = 0;
type2 = zeros(neurons_per_channel,1);
type2(syn(:,2) > 0 & syn(:,2) <= proportions(2,1)) = 1; % distal
type2(syn(:,2) > proportions(2,1) & syn(:,2) <= (proportions(2,1) + proportions(2,2))) = 2; % proximal
type2(syn(:,2) > (proportions(2,1) + proportions(2,2))) = 3; % somatic
link_t{SD2(out)} = uint32(type2);
link_n{GPe(out)} = GPe_out;
link_wGABAa{GPe(out)} = [GPe_STN_w_GABAa(in) .* ai_GABAa(STN(in))' .* blnIn(:,3)',...
GPe_GPi_w_GABAa(in) .* ai_GABAa(GPi(in))' .* blnIn(:,4)',...
GPe_GPe_w_GABAa .* ai_GABAa(GPe)' .* blnDiffuseIn(:,3)'];
link_wAMPA{GPe(out)} = 0;
link_wNMDA{GPe(out)} = 0;
type3 = zeros(neurons_per_channel,1); % TO STN
type3(syn(:,3) > 0 & syn(:,3) <= proportions(3,1)) = 1; % distal
type3(syn(:,3) > proportions(3,1) & syn(:,3) <= (proportions(3,1) + proportions(3,2))) = 2; % proximal
type3(syn(:,3) > (proportions(3,1) + proportions(3,2))) = 3; % somatic
type4 = zeros(neurons_per_channel,1); % TO GPI
type4(syn(:,4) > 0 & syn(:,4) <= proportions(4,1)) = 1; % distal
type4(syn(:,4) > proportions(4,1) & syn(:,4) <= (proportions(4,1) + proportions(4,2))) = 2; % proximal
type4(syn(:,4) > (proportions(4,1) + proportions(4,2))) = 3; % somatic
type5 = zeros(neurons_per_nucleus,1); % TO ITSELF
type5(syn_diffuse(:,1) > 0 & syn_diffuse(:,1) <= proportions(5,1)) = 1; % distal
type5(syn_diffuse(:,1) > proportions(5,1) & syn_diffuse(:,1) <= (proportions(5,1) + proportions(5,2))) = 2; % proximal
type5(syn_diffuse(:,1) > (proportions(5,1) + proportions(5,2))) = 3; % somatic
link_t{GPe(out)} = [uint32(type3)' uint32(type4)' uint32(type5)'];
link_n{GPi(out)} = GPi_out;
link_wGABAa{GPi(out)} = GPi_GPi_w_GABAa .* ai_GABAa(GPi)' .* blnDiffuseIn(:,4)';
link_wAMPA{GPi(out)} = 0;
link_wNMDA{GPi(out)} = 0;
type6 = zeros(neurons_per_nucleus,1);
type6(syn_diffuse(:,2) > 0 & syn_diffuse(:,2) <= proportions(6,1)) = 1; % distal
type6(syn_diffuse(:,2) > proportions(6,1) & syn_diffuse(:,2) <= (proportions(6,1) + proportions(6,2))) = 2; % proximal
type6(syn_diffuse(:,2) > (proportions(6,1) + proportions(6,2))) = 3; % somatic
link_t{GPi(out)} = uint32(type6)';
% keyboard
% diffuse STN
link_n{STN(out)} = STN_out;
link_wAMPA{STN(out)} = [STN_GPe_w_AMPA .*ae_AMPA(GPe)'.*blnDiffuseIn(:,1)',...
STN_GPi_w_AMPA .*ae_AMPA(GPi)'.*blnDiffuseIn(:,2)'];
link_wNMDA{STN(out)} = [STN_GPe_w_NMDA .*ae_NMDA(GPe)'.*blnDiffuseIn(:,1)',...
STN_GPi_w_NMDA .*ae_NMDA(GPi)'.*blnDiffuseIn(:,2)'];
link_wGABAa{STN(out)} = 0;
link_t{STN(out)} = zeros(neurons_per_nucleus*2,1); % STN output excitatory
% keyboard
% discrete STN
% link_n{STN(out)} = STN_out;
% link_w{STN(out)} = [STN_GPew.*ae(GPe(in))'.*blnSTN_in(:,1)', STN_GPiw.*ae(GPi(in))'.*blnSTN_in(:,2)'];
% link_t{STN(out)} = zeros(neurons_per_nucleus*2,1); % STN output excitatory
link_n{EXT(out)} = EXT_out;
link_wAMPA{EXT(out)} = [EXT_w_AMPA(in) .* ae_AMPA(SD1(in))' .* blnExt(2,:),...
EXT_w_AMPA(in+neurons_per_nucleus) .* ae_AMPA(SD2(in))' .* blnExt(3,:),...
EXT_w_AMPA(in+2*neurons_per_nucleus) .* ae_AMPA(STN(in))' .* blnExt(1,:)];
link_wNMDA{EXT(out)} = [EXT_w_NMDA(in) .* ae_NMDA(SD1(in))' .* blnExt(2,:),...
EXT_w_NMDA(in+neurons_per_nucleus) .* ae_NMDA(SD2(in))' .* blnExt(3,:),...
EXT_w_NMDA(in+2*neurons_per_nucleus) .* ae_NMDA(STN(in))' .* blnExt(1,:)];
link_wGABAa{EXT(out)} = 0;
link_t{EXT(out)} = zeros(neurons_per_channel*3,1); % no shunting for external drive
% collate synapse types - divide by weight to leave only on
% numerator so that they have effect!
% fast_weights = [SD1_w
% SD2_w
% STN_GPew
% STN_GPiw
% GPe_STNw
% GPe_GPiw
% GPe_GPew
% GPi_GPiw
% EXT_w
% STN_ext_ratio];
if fast_weights(5) ~= 0
max_prox(STN(in)) = max_prox(STN(in)) + (type3 == 2) .* (GPe_STN_w_GABAa(in)' ./ ( abs(fast_weights(5)) .* (1-stnda(2)*dop2) ) ) .* ai_GABAa(STN(in)); % maximum weight value possible
max_soma(STN(in)) = max_soma(STN(in)) + (type3 == 3) .* (GPe_STN_w_GABAa(in)' ./ ( abs(fast_weights(5)) .* (1-stnda(2)*dop2) ) ) .* ai_GABAa(STN(in)); % maximum weight value possible
end
if fast_weights(2) ~= 0
max_prox(GPe(in)) = max_prox(GPe(in)) + (type2 == 2) .* (SD2_w_GABAa(in)' ./ (abs(fast_weights(2)).* (1-gpeda(1)*dop2) ) ) .* ai_GABAa(GPe(in));
max_soma(GPe(in)) = max_soma(GPe(in)) + (type2 == 3) .* (SD2_w_GABAa(in)' ./ (abs(fast_weights(2)) .* (1-gpeda(1)*dop2) ) ) .* ai_GABAa(GPe(in));
end
if fast_weights(1) ~= 0
max_prox(GPi(in)) = max_prox(GPi(in)) + ((type1 == 2) .* (SD1_w_GABAa(in)' ./ abs(fast_weights(1))) .* ai_GABAa(GPi(in)));
max_soma(GPi(in)) = max_soma(GPi(in)) + ((type1 == 3) .* (SD1_w_GABAa(in)' ./ abs(fast_weights(1))) .* ai_GABAa(GPi(in)));
end
if fast_weights(6) ~=0
max_prox(GPi(in)) = max_prox(GPi(in)) + ((type4 == 2) .* (GPe_GPi_w_GABAa(in)'./ abs(fast_weights(6))) .* ai_GABAa(GPi(in)));
max_soma(GPi(in)) = max_soma(GPi(in)) + ((type4 == 3) .* (GPe_GPi_w_GABAa(in)'./ abs(fast_weights(6))) .* ai_GABAa(GPi(in)));
end
%keyboard
% collate synapse types
% pad = zeros(neurons_per_channel*(n_channels-1),1); % use to pad all those for which type is specified on a per channel basis
% max_prox(STN) = max_prox(STN) + [(type3 == 2); pad] .* GPe_STNw' .* ai(STN); % maximum weight value possible
% max_prox(GPe) = max_prox(GPe) + ([(type2 == 2); pad] .* SD2_w' .* ai(GPe)) + ((type5 == 2) .* GPe_GPew' .* ai(GPe));
% max_prox(GPi) = max_prox(GPi) + ([(type1 == 2); pad] .* SD1_w' .* ai(GPi)) + ([(type4 == 2); pad] .* GPe_GPiw' .* ai(GPi)) + ((type6 == 2) .* GPi_GPiw' .* ai(GPi));
%
% max_soma(STN) = max_soma(STN) + [(type3 == 3); pad] .* GPe_STNw' .* ai(STN); % maximum weight value possible
% max_soma(GPe) = max_soma(GPe) + ([(type2 == 3); pad] .* SD2_w' .* ai(GPe)) + ((type5 == 3) .* GPe_GPew' .* ai(GPe));
% max_soma(GPi) = max_soma(GPi) + ([(type1 == 3); pad] .* SD1_w' .* ai(GPi)) + ([(type4 == 3); pad] .* GPe_GPiw' .* ai(GPi)) + ((type6 == 3) .* GPi_GPiw' .* ai(GPi));
end
end
% add diffuse inhibitory output to max prox/soma calculations
GG_sum_type = zeros(3,neurons_per_nucleus);
SS_sum_type = zeros(3,neurons_per_nucleus); % use S prefix fo SNr-SNr collaterals (named GPi-GPi above)
GPe_idx = double(bounds(GPe(1),3))+1:double(bounds(GPe(1),4));
GPi_idx = double(bounds(GPi(1),1))+1:double(bounds(GPi(1),2));
for loop = 1:neurons_per_nucleus
out_type = double(link_t{GPe(loop)}(GPe_idx)); % get output type for non-summed connection
GG_sum_type(1,:) = GG_sum_type(1,:) + (out_type == 1); % distal
GG_sum_type(2,:) = GG_sum_type(2,:) + (out_type == 2); % proximal
GG_sum_type(3,:) = GG_sum_type(3,:) + (out_type == 3); % somatic
out_type = double(link_t{GPi(loop)}(GPi_idx));
SS_sum_type(1,:) = SS_sum_type(1,:) + (out_type == 1); % distal
SS_sum_type(2,:) = SS_sum_type(2,:) + (out_type == 2); % proximal
SS_sum_type(3,:) = SS_sum_type(3,:) + (out_type == 3); % somatic
end
% weights = [SD1_w
% SD2_w
% STN_GPew
% STN_GPiw
% GPe_STNw
% GPe_GPiw
% GPe_GPew
% GPi_GPiw
% EXT_w
% STN_ext_ratio];
if fast_weights(7) ~= 0
max_prox(GPe) = max_prox(GPe) + GG_sum_type(2,:)' .* ai_GABAa(GPe) .* (GPe_GPe_w_GABAa' ./ abs(fast_weights(7))); % NOTE: this will have to change should weights ever be output specific
max_soma(GPe) = max_soma(GPe) + GG_sum_type(3,:)' .* ai_GABAa(GPe) .* (GPe_GPe_w_GABAa' ./ abs(fast_weights(7)));
end
if fast_weights(8) ~= 0
max_prox(GPi) = max_prox(GPi) + SS_sum_type(2,:)' .* ai_GABAa(GPi) .* (GPi_GPi_w_GABAa' ./ abs(fast_weights(8))); % NOTE: this will have to change should weights ever be output specific
max_soma(GPi) = max_soma(GPi) + SS_sum_type(3,:)' .* ai_GABAa(GPi) .* (GPi_GPi_w_GABAa' ./ abs(fast_weights(8)));
end
% scale max_prox and max_soma - to be the sum of simultaneous IPSPs or
% less..
max_prox = max_prox .* max_scale;
max_soma = max_soma .* max_scale;
% remove all zeros - necessary to prevent divide by zero in C code
% is safe because all corresponding link_t{} values are also zero and therefore inhibtory currents
% for these cells never change!
max_prox(max_prox==0) = 1;
max_soma(max_soma==0) = 1;
% find the median value
M_hat = median([max_prox(max_prox<1)' max_soma(max_soma<1)']);
% hist([max_prox(max_prox<1)' max_soma(max_soma<1)'],25)