# Using ODE for spiking neural networks

13 views (last 30 days)

Show older comments

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?

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

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)));

##### 0 Comments

### Answers (1)

Star Strider
on 7 Aug 2019

##### 0 Comments

### See Also

### Community Treasure Hunt

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

Start Hunting!