Repository navigation
Expand file tree
/
Copy pathDISINHIBITION_MODEL_2.m
More file actions
135 lines (109 loc) · 4.62 KB
/
Copy pathDISINHIBITION_MODEL_2.m
File metadata and controls
135 lines (109 loc) · 4.62 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
%% disinhibition model 2
clear all
clc
%close all
rng(2) %fix RNG just for refersion interneurons down the road (see documentation), user can change this otherwise.
%% Neuronal Parameters
NE = 200; %number of excitatory neurons
NI = 200; %number of inhibitory neurons
mu = NE/(NE+NI); %ratio of E to total
N = NE + NI; %total number of neurons
T = 900; %total simulation time (s)
dt = 0.00005; %time step (s)
nt = round(T/dt);
tref = 0.002; %Refractory time constant in seconds
tm = 0.01; %Membrane time constant
vreset = -65; %Voltage reset
vpeak = -40; %Voltage peak.
td = 0.02; %decay time;
tr = 0.002; %rise time;
nq = 20; %only store every nq time steps
tdo = 0.8; %decay time of synthetic calcium
tro = 0.02; %rise time of synthetic calcium
sig = 0;
%% Coupling Weight Matrix
disp('Initializing Weight Matrix')
G = 0.05; %Static Weight Magnitude
p = 0.1; %sparsity of inhibition
CE = round(p*NE); %Number of E connections
CI = round(p*NI); %Number of I connectiosn
OMEGA(:,NE+1:N) = -G*(2*rand(N,NI)).*(rand(N,NI)<p); %/sqrt(p*NI);
OMEGA(NE+1:N,1:NE) = 0.5*(rand(NI,NE)<0.01);
%% initialization parameters for the network
IPSC = zeros(N,1); %post synaptic current storage variable
h = zeros(N,1); %Storage variable for filtered firing rates
r = zeros(N,1); %second storage variable for filtered rates
hr = zeros(N,1); %Third variable for filtered rates
ro = r; %used for integrating calcium traces
hro = hr; %used for integrating calcium traces,
JD = 0*IPSC; %storage variable required for integrating connections
tspike = zeros(10*nt,2); %store spike times in tspike
ns = 0; %Number of spikes, counts during simulation
v = vreset + rand(N,1)*(10); %Initialize neuronal voltage with random distribtuions
v_ = v; %v_ is the voltage at previous time steps
tlast = zeros(N,1); %variable sed to implement refractory period
kd = 0; %used to initialize storage every nq time steps, kd = kd + 1 if mod(i,nq)==1;
nq = 20; %store every nq time steps
REC2 = zeros(nt/nq,N); %store calcium traces in REC2
REC = zeros(nt,20); %store voltage traces in REC
BIAS(1:NE,1)= -35;
BIAS(NE+1:N,1)= -25; %Background current to Inhibitory neurons
sig = 0.0; %white noise standard deviation
idstore = randperm(N,20); %store a random subset of neurons
%% set input weights for the burst/home-cage/novel environment signal
WNV = -7*(ones(N,1));
WNV(1:NE,:) = 0;
NOV = 1; %amount of increase in current due to novel environment
%% generate bursts during home-cage sessions
tcrit1 = 300; %transfer from home-cage to novel environment
tcrit2 = 600; %transfer from novel-environment back to home-cage
dfs = 0.05; %approximate burst duration
time = ((1:nt)*dt)';
dtimes = [tcrit1*rand(10,1);tcrit2 + (T-tcrit2)*rand(10,1)]; %generate random burst times
dtimes = sort(dtimes);
I_INT = zeros(nt,1); %burst storage variable
for i = 1:length(dtimes)
I_INT = I_INT + exp(-((time - dtimes(i)).^2)/dfs);
end
I_INT = 2.5*I_INT;
%% START INTEGRATION
for i = 1:1:nt
I = IPSC + BIAS + WNV*NOV*(dt*i>tcrit1)*(dt*i<tcrit2) + WNV*I_INT(i); %implement + home cage changes
dv = (dt*i>tlast + tref).*( (-v+I)/tm) ; %Voltage equation with refractory period
v = v + dt*(dv) + sqrt(dt)*randn(N,1)*sig/sqrt(tm);
index = find(v>=vpeak); %Find the neurons that have spiked
%Store spike times, and get the weight matrix column sum of spikers
if length(index)>0
JD = sum(OMEGA(:,index),2); %compute the increase in current due to spiking
tspike(ns+1:ns+length(index),:) = [index,0*index+dt*i]; %store spike times
ns = ns + length(index); % total number of spikes so far
end
tlast = tlast + (dt*i -tlast).*(v>=vpeak); %Used to set the refractory period of LIF neurons
% Code if the rise time is 0, and if the rise time is positive
if tr == 0
IPSC = IPSC*exp(-dt/td)+ JD*(length(index)>0)/(td);
r = r *exp(-dt/td) + (v>=vpeak)/td;
else
IPSC = IPSC*exp(-dt/tr) + h*dt;
h = h*exp(-dt/td) + JD*(length(index)>0)/(tr*td); %Integrate the current
r = r*exp(-dt/tr) + hr*dt;
hr = hr*exp(-dt/td) + (v>=vpeak)/(tr*td);
%%integrate calcium resposnes
ro = ro*exp(-dt/tro) + hro*dt ;
hro = hro*exp(-dt/tdo) + (v>=vpeak)/(tro*tdo);
end
v = v + (30 - v).*(v>=vpeak); %rest the voltage and apply a cosmetic spike.
REC(i,:) = [v(1:10);v(NE+1:NE+10)];
v = v + (vreset - v).*(v>=vpeak); %reset spike time
% Store filtered spike trains for 10 neurons.
if mod(i,nq)==1
kd = kd + 1;
REC2(kd,:) = ro;
end
%% plotting results
if mod(i,round(10/dt))==1
prog = dt*i/T
end
end
%%
save disinhibition_model_2.mat -v7.3