Resonate-and-Fire (RAF) Neuron Model Cortical Modelling Assignment

%% Simulator - Resonate-and-Fire (RAF) Neuron Model 
%% D. Bernal-Casas 
%% Version 10: 2019 - Adapted from the literature 
%% Main Sources: 
%% [1] Trappenberg  
%% [2] Dayan and Abbot 
%% [3] Churchland and Sejnowski 
%% [4] Rolls and Deco 
%% [5] Izhikevich 
%% [6] Wilson 
%% [7] Gerstner 
 
%% Numerical Integration with the Euler method 
 
clear all; 
close all; 
hold off;  
clc; 
 
%% Setting parameters 
 
%Samples 
Nsamp = 300; 
%Damping factor 
damp_1 = 2.5; 
%Natural frequency, make oscillations of 1.1Hz 
w_1 = 1.1*2*pi; 
%Time-step 
dt = 0.1;  
 
times = 0:dt:(Nsamp*dt - dt); 
 
%% Set up neurons 
 
dv_1 = zeros(2, Nsamp); 
v_1 = zeros(2, Nsamp); 
 
dv_2 = zeros(2, Nsamp); 
v_2 = zeros(2, Nsamp); 
 
dv_3 = zeros(2, Nsamp); 
v_3 = zeros(2, Nsamp); 
 
%% Initialize weights 
%Pos weight 
n1 = randi([1 20],1,1); 
n2 = randi([1 20],1,1); 
n3 = randi([1 20],1,1); 
n4 = randi([1 20],1,1); 
 
 
 
%Neg weight 
n5 = (randi([1 20],1,1))*-1; 
n6 = (randi([1 20],1,1))*-1; 
 
%Excitatory Neuron one to three and one to two excites: 
w12 = n1; 
w13 = n2; 
 
%Excitatory Neuron two to one excites and two to three excites: 
w21 = n3; 
w23 = n4; 
 
%Inhibitory: Neuron three to one inhibits and three to two inhibits:  
w31 = n5; 
w32 = n6; 
 
 
 
 
 
%% Inputs 
I_0 = zeros(1, Nsamp); 
%With a single input, the neuron does not fire 
%I_0(1,10) = 37; 
%With two-consecutive inputs, you have a coincidence detector 
%I_0(1,10) = 37; 
% I_0(1,11) = 37; 
 
spike_1 = zeros(1, Nsamp); 
spike_2 = zeros(1, Nsamp); 
spike_3 = zeros(1, Nsamp); 
 
I_0 = zeros(1, Nsamp); 
for in = 10:10:Nsamp 
I_0(1,in) = 37; 
%if in < 290 
%I_0(1,in+1) = 37; 
%end 
end 
 
 
 
 
 
%Neuron 1 
 
for it=2:Nsamp 
%1 Neuron 
dv_1(1,it) = v_1(2,it-1); 
dv_1(2,it) = -(w_1^2)*v_1(1,it-1) - 2*damp_1*v_1(2,it-1) + I_0(1,it-1) + w21*spike_2(1,it-1)+w31*spike_3(1,it-1); 
v_1(1,it) = v_1(1,it-1) + dt*dv_1(1,it); 
v_1(2,it) = v_1(2,it-1) + dt*dv_1(2,it); 
if (v_1(1,it) >= 1.0) 
spike_1(1,it) = 1; 
else 
spike_1(1,it) = 0; 
end 
 
 
 
%Neuron 2 
 
 
 
%2 Neuron 
dv_2(1,it) = v_2(2,it-1); 
dv_2(2,it) = -(w_1^2)*v_2(1,it-1) - 2*damp_1*v_2(2,it-1) + I_0(1,it-1) + w12*spike_1(1,it-1)+w32*spike_3(1,it-1); 
v_2(1,it) = v_2(1,it-1) + dt*dv_1(1,it); 
v_2(2,it) = v_2(2,it-1) + dt*dv_1(2,it); 
if (v_2(1,it) >= 1.0) 
spike_2(1,it) = 1; 
else 
spike_2(1,it) = 0; 
end 
 
 
 
%Neuron 3 
 
 
 
%3 Neuron 
dv_3(1,it) = v_3(2,it-1); 
dv_3(2,it) = -(w_1^2)*v_3(1,it-1) - 2*damp_1*v_3(2,it-1) + I_0(1,it-1) + w23*spike_2(1,it-1)+w13*spike_1(1,it-1); 
v_3(1,it) = v_3(1,it-1) + dt*dv_3(1,it); 
v_3(2,it) = v_3(2,it-1) + dt*dv_3(2,it); 
if (v_3(1,it) >= 1.0) 
spike_3(1,it) = 1; 
else 
spike_3(1,it) = 0; 
end 
end 
 
 
 
%% 
figure('name','Excitatory neuron 1, Spikes, Membrane Potential, and Input Current') 
 
subplot(3, 1, 1) 
stem(times, spike_1(1,:)) 
 
subplot(3, 1, 2) 
plot(times, v_1(1,:)) 
 
hold on 
plot(times, ones(1,Nsamp), 'r') 
subplot(3, 1, 3) 
stem(times, I_0(1,:)) 
 
%%  
 
figure('name','Excitatory neuron 2, Spikes, Membrane Potential, and Input Current') 
title('Excitatory Neuron 2') 
 
subplot(3, 1, 1) 
stem(times, spike_2(1,:)) 
 
subplot(3, 1, 2) 
plot(times, v_2(1,:)) 
 
hold on 
plot(times, ones(1,Nsamp), 'g') 
subplot(3, 1, 3) 
stem(times, I_0(1,:)) 
 
%% 
figure('name','Inhibitory neuron 1, Spikes, Membrane Potential, and Input Current') 
 
subplot(3, 1, 1) 
stem(times, spike_3(1,:)) 
 
subplot(3, 1, 2) 
plot(times, v_3(1,:)) 
 
hold on 
plot(times, ones(1,Nsamp), 'b') 
subplot(3, 1, 3) 
stem(times, I_0(1,:)) 
 
%% 
figure('name','Combined interactions, Spikes, Membrane Potential, and Input Current') 
subplot(3, 1, 1) 
hold on 
stem(times, spike_1(1,:),'r') 
stem(times, spike_2(1,:),'g') 
stem(times, spike_3(1,:),'b') 
subplot(3, 1, 2) 
 
hold on 
plot(times, v_1(1,:)) 
plot(times, v_2(1,:)) 
plot(times, v_3(1,:)) 
 
hold on 
plot(times, ones(1,Nsamp), 'r') 
 
subplot(3, 1, 3) 
stem(times, I_0(1,:)) 
%stem(times, I_0(1,:)) 
%stem(times, I_0(1,:))