% This code computes certain properties of a 'networkbased' segmental
% oscillator model of the lamprey spinal chord proposed by [Williams, 1990].
%
%The following plots are generated:
% 1. the dynamics of an oscillator (Figure 2, top of paper).
%
% 2. The PRC of an oscillator for perturbations of E,L, and C cells (Figure
% 2, bottom of paper). The direct method is used, i.e. the effect of small finite
% perturbations is simulated numerically.
%
% 3. Coupling function of two unidirectionally coupled segmental networkbased oscillators (Figure 6 of paper),
% This is calculated as the integral of
% [presynaptic activity] x [difference between reversion potential and current postsynaptic potential] x [postsynaptic PRC]
% over a period summed for all connections, at given phaseshift between the connected oscillators
%
%4. the stable (decreasing) roots of the coupling functions (Figure 7 of paper)
%
%All calculations are done for several values of tonic drive to E cells,
%which can be specified determined by user.
%levels of tonic drive to E cells
disp('Insert level(s) of tonic drive to E cells ');
disp('please write number in [...] divided by blank spaces');
disp('or press return for default [.005 .0075 .01 .02 .04 .06 .07]');
ves=input('>>>:');
if length(ves)== 0
ves=[.005 .0075 .01 .02 .04 .06 .07];
disp(ves);
end
%ve=[.005 .01 .07];
%%%%%%%%%%%%%%%%%%%%%%%%%determination of coupling functions;
%%%%elements of Fig 2 are plotted during evaluation of @couplingfuncion
for i=1:length(ves)
disp('Current level of tonic drive to E cells is')
disp(ves(i)); disp('Please wait!');
disp('****************************');
[shift{i}, H{i}, fi{i}, PRC{i}, T(i)]=couplingfunction(ves(i),i,length(ves));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%plotting Fig 6
figure(6);
for i=1:length(shift)
plot(shift{i},H{i});
title('***********FIGURE 6*************');
hold on;
end
%%%%%%%%%%%%%%%%%stable zeros of coupling function
disp('calculating roots of copupling functions');
for i=1:length(shift)
for j=1:length(H{i})1
if H{i}(j+1)*H{i}(j)<=0
if H{i}(j+1)>H{i}(j)
stab(i)=(shift{i}(j)*H{i}(j+1)shift{i}(j+1)*H{i}(j))/(H{i}(j+1)H{i}(j))/2/pi;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%plotting Fig 7
figure(7);
plot(1./T,stab,'s');
title('***********FIGURE 7*************');
disp('Ready!');
