% function dnf
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% Dynamic Neural Field Model (1D)
%
% G. Blohm (2010)
% adapted from Trappenberg (2010)
%-------------------------------------------------------------------------

clear; 

%% initialize
nn = 100; 
dx = 2*pi/nn; % node distance
sig = 2*pi/10; % node SD
C = 0.5; % weight inhibition constant
tall = []; % time vector to fill
rall = []; % activation matrix to fill

%% generate weight matrix (center-surround)
for loc = 1:nn;
    i = (1:nn)'; 
    dis = min(abs(i-loc),nn-abs(i-loc)); % distance between nodes
    pat(:,loc) = exp(-(dis*dx).^2/(2*sig^2)); % connection weight between nodes
end
w = pat*pat'; % bi-directional connection weights
w = w/w(1,1); % normalize
w = 4*(w-C); % shift Gaussian weights: center excitation - sourround inhibition

%% Update with localised input (I_ext >0)
I_ext = zeros(nn,1); 
I_ext(nn/2-floor(nn/10):nn/2+floor(nn/10)) = 1; % set input for 2*nn/10 neurons at center (nn/2) to 1 (other inputs are 0)
[t,u] = ode45('rnn_ode',[0 10],zeros(1,nn),[],nn,dx,w,I_ext); % Matlab numerical integration of differential equations in rnn_ode.m
r = 1./(1+exp(-u)); % network output after passing through sigmoid
tall = [tall;t]; 
rall = [rall;r];

%% Update without input (I_ext = 0): system relaxes here
I_ext = zeros(nn,1);
[t,u] = ode45('rnn_ode',[10 20],u(size(u,1),:),[],nn,dx,w,I_ext);
r = 1./(1+exp(-u)); 
tall = [tall;t]; 
rall = [rall;r];

%% Plotting results
figure
clf; 
hold on;
surf(tall',1:nn,rall','linestyle','none'); 
xlabel('time')
ylabel('node #')
zlabel('activity')
view(45,45);