# Using ODE for spiking neural networks

13 views (last 30 days)
shani maoz on 7 Aug 2019
Answered: Star Strider on 7 Aug 2019
Hi,
I'm trying to solve a set of vectorial differential equations, aiming to simulate a spiking neural network (meaning that I have some variable which represents a voltage - V, and as it crosses some threshold a dirac delta function is added to another vector - S, which represents the "spikes" of the neurons.
I don't know how to do it with the ODE of matlab, so I have tried the regular Euler method, with a descrete delta correspondig to 1/dt... I'm not sure why, but it's not working for me and I suspect that this method is not accurate enogh... (meaning that this code is working, but the outcome is not good enogh)
can anybody help me convert this code to ODE?
% the network contains N neurons
% threshold = some pre- dfined scalar
% A = a matrix generaly, here it is 0.
% OmegaSlow,OmegaFast - N by N matrixes of connectivity between different neurons.
% c - vector of external input, it's length is the length of tspan, the time of the experiment.
% OutpMatrix - in this simulation it's an N by 1 vector. defines some weights for the external input c.
% lambda_d, lambda_v = pre-defined scalars.
% Run experiment:
%Initial conditions are all 0.
x(:,1) = Xvec0; % x - some hidden time changing variable which the network is trying to follow without "knowing" explicitly.
r(:,1) = fireRate0; % xHat - the final output of the network. should track x.
xHat(:,1) = xHat0;
V(:,1) = V0; % voltage
%timesteps and experiment time:
dt = 0.1*10^-3;
tspan = 0:dt:1;
idexs = [];
Svec = zeros(N,length(tspan));%the name of this variable is confusing. actually, each colomn in this matrix is the spins vector.
for i = 1:length(tspan)-1
x(:,i+1) = x(:,i) + dt*(A*x(:,i) + c(:,i));
xHat(:,i+1) = xHat(:,i) + dt*(-lambda_d*xHat(:,i) + OutpMatrix*Svec(:,i));
r(:,i+1) = r(:,i) + dt*(lambda_d*(-r(:,i) + Svec(:,i)));
V(:,i+1) = V(:,i) + dt*(-lambda_v*V(:,i) + (1/lambda_d)*OmegaSlow*(r(:,i)) - OmegaFast*Svec(:,i) + OutpMatrix'*c(:,i) + noise(:,i));
idexs = find(V(:,i+1) > threshold);
if length(idexs) >= 1
tau = tspan(i+1);
for idx = 1:length(idexs)
Ovec(idexs(idx),:) = dirac(tspan - tau);
j = Svec == Inf; % find Inf
Svec(j) = 1/dt; % set Inf to finite value
end
idexs = [];
end
end
In order to understand the idea this code is enogh, but in order to run it you need the exact parameters, which I'm adding in the code below. I may also add my outcome: %defining parameters of the system:
N = 1000; J = 1; lambda_d = 10; lambda_v = 20; miu = 0.2/N^2; ni = 0.7/N^2; A = 0; OutpMatrix = randn(J,N);
%normalizing matrix OutpMatrix
for n=1:size(OutpMatrix,2)
vector = OutpMatrix(:,n);
normalization = 1/(N*sqrt(vector'*vector));
OutpMatrix(:,n) = normalization*vector;
end
I1 = eye(J,J); I2 = eye(N,N);
OmegaSlow = OutpMatrix'*(A+lambda_d*I1)*OutpMatrix;
OmegaFast = OutpMatrix'*OutpMatrix +miu*(lambda_d)^2*I2;
%external input:
c=(N/2*(heaviside(tspan-0.5)-heaviside(tspan-0.8))-N/2*(heaviside(tspan-1)-heaviside(tspan-1.5)));
sigma_v = 10^-3;%noise level
noise = sigma_v*normrnd(0,1/sqrt(dt),[J,length(tspan)]);%noise shape.
%initial conditions:
V0 = zeros(N,1); fireRate0 = zeros(N,1); Xvec0 = 0; xHat0 = Xvec0;
%definition of the threshold:
threshold = (1/2)*(ni*lambda_d + miu*(lambda_d)^2 + sqrt(OutpMatrix(:,1)'*OutpMatrix(:,1)));

Star Strider on 7 Aug 2019
See The Hodgkin-Huxley Model - McGill University paper for a clear description of the various ways to solve such problems, complete with MATLAB code.