/
solveODEreduc.m
142 lines (119 loc) · 3.76 KB
/
solveODEreduc.m
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
136
137
138
139
function [synchr_index, tConv, in_freq_absPow,delayPh]=solveODEreduc(plotflag, par)
% input set of parameters:
% par = [P w1 w2 w3 w4 w5 w6 w7 q AMP FREQ phase]
%
% output:
% synchronization index
% time needed for convergence
% absolute power of input frequency
% phase delay
% Settings for the solver
options = odeset('InitialStep',1e-03,'MaxStep',0.005);
% Define simulation time
t_end=20;
timestep=0.002;
tRange= 0 : timestep : t_end-timestep;
Fs=1/timestep;
x_ini= [0 0 0];
phase = par(12);
in_freq= par(11);
ampl=par(10);
offset= par(1);
thalf = t_end/2; % first half of simulation where no transition allowed
%prepare constant drive - parameter offset used
constantDrive = offset * ones(size(tRange));
% simulate with a constant input first
% the purpose is to find the transition time: when the constant input
% should become oscillating input
parameters={tRange, constantDrive,par(2),par(3),par(4),par(5),par(6),par(7),par(8),par(9)};
[t,x]=ode45(@modelFunqVarInput,tRange,x_ini,options,parameters);
% find the peaks in both the signal and its first derivative
% the peaks represent the start of phases 0, pi/2, pi, 3pi/2
X = x(:,1);
XX = [0; diff(X)];
switch phase
case 0
[~,IND] = findpeaks(XX);
case 1
[~,IND] = findpeaks(X);
case 2
[~,IND] = findpeaks(-XX);
case 3
[~,IND] = findpeaks(-X);
otherwise
error('check phase value')
end
clear X XX
count=1;
for i=length(IND)-1:-1:1 % cleaning up the phase indices
if IND(i)==IND(i+1)-count
IND(i)=[];
count= count+1;
else
count=1;
end
end
transIND = IND(find(t(IND)>thalf,1)); % transition index, first after thalf
transitionT = t(transIND); % transition time
% prepare the oscillatory input - starts at transitionT with phase 0
oscillDrive=offset+ampl*sin(2*pi*in_freq*tRange);
oscillDrive = circshift(oscillDrive,transIND);
oscillDrive(1:transIND) = offset*ones(1,transIND);
% simulate again , now with oscillatory input
parameters={tRange, oscillDrive,par(2),par(3),par(4),par(5),par(6),par(7),par(8),par(9)};
[t,y]=ode45(@modelFunqVarInput,tRange,x_ini,options,parameters);
% make up the attractor set to detect convergence in all 3 dimensions
tolerance = 1e-3;
indP = round(17*size(y,1)/20);
attrSet = uniquetol(y(indP:end,:), tolerance,'ByRows',true);
% find the timpoint of convergence
for j = indP :-1 :1
if ~ismembertol(y(j, :), attrSet, 10*tolerance, 'ByRows',true)
convergedAt = tRange(j+1);
break
end
end
% time needed for convergence
tConv = convergedAt - transitionT;
tConv(tConv<0) = 0;
% find the delay of between y1 and oscillDrive, with oscillDrive being the
% reference
delay = finddelay(oscillDrive(indP:end), y(indP:end,1)); % this is in timesteps
delayT = delay*timestep;
% delay in phase, phase difference
delayPh = 2*pi*(delayT*in_freq);
if delayPh < -pi
delayPh = delayPh + 2*pi;
end
if delayPh > pi
delayPh = delayPh - 2*pi;
end
% plot the response of the system and indicate the timepoint of
% convergence
if plotflag
figure(1)
subplot(2,1,1)
plot(t,oscillDrive);
hold on
stem(convergedAt, 0.5)
legend('external drive P(t)', 'convergence timepoint')
ylabel('P(t)')
subplot(2,1,2)
plot(t,y(:,1))
hold on
stem(convergedAt, 0.5)
legend( 'excitatory population E(t)', 'convergence timepoint')
xlabel('time (s)')
ylabel('E(t)')
hold off
end
% run spectral analysis on the last third of system response
ypart = y(round(2*end/3):end,1); %last third of system response
[out_freq, powers, in_freq_absPow]=freq_analys(ypart,Fs,plotflag,in_freq);
synchr_index=0;
for i=1:length(out_freq)
if out_freq(i)>=in_freq-0.2 && out_freq(i)<=in_freq+0.2
synchr_index=powers(i);
break
end
end