%modified Kuramoto model % Eric Lowet,2013 ringe-network architecture with sinusodial intrinsic % frequency variation tic clear all number_of_oscillators=160; % --> increasing the oscillator number will increase the number of coupled oscillators... initial_phase= (rand(number_of_oscillators,21)*1*(2*pi)); % Initial phases simulation_time=5 %in sec dt=0.002; %time step (here 2ms) phases = zeros(number_of_oscillators,simulation_time./dt); phases(:,1:21)= initial_phase; %%%%Connectivity matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear dist op=-pi:(2*pi)/(number_of_oscillators-1):pi; for ind2=1:number_of_oscillators dist(ind2,:) = abs( angle( exp(1i*op(ind2) )./exp(1i*op)) ); end dist(dist<0.001) = NaN; % avoid connections with itself s=0.4; %scaling constant for connectivity C=1.65; % strength of connectivity K= (C*exp(-dist./s)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% DEFINING INTRINSIC FREQUENCIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Amplitude=3; Mean_intrinsic_frequency =35; W=(sin((-(1*pi):(2*pi)/((number_of_oscillators./1)-1):(1*pi)))).*Amplitude+Mean_intrinsic_frequency; W=W.*(2*pi); %sum of radians per sec %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear noiseterm %%% pink noise (intrinsic/ dynamical noise) for ind=1:number_of_oscillators %%%%%%%%%%% from 2008.Little MA et al. (2007), "Exploiting nonlinear recurrence and fractal % scaling properties for voice disorder detection", Biomed Eng Online, 6:23 N=(simulation_time./dt);alpha=1; N2 = floor(N/2)-1; f = (2:(N2+1))'; A2 = 1./(f.^(alpha/2)); p2 = (rand(N2,1)-0.5)*2*pi; d2 = A2.*exp(i*p2); d = [1; d2; 1/((N2+2)^alpha); flipud(conj(d2))]; x = real(ifft(d)); %%%%%%%%%%%%%% noiseterm(ind,:)= ((x-mean(x))/std(x))./50; % phase detuning of std=0.02 end %% noise is important to be able to distinguih true synchrony from false %% synchrony (= insufficient dephasing) %%%%%%%%%%%%%% MAIN LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('start simulation') for time=21:simulation_time./dt if mod(time,50) ==0 disp([ num2str(time/(0.001/dt) ) 'ms of ' num2str(simulation_time*1000 ) 'ms']) end for ind=1:number_of_oscillators interact=(sin((phases(ind,time-1) -phases(:,time-1)))) ; phases(ind,time)= phases(ind,time-1) + (dt*W(ind)+ nansum(dt.*K(:,ind).* -interact ) ) + (noiseterm(ind,time)); end end disp('done simulation') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ph= (mod(phases',2*pi))'; %phases xx=(exp(1i*ph(:,100:1:end))'); ab=xx' *xx; allcoh=abs(ab);allcoh=allcoh./max(allcoh(:)); %phase locking matrix allph=angle(ab); %phase relation matrix %% Plotting the phase locking matrix figure('COlor','w','Position',[300 300 300 200]),subplot(1,1,1,'Fontsize',17);imagesc(allcoh) set(gca,'Clim',[ 0 1]) colormap('hot') set(gca,'xticklabel',[],'yticklabel',[]);colorbar axis xy %% Plotting the phase relation matrix coh_thres=0.25; % theshold for illustration figure('COlor','w','Position',[640 300 300 200]), subplot(1,1,1,'Fontsize',17);h=imagesc(allph) acoh=allcoh; tt=(acoh)>coh_thres; % phase-lcoking threshold set(h,'AlphaData',tt ); set(h, 'AlphaDataMapping', 'scaled'); set(gca,'xticklabel',[],'yticklabel',[]) set(gca,'Clim',[ -pi/2 pi/2]) colorbar axis xy toc