% Derrick Roos % % LF_network.m - This function rules the Lansner-Fransen type simulations as % detailed in their paper - "Modelling Hebbian cell assemblies comprised of % cortical neurons". There are an equal number of excitory and inhibitory % cells. Each excitory cell has a companion inhibitory cell which synapses % onto it. The network is trained with different patterns loaded in. Then % its calls a code which produces connection strengths. Positive strengths % form Excitory to Excitory connections. Negative strengths form % connections from Excitory to Inhibitory cells. % % v = LF_network(dt,Tcutoff,Tfin,I0,stim_cells,Patterns) % where: dt = time step % Tcutoff = time cutoff for stimulated cells % Tfin = time end simulation % I0 = stimulation size for cells % stim_cells = vector of cells recieving stimulus % Patterns = matrix where each row is a pattern the network is % trained with % % returns: v = Voltages for excitory cells, indexed as % v(cell,timestep,compartment) % % % ex. p = pattern_maker(40,8,8) % v = LF_network(.01,25,45,1,[1 6 8 10 25],p); function [v,Fsv] = LF_network(dt,Tcutoff,Tfin,I0,stim_cells,Patterns) % find out number of patterns and number of cells from Pattern matrix [n_patterns,n_cells] = size(Patterns); % Patterns = pattern_maker(n_cells,2,2); % Use Tony Kellems Code to produce Wieght matrix, each pattern is given equal signifigance W = make_W(Patterns,ones(n_patterns,1)); % Use D.Roos code to produce G connection/conduction matrix ( splits into new network with % excitory cells and inhibitory) hence all entries in G are positive. G(i,j) represents % strength from cell j to cell i. First half of rows are connections to % excitory cells. Second half of rows are inhibitory cells. Constants % entered represent the ideal average value. See make_G for more detail G = make_G(W,0.01,0.007,0.03); % time that excitory synapse stays on. Time and strength where tested on % p_cell_syn to achieve goal in paper of about 3.6 mV EPSP and 20 ms half % width exc_syn_time = 3; % time that inhibitory synapse stays on. Time and strength where tested on % p_cell_syn to achieve goal in paper of about -5 mV IPSP and 10 ms half % width in_syn_time = 7; % FS-cell parameters, FS - Fast spiking cells, the model used for all % inhibitory cells. Fs.comps = 2; Fs.Gm = 0.0016; Fs.Cm = 0.016; Fs.Gm(2:Fs.comps) = 0.0096; Fs.Cm(2:Fs.comps) = 0.288; Fs.Gcore = 0.0638; Fs.GNa = 1; % mewS Fs.vNa = 50; Fs.vK = -90; Fs.GK = 1; Fs.vCa = 150; Fs.GCa = 0; Fs.GKCA = 0.01; Fs.Gleak = Fs.Gm; Fs.vLeak = -70; Fs.APinflux = 0.013; Fs.APdecay = 20; Fs.v1 = -70; % p-cell parameters, Model used for all excititory cells comps = 4; Gm = 0.0032; %soma Cm = 0.032; % soma Gm(2:comps) = 0.0096; Cm(2:comps) = 0.288; Gcore = 0.04; GNa = 1; % mewS vNa = 40; vK = -70; GK = .5; vCa = 150; GCa = 0; GKCA = 0.0017; Gleak = Gm; vLeak = -50; APinflux = 4; APdecay = 30; v1 = -50; % number of timesteps N = ceil(Tfin/dt); t = zeros(N,1); % allocate space for long vectors % FS - preallocate ion gating variables, voltage and more Fs.v = zeros(n_cells,N,Fs.comps); Fs.n = zeros(n_cells,N); Fs.m = Fs.n; Fs.h = Fs.n; Fs.q = Fs.n; Fs.CaAP = Fs.n; Fs.ISyn = zeros(n_cells,Fs.comps); % p-cells - preallocate ion gating variables, voltage and more v = zeros(n_cells,N,comps); n = zeros(n_cells,N); m = n; h = n; q = n; CaAP = n; Istim = zeros(n_cells,1); ISyn = zeros(n_cells,comps); % find intial points for gating variables Fs.n(:,1) = Fan(Fs.v1)/(Fan(Fs.v1)+Fbn(Fs.v1)); Fs.m(:,1) = Fam(Fs.v1)/(Fam(Fs.v1)+Fbm(Fs.v1)); Fs.h(:,1) = Fah(Fs.v1)/(Fah(Fs.v1)+Fbh(Fs.v1)); Fs.q(:,1) = Faq(Fs.v1)/(Faq(Fs.v1)+Fbq(Fs.v1)); Fs.CaAP(1:n_cells,1)=0; Fs.v(1:n_cells,1,:) = Fs.v1; n(1:n_cells,1) = an(v1)/(an(v1)+bn(v1)); m(1:n_cells,1) = am(v1)/(am(v1)+bm(v1)); h(1:n_cells,1) = ah(v1)/(ah(v1)+bh(v1)); q(1:n_cells,1) = aq(v1)/(aq(v1)+bq(v1)); CaAP(1:n_cells,1)=0; v(1:n_cells,1,:) = v1; % vector containing time remaining for cells synapse to remain on S = zeros(2*n_cells,1); %%% ENTER MAIN LOOP THAT TIMESTEPS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=2:N, t(j) = (j-1)*dt; %%%%%%%% DO FAST SPIKING/ INHIBITORY CELL UPDATE %%%%%%%%%%%%% % Update Soma For Fast Spiking Cells, use hybrid euler scheme Fs.n(:,j) = ( Fs.n(:,j-1) + dt*Fan(Fs.v(:,j-1,1)) )./(1 + dt*(Fan(Fs.v(:,j-1,1))+Fbn(Fs.v(:,j-1,1))) ); Fs.m(:,j) = ( Fs.m(:,j-1) + dt*Fam(Fs.v(:,j-1,1)) )./(1 + dt*(Fam(Fs.v(:,j-1,1))+Fbm(Fs.v(:,j-1,1))) ); Fs.h(:,j) = ( Fs.h(:,j-1) + dt*Fah(Fs.v(:,j-1,1)) )./(1 + dt*(Fah(Fs.v(:,j-1,1))+Fbh(Fs.v(:,j-1,1))) ); Fs.q(:,j) = ( Fs.q(:,j-1) + dt*Faq(Fs.v(:,j-1,1)) )./(1 + dt*(Faq(Fs.v(:,j-1,1))+Fbq(Fs.v(:,j-1,1))) ); Fs.CaAP(:,j) = ( dt*(Fs.vCa - Fs.v(:,j-1,1)).*Fs.APinflux.*Fs.q(:,j).^5 + Fs.CaAP(:,j-1,1) )./( 1+dt*Fs.APdecay ); Fs.top = Fs.Cm(1)*Fs.v(:,j-1,1)+dt*(Fs.Gleak(1)*Fs.vLeak + Fs.GNa*Fs.m(:,j).^3.*Fs.h(:,j)*Fs.vNa + Fs.GK*Fs.n(:,j).^4*Fs.vK + Fs.GCa*Fs.q(:,j).^5*Fs.vCa +Fs.GKCA*Fs.vK*Fs.CaAP(:,j) + Fs.ISyn(:,1)); Fs.bot = Fs.Cm(1) + dt*(Fs.Gleak(1) + Fs.GNa*Fs.m(:,j).^3.*Fs.h(:,j) + Fs.GK*Fs.n(:,j).^4 + Fs.GCa*Fs.q(:,j).^5 + Fs.GKCA*Fs.CaAP(:,j)); % Add compartmental current Fs.top = Fs.top + dt*Fs.Gcore*Fs.v(:,j-1,2); Fs.bot = Fs.bot + dt*Fs.Gcore; Fs.v(:,j,1) = Fs.top./Fs.bot; for k = 2:Fs.comps Fs.top = Fs.Cm(k)*Fs.v(:,j-1,k)+dt*(Fs.Gleak(k)*Fs.vLeak + Fs.ISyn(:,k)); Fs.bot = Fs.Cm(k) + dt*Fs.Gleak(k); % check for right neighbor if k < Fs.comps Fs.top = Fs.top + dt*Fs.Gcore*Fs.v(:,j-1,k+1); Fs.bot = Fs.bot + dt*Fs.Gcore; end % left neighbor compartment Fs.top = Fs.top + dt*Fs.Gcore*Fs.v(:,j-1,k-1); Fs.bot = Fs.bot + dt*Fs.Gcore; Fs.v(:,j,k) = Fs.top./Fs.bot; end %%%%%%%% DO P-CELL/ EXCITITORY CELL UPDATE %%%%%%%%%%%%% % update injected stimulus Istim(stim_cells) = I0*(t(j)>2)*(t(j)