|
- %% CANN + SFA - spiking version
- clear
- tic
-
- rng('shuffle');
- neural_pops_list = {'pyr_popA' 'inh_popA'};
- n_neural_pops = length(neural_pops_list);
-
- pyr_popA.Nneurons = 128; %1024;%512;%128;
- inh_popA.Nneurons = pyr_popA.Nneurons/4;
- norm_factor = 8; %1024/pyr_popA.Nneurons;
-
- t_end = 6000;%10000;%5000;
- delta_t = 0.1;%0.02; %0.05;%0.1;%0.02;
-
- % *********************
- % Membrane Parameters
- % *********************
-
- % ** Pyramidal Cells **
- pyrCellParams.Cm = 0.5; % [nF] membrane capacitance
- pyrCellParams.g_leak = 0.025; % [microSiemens] membrane leak conductance
- pyrCellParams.Vm_leak = -70; % [mV] resting potential
- pyrCellParams.Vm_thresh = -50; % [mV] firing threshold
- pyrCellParams.Vm_reset = -60; % [mV] reset potential
- pyrCellParams.tau_ref = 2; %%2; % [ms] refractory period
-
- % ** Interneurons **
- inhCellParams.Cm = 0.2; % [nF]
- inhCellParams.g_leak = 0.020; % [microS]
- inhCellParams.Vm_leak = -70; % [mV]
- inhCellParams.Vm_thresh = -50; % [mV]
- inhCellParams.Vm_reset = -60; % [mV]
- inhCellParams.tau_ref = 1; % [ms]
-
- pyr_popA.cellParams = pyrCellParams;
- inh_popA.cellParams = inhCellParams;
-
- % Transmission delays: 5
- %pyr_popA.transmission_delays = zeros(pyr_popA.Nneurons,1);
- %inh_popA.transmission_delays = zeros(inh_popA.Nneurons,1);
-
- % *********************
- % Synapese parameters
- % *********************
-
- % ** AMPA **
- synParams.Ve_AMPA = 0; % [mV]
- synParams.tau_s_AMPA = 2; % [ms]
- synParams.alpha_s_AMPA = 1;
-
- % ** X ** NMDA?
- synParams.Ve_X = 0; % [mV]
- synParams.tau_s_X = 30;%40; %%30; % [ms]
- synParams.alpha_s_X = 1;
-
- varying_tauX = 0;
- if varying_tauX == 1 % tauX decreases once every 10 ms
- %dtauX = (synParams.tau_s_X - 10)/(t_end/10); %tau_s_X=100; t_end=5000
- dtauX = (synParams.tau_s_X)/(t_end/10);
- end
-
- % ** GABA **
- synParams.Ve_GABA = -70; % [mV]
- synParams.tau_s_GABA = 10; % [ms]
- synParams.alpha_s_GABA = 1;
-
- % SFA
- synParams.tau_sfaG = 90; %%90; %%80; %[ms]
- synParams.Ve_sfa = -80; %[mV]
- %synParams.d_sfaG = 0.0024;%0.0027;%0.0025;%0.002;%0.001;%0.003; %%0.009; %[microS]
- synParams.d_sfaG = 1/synParams.tau_sfaG*80*0.004;%0.004;%0.005;%0.08; %%0.003; %%0.0025;
- synParams.a_sfaG = 1;%3;%2;%1; %multiplier: Hsfa-2, HHsfa-3
-
- %pyr_popA.sfaG = zeros(pyr_popA.Nneurons,1);
- pyr_popA.sfa_flag = 1;%1;
- %inh_popA.sfa_flag = 0;
- pyr_popA.Wdb_flag = 0;
- conn_plot = 0;
- I_sfa_rec_flag = 1;
-
- inStr_ori = 1;%0.8;%1;%2;%5;%10;%2;%1.3;%0.8; %%0.8; %%0.3; %%1.2; %%1.3;
- stmSpeed = 70;%88; %%70;
-
- pyr_popA.synParams = synParams;
- inh_popA.synParams = synParams;
-
-
- % ****************
- % External Noise 0.8, 0.8
- % ****************
-
- pyr_popA.ext_noise.rate = 0; %%0.1; %0.2;%0.8; %%0.2; %0.15;%0.2;%0.23;%0.25; %%0.8; % [1/ms];
- inh_popA.ext_noise.rate = 0; %%0.1; %0.2;%0.8; %%0.2; %0.15;%0.2;%0.23;%0.25; %%0.8;
-
- % External (noise) to pyramidal cells and interneurons:
- % To adjust spontaneous activities. (20190425)
- pyr_popA.ext_noise.g_AMPA = 0.0031*1;%1.5; % [microS]
- inh_popA.ext_noise.g_AMPA = 0.00238*1;%1.5; % [microS]
-
-
- % ****************
- % Connectivity Profiles
- % Maximum Conductances
- % ****************
-
- J0 = 4;
- sigma_pyr2pyr = 30; %50;%40;%30;%50; %%30; %%15; %tuning curve 30-50
- %popA.J_plus_pyr2pyr = 4.02;
-
- % Difference in preferred angle between two neurons;
- dtheta_pop = 360/pyr_popA.Nneurons;
- % Vector of preferred directions of the neural population:
- pref_dirs_pop = [0:dtheta_pop:(360-dtheta_pop)]';
- pref_dir_diff = min(pref_dirs_pop,360-pref_dirs_pop);
- gauss_con = exp(-0.5*pref_dir_diff.^2/sigma_pyr2pyr^2);
- con_vec_not_normalized = J0/sqrt(2*pi*sigma_pyr2pyr/360)*gauss_con;
- % Normalize the vector:
- con_vec = con_vec_not_normalized/sum(con_vec_not_normalized);
- W_pyr2pyr_tmp = con_vec*1;%1;
-
- if pyr_popA.Wdb_flag == 1
- w = W_pyr2pyr_tmp;
- nbits = 3; %4
- %wstep = (max(w) - min(w))/2^nbits;
- edges = linspace(min(w), max(w), 2^nbits+1);
- values = edges(1:end-1);
- %values = (edges(1:end-1) + edges(2:end))/2;
- digw = discretize(w, edges, values);
- W_pyr2pyr_tmp = digw;
- end
- % unstructured connectvity also leads to bump activity!
- %W_pyr2pyr_tmp = mean(W_pyr2pyr_tmp)*ones(size(W_pyr2pyr_tmp));
- if conn_plot == 1
- figure;
- plot(circshift(W_pyr2pyr_tmp, pyr_popA.Nneurons/2), 'o');
- xlim([0, pyr_popA.Nneurons]);
- end
- popA.W_pyr2pyr_fft = fft(W_pyr2pyr_tmp); %connectivity_profile
-
- popA.W_pyr2inh = 1/pyr_popA.Nneurons;
- popA.W_inh2pyr = 1/inh_popA.Nneurons;
- popA.W_inh2inh = 1/inh_popA.Nneurons;
-
-
- % Maximum Conductances
- inh2pyr_ratio = 3.507; % inh2pyr_ratio=
- rampa = 0.25; %datanum 25
-
- % for inStr_ori = 0.3 (or 1 ---> 0.05)
- popA.G_pyr2pyr_X = norm_factor*0.000381*(1-rampa)*25; %%25; %30;%25;%20; %%20; (32/30, 25, 20 for sigma_pyr2pyr = 50, 30, 15)
- popA.G_pyr2inh_X = norm_factor*0.000292*(1-rampa)*10; %%10;
- popA.G_inh2pyr_GABA = norm_factor*inh2pyr_ratio*0.000381*15;%17;%15; %%15; %20-narrower bump
- popA.G_inh2inh_GABA = norm_factor*inh2pyr_ratio*0.000292;%*1;
-
-
- % ****************
- % Initial Values
- % ****************
-
- pyr_popA.Vm_prev = -51*ones(pyr_popA.Nneurons,1); %-51
- pyr_popA.Vm_new = -51*ones(pyr_popA.Nneurons,1); %-51
- pyr_popA.ext_noise.s_AMPA = 0.6*ones(pyr_popA.Nneurons,1); %0.6
- pyr_popA.gates.s_X = 0.1*ones(pyr_popA.Nneurons,1); %0.1
-
- pyr_popA.sfaG = zeros(pyr_popA.Nneurons,1);
- I_sfa_rec = zeros(pyr_popA.Nneurons, floor(t_end/delta_t)+1);
-
- inh_popA.Vm_prev = -51*ones(inh_popA.Nneurons,1); %-51
- inh_popA.Vm_new = -51*ones(inh_popA.Nneurons,1); %-51
- inh_popA.ext_noise.s_AMPA = 0.6*ones(inh_popA.Nneurons,1); %0.6
- inh_popA.gates.s_GABA = 0.5*ones(inh_popA.Nneurons,1); %0.5
-
- pyr_popA.spikes.WhoFired = [];
- pyr_popA.spikes.SpikeTimes = [];
- pyr_popA.spikes.LastTimeEachFired = -pyr_popA.cellParams.tau_ref(ones(pyr_popA.Nneurons,1));
-
- inh_popA.spikes.WhoFired = [];
- inh_popA.spikes.SpikeTimes = [];
- inh_popA.spikes.LastTimeEachFired = -inh_popA.cellParams.tau_ref(ones(inh_popA.Nneurons,1));
-
- pyr_popA.spikes.Rec = [];
- inh_popA.spikes.Rec = [];
-
- %% ****************
- % Run simulation
- % ****************
-
- for current_time = 0:delta_t:t_end
-
- if mod(current_time,500)==0
- current_time
- end
-
- if varying_tauX == 1 && mod(current_time, 10)==0 && current_time>0
- pyr_popA.synParams.tau_s_X = pyr_popA.synParams.tau_s_X - dtauX;
- pyr_popA.synParams.tau_s_X = max(pyr_popA.synParams.tau_s_X, 30); %30;%1
- end
-
- % ****************
- % pyr_popA
- % ****************
-
- pyr_popA = apply_inputs_changeT1(pyr_popA, current_time, inStr_ori, stmSpeed);
-
- pyr_popA.Vm_prev = pyr_popA.Vm_new;
-
- % Background AMPA current
- rand_vec = rand(pyr_popA.Nneurons,1) < delta_t*pyr_popA.ext_noise.rate;
- pyr_popA.ext_noise.s_AMPA = pyr_popA.ext_noise.s_AMPA + delta_t*(-pyr_popA.ext_noise.s_AMPA/pyr_popA.synParams.tau_s_AMPA) + ...
- pyr_popA.synParams.alpha_s_AMPA*rand_vec;
- pyr_popA.currents.I_noise = pyr_popA.ext_noise.g_AMPA*pyr_popA.ext_noise.s_AMPA.*(pyr_popA.Vm_prev-pyr_popA.synParams.Ve_AMPA);
- % leaky current
- pyr_popA.currents.I_leak = pyrCellParams.g_leak*(pyr_popA.Vm_prev-pyrCellParams.Vm_leak);
-
- % SFA for pyramid neurons
- pyr_popA.sfaG(pyr_popA.spikes.WhoFired) = pyr_popA.sfaG(pyr_popA.spikes.WhoFired) + pyr_popA.synParams.d_sfaG;
- pyr_popA.sfaG = pyr_popA.sfaG + delta_t*(-pyr_popA.sfaG/pyr_popA.synParams.tau_sfaG);
- pyr_popA.currents.I_sfa = pyr_popA.synParams.a_sfaG*pyr_popA.sfaG.*(pyr_popA.Vm_prev - pyr_popA.synParams.Ve_sfa);
- if I_sfa_rec_flag == 1
- I_sfa_rec(:, floor(current_time/delta_t)+1) = pyr_popA.currents.I_sfa;
- end
-
- % recurrent synaptic inputs
- pyr_popA.gates.s_X = pyr_popA.gates.s_X + delta_t*(-pyr_popA.gates.s_X/pyr_popA.synParams.tau_s_X);
- pyr_popA.gates.s_X(pyr_popA.spikes.WhoFired) = pyr_popA.gates.s_X(pyr_popA.spikes.WhoFired) + pyr_popA.synParams.alpha_s_X;
-
- inh_popA.gates.s_GABA = inh_popA.gates.s_GABA + delta_t*(-inh_popA.gates.s_GABA/pyr_popA.synParams.tau_s_GABA);
- inh_popA.gates.s_GABA(inh_popA.spikes.WhoFired) = inh_popA.gates.s_GABA(inh_popA.spikes.WhoFired) + inh_popA.synParams.alpha_s_GABA;
-
- pyr_popA.X_fft = fft(pyr_popA.gates.s_X);
- w_dot_s_X = ifft(popA.W_pyr2pyr_fft.*pyr_popA.X_fft);
- %pyr_popA.currents.I_X = popA.G_pyr2pyr_X*pyr_popA.Nneurons*w_dot_s_X.*(pyr_popA.Vm_prev-pyr_popA.synParams.Ve_X);
- pyr_popA.currents.I_X = popA.G_pyr2pyr_X*w_dot_s_X.*(pyr_popA.Vm_prev-pyr_popA.synParams.Ve_X);
-
- %pyr_popA.currents.I_GABA = popA.G_inh2pyr_GABA*sum(inh_popA.gates.s_GABA)*(pyr_popA.Vm_prev - pyr_popA.synParams.Ve_GABA);
- pyr_popA.currents.I_GABA = popA.G_inh2pyr_GABA*popA.W_inh2pyr*sum(inh_popA.gates.s_GABA)*(pyr_popA.Vm_prev - pyr_popA.synParams.Ve_GABA);
-
- % total input
- if pyr_popA.sfa_flag == 1
- pyr_popA.currents.I_total = - pyr_popA.currents.I_leak - pyr_popA.currents.I_X - pyr_popA.currents.I_sfa...
- - pyr_popA.currents.I_GABA - pyr_popA.currents.I_noise + pyr_popA.currents.I_applied;
- else
- pyr_popA.currents.I_total = - pyr_popA.currents.I_leak - pyr_popA.currents.I_X ...
- - pyr_popA.currents.I_GABA - pyr_popA.currents.I_noise + pyr_popA.currents.I_applied;
- end
-
- % Membrate voltage:
- pyr_popA.Vm_new = pyr_popA.Vm_new + 1/pyrCellParams.Cm*delta_t*pyr_popA.currents.I_total;
- pyr_popA.Vm_new((current_time-pyr_popA.spikes.LastTimeEachFired) < pyrCellParams.tau_ref) = pyrCellParams.Vm_reset;
- pyr_popA.spikes.WhoFired = find(pyr_popA.Vm_new > pyrCellParams.Vm_thresh);
- pyr_popA.spikes.SpikeTimes = current_time + delta_t*((pyrCellParams.Vm_thresh - pyr_popA.Vm_prev(pyr_popA.spikes.WhoFired))./ ...
- (pyr_popA.Vm_new(pyr_popA.spikes.WhoFired) - pyr_popA.Vm_prev(pyr_popA.spikes.WhoFired))); % neural_pop.spikes.SpikeTimes
- pyr_popA.Vm_new(pyr_popA.spikes.WhoFired) = pyrCellParams.Vm_reset;
-
- pyr_popA.spikes.LastTimeEachFired(pyr_popA.spikes.WhoFired) = pyr_popA.spikes.SpikeTimes;
-
- tSpkInfo = [pyr_popA.spikes.WhoFired, pyr_popA.spikes.SpikeTimes];
- pyr_popA.spikes.Rec = [pyr_popA.spikes.Rec; tSpkInfo];
-
- % ****************
- % inh_popA
- % ****************
-
- inh_popA.Vm_prev = inh_popA.Vm_new;
-
- % Background AMPA current
- rand_vec = rand(inh_popA.Nneurons,1) < delta_t*inh_popA.ext_noise.rate;
- inh_popA.ext_noise.s_AMPA = inh_popA.ext_noise.s_AMPA + delta_t*(-inh_popA.ext_noise.s_AMPA/inh_popA.synParams.tau_s_AMPA) + ...
- inh_popA.synParams.alpha_s_AMPA*rand_vec;
- inh_popA.currents.I_noise = inh_popA.ext_noise.g_AMPA*inh_popA.ext_noise.s_AMPA.*(inh_popA.Vm_prev-inh_popA.synParams.Ve_AMPA);
- % leaky current
- inh_popA.currents.I_leak = inhCellParams.g_leak*(inh_popA.Vm_prev-inhCellParams.Vm_leak);
- %recurrent input
- %inh_popA.currents.I_X = popA.G_pyr2inh_X*sum(pyr_popA.gates.s_X)*(inh_popA.Vm_prev - inh_popA.synParams.Ve_X);
- %inh_popA.currents.I_GABA = popA.G_inh2inh_GABA*sum(inh_popA.gates.s_GABA)*(inh_popA.Vm_prev - inh_popA.synParams.Ve_GABA);
- inh_popA.currents.I_X = popA.G_pyr2inh_X*popA.W_pyr2inh*sum(pyr_popA.gates.s_X)*(inh_popA.Vm_prev - inh_popA.synParams.Ve_X);
- inh_popA.currents.I_GABA = popA.G_inh2inh_GABA*popA.W_inh2inh*sum(inh_popA.gates.s_GABA)*(inh_popA.Vm_prev - inh_popA.synParams.Ve_GABA);
-
- % total input
- inh_popA.currents.I_total = - inh_popA.currents.I_leak - inh_popA.currents.I_X - inh_popA.currents.I_GABA - inh_popA.currents.I_noise;
-
- % Membrate voltage:
- inh_popA.Vm_new = inh_popA.Vm_new + 1/inhCellParams.Cm*delta_t*inh_popA.currents.I_total;
- inh_popA.Vm_new((current_time-inh_popA.spikes.LastTimeEachFired) < inhCellParams.tau_ref) = inhCellParams.Vm_reset;
- inh_popA.spikes.WhoFired = find(inh_popA.Vm_new > inhCellParams.Vm_thresh);
- inh_popA.spikes.SpikeTimes = current_time + delta_t*((inhCellParams.Vm_thresh - inh_popA.Vm_prev(inh_popA.spikes.WhoFired))./ ...
- (inh_popA.Vm_new(inh_popA.spikes.WhoFired) - inh_popA.Vm_prev(inh_popA.spikes.WhoFired))); % neural_pop.spikes.SpikeTimes
- inh_popA.Vm_new(inh_popA.spikes.WhoFired) = inhCellParams.Vm_reset;
-
- inh_popA.spikes.LastTimeEachFired(inh_popA.spikes.WhoFired) = inh_popA.spikes.SpikeTimes;
-
- tSpkInfo = [inh_popA.spikes.WhoFired, inh_popA.spikes.SpikeTimes];
- inh_popA.spikes.Rec = [inh_popA.spikes.Rec; tSpkInfo];
-
- end
-
- save('cann_spk_data');
-
- toc
-
-
-
- %% raster plot
- %t_end = 5000; pyr_popA = pyr_popA_rec{5};
- fireA = pyr_popA.spikes.Rec;
-
- nA = zeros(t_end+1,1);
- time_window = 50;
- for i = 0:t_end
- delta_t = min(t_end, i+time_window/2)- max(0, i-time_window/2);
- tnA = length(find(fireA(:,2) > i-time_window/2 & fireA(:,2) <= i+time_window/2)) / delta_t / pyr_popA.Nneurons * 1000;
- nA(i+1) = tnA;
- end
-
- figure;
- subplot(211)
- plot(fireA(:, 2), fireA(:, 1), '.k'); grid on;
- ylim([0, pyr_popA.Nneurons]); xlim([0, t_end]);
- ylabel('Neuron ID'); xlabel('Time (ms)');
- subplot(212);
- b=0:t_end;
- plot(b,nA,'k');
|