Using ODE for spiking neural networks

13 views (last 30 days)
shani maoz
shani maoz on 7 Aug 2019
Answered: Star Strider on 7 Aug 2019
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?
Many thanks in advance! :)
% 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
idexs = [];
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;
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:
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)));

Answers (1)

Star Strider
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!