function ncmodel()
 
%  ---------------------- Computational model  ------------------------------------------------------
%
%   Code author: Chao Li
%   Date:      11/23/2014
%   Please cite: Saha D., Li C. et al. Behavioural correlates of combinatorial versus
%                temporal features of odour codes. Nat. Commun. 6:6953 doi: 10.1038/ncomms7953 (2015). 
%   This code passed the test on MATLAB 8.3.0.532 (R2014a) installed with Statistics Toolbox
 
%   Running this code will generate five figures that correspond to
%   all the panels of Fig. 6 in the main manuscript
%   the figure(1) : Fig. 6a
%   the figure(2) : Fig. 6c left panel
%   the figure(3) : Fig. 6b
%   the figure(4) : Fig. 6c right panel
%   the figure(5) : Fig. 6d
%
%   Note that since several parameters are obtained by random number,
%   the figures created here may not be exactly the same as in Fig.6.
%  -------------------------------------------------------------------------------------------------
    clear,  clc,  close all;
 
%% Parameter set obtained from physiological data to constrain the model 
    
    % time-to-peak distributions derived from real data {see Supplementary Fig. 3
    %  hex, hex(hex); 2oct, 2oct(hex) }
    scale_cnt = 100; % scale constant
    % hex
    peakIndex_hex=scale_cnt*[7,6,7,9,8,10,5,5,6,7,9,7,6,12,7,8,7,6,11,11,25,9,5,7,9,44,... 
        26,18,8,9,8,7,10,8,6,9,11,9,9,10,6,13,7,11,13,7,9,14,10,25,6,7,7,6,6,18,24,5];
    % hex(hex)    
    peakIndex_hexhex=scale_cnt*[3,50,54,68,10,7,23,11,72,52,42,56,36,10,6,39,10,20,8,16,51,... 
        6,30,36,18,15,10,11,41,10,9,27,52,8,42,50,31,49,53,29,31,57,34,6,8,77,];
    % 2oct
    peakIndex_oct = scale_cnt*[8,34,13,8,6,7,7,12,6,6,8,8,11,12,15,6,13,6,7,7,6,7,20,7,... 
        6,9,9,8,10,8,10,24,8,14,7,20,24,5,16,3,5,9,6,6,7,8];
    % 2oct(hex)
    peakIndex_octhex = scale_cnt*[35,74,36,21,25,16,11,50,35,47,53,15,15,17,76,10,73,53,28,...
        7,68,10,9,68,8,21,16,11,16,63,60,39];
    
    
%  -------------------------------------------
 
    
%% define alpha function to characterize neural response profile
    t = 0.1:.1:1000; % time point series
    neuron_number = 100; % number of neurons
    
    spikeActivity_1= []; % spiking activity by Odor 1
    spikeActivity_2= []; % spiking activity by Odor 2
    
    activationSet_1 = 65; %  65% of PNs activated by odor 1 (hex)
    activationSet_2 = 50; %  50%             .... by odor 2 (2oct)
    
    Amplitude_1 = [exprnd(10,[1,activationSet_1]), zeros(1,100-activationSet_1)]; % The distributed Amplitude value for odor 1
    Amplitude_2 = [zeros(1,30),exprnd(10,[1,activationSet_2]),zeros(1,20)]; % ............. for odor 2
    
    tau = ceil(rand(1,1000).*250); % time constant of each response by odor 1
    tau2 = ceil(rand(1,1000).*250); % ....... by odor 2
    
    % alpha response shape
    for i = 1: neuron_number
        spikeActivity_1(i,:) = (t.*exp(-t.^1.25/tau(i)));
        spikeActivity_1(i,:) = (spikeActivity_1(i,:)./max(spikeActivity_1(i,:))).*Amplitude_1(i);
        
        spikeActivity_2(i,:) = (t.*exp(-t.^1.25/tau2(i)));
        spikeActivity_2(i,:) = (spikeActivity_2(i,:)./max(spikeActivity_2(i,:))).*Amplitude_2(i);
    end
    
    spikeActivity_1 = [zeros(100,2000),spikeActivity_1]; % add baseline period
    spikeActivity_2 = [zeros(100,2000),spikeActivity_2]; % add baseline period
%  -------------------------------------------------
 
 
%% Simluate PN response for Odor 1 and jittered version of Odor 1
    duration = size(spikeActivity_1,2);
    sampleRate = 2400; % sampling rate 
    timelength = round(duration/sampleRate);
    
    
    % add jitter for the peak position
    % for hex condition
    jitter_hex = generateJitter(peakIndex_hex, neuron_number); % generate jitter level
    [neuralActivity_odor1, peakIndex_1] = addJitter(spikeActivity_1,jitter_hex,neuron_number); % add jitter to neural activity
    peakIndex_1 = peakIndex_1/sampleRate; % convert time into second
    % for hex(hex) condition
    jitter_hexhex = generateJitter(peakIndex_hexhex, neuron_number); % generate jitter level
    [neuralActivity_odor1_jitter, peakIndex_1_jitter] = addJitter(spikeActivity_1,jitter_hexhex,neuron_number);% add jitter to neural activity
    peakIndex_1_jitter = peakIndex_1_jitter/sampleRate; % convert time into second
    
    % -----  figure plot -------------------------------------------------
    
    % color parameter
    color_PSTH = [  1 0 0;  0 0 0]; % color for PSTH
    color_Shade = [1 .8 .8; .8 .8 .8]; % shade color
    color_map = redblue(100);
    color = color_map(1*size(color_map,1)/2:end,:);
    
    % plot Fig. 6a
    timex = 1/sampleRate:1/sampleRate:timelength; % time series in second
    figure(1); 
    subplot(2,2,1);imagesc( timex ,1:100,neuralActivity_odor1); colormap(color); caxis([0 25]); 
    title 'Odor 1'; xlabel 'Time (s)'; ylabel 'Cell index'
    subplot(2,2,2);imagesc( timex ,1:100, neuralActivity_odor1_jitter); colormap(color); caxis([0 25]); 
    title 'Jittered Odor 1'; xlabel 'Time (s)'; ylabel 'Cell index'
    
    tt = [0:600/sampleRate:12000/sampleRate]; % time bin center for histogram
    distribution_1 = hist(peakIndex_1,tt); 
    distribution_1_jitter = hist(peakIndex_1_jitter,tt);
    distribution_1 = distribution_1./activationSet_1; 
    distribution_1_jitter = distribution_1_jitter./activationSet_1;
    subplot(2,2,3);
    plot(tt,distribution_1,'r-o','MarkerFaceColor',[1 0 0]); axis([0 12000/sampleRate 0 .6]); 
    title 'Peak time distribution'; xlabel 'Time to peak'; ylabel 'Probability'
    subplot(2,2,4);
    plot(tt,distribution_1_jitter,'r-o','MarkerFaceColor',[.3 0 0]); axis([0 12000/sampleRate 0 .6]);
    title 'Peak time distribution'; xlabel 'Time to peak'; ylabel 'Probability'
    
    % plot Fig. 6c left panel
    figure(2);     
    fill([timex,timex(end:-1:1)],[mean(neuralActivity_odor1)-std(neuralActivity_odor1)/sqrt(100),...
        mean(neuralActivity_odor1(:,duration:-1:1))+std(neuralActivity_odor1(:,duration:-1:1))/sqrt(100)],color_Shade(1,:),'lineStyle','none'); hold on;
    fill([timex,timex(end:-1:1)],[mean(neuralActivity_odor1_jitter)-std(neuralActivity_odor1_jitter)/sqrt(100),...
        mean(neuralActivity_odor1_jitter(:,duration:-1:1))+std(neuralActivity_odor1_jitter(:,duration:-1:1))/sqrt(100)],color_Shade(2,:),'lineStyle','none');
    plot(timex,mean(neuralActivity_odor1),'color',color_PSTH(1,:),'linewidth',2);
    plot(timex,mean(neuralActivity_odor1_jitter),'color',color_PSTH(2,:),'linewidth',2); %ylim([0 .25]);
    legend('Odor 1','Jittered Odor 1');
    title 'PSTH'; xlabel 'Time (s)'
%  -------------------------------------------
 
%% Simulate responses evoked by Odor 2 and jittered version of Odor 2

    % add jitter for the peak position
    % for oct condition
    jitter_oct = generateJitter(peakIndex_oct, neuron_number); % generate jitter level
    [neuralActivity_odor2, peakIndex_2] = addJitter(spikeActivity_2,jitter_oct,neuron_number); % add jitter to neural activity
    peakIndex_2 = peakIndex_2/sampleRate; % convert time into second
    % for oct(hex) condition
    jitter_octhex = generateJitter(peakIndex_octhex, neuron_number); % generate jitter level
    [neuralActivity_odor2_jitter, peakIndex_2_jitter] = addJitter(spikeActivity_2,jitter_octhex,neuron_number); % add jitter to neural activity
    peakIndex_2_jitter = peakIndex_2_jitter/sampleRate; % convert time into second
    
    % -----  figure plot -------------------------------------------------
    color = color_map(1*size(color_map,1)/2:-1:1,:);% color map
    % plot Fig. 6b
    figure(3);
    subplot(2,2,1);imagesc(timex ,1:100, neuralActivity_odor2); colormap(color); caxis([0 25]);
    title 'Odor 2'; xlabel 'Time (s)'; ylabel 'Cell index'
    subplot(2,2,2);imagesc(timex ,1:100, neuralActivity_odor2_jitter); colormap(color); caxis([0 25]); 
    title 'Jittered Odor 2'; xlabel 'Time (s)'; ylabel 'Cell index'
    distribution_2 = hist(peakIndex_2,tt); 
    distribution_2_jitter = hist(peakIndex_2_jitter,tt);
    distribution_2 = distribution_2./activationSet_2; 
    distribution_2_jitter = distribution_2_jitter./activationSet_2;
    subplot(2,2,3);
    plot(tt,distribution_2,'b-o','MarkerFaceColor',[0 1 1]); axis([0 12000/sampleRate 0 .6]);
    title 'Peak time distribution'; xlabel 'Time to peak'; ylabel 'Probability'
    subplot(2,2,4);
    plot(tt,distribution_2_jitter,'b-o','MarkerFaceColor',[0 0 .3]);axis([0 12000/sampleRate 0 .6]);
    title 'Peak time distribution'; xlabel 'Time to peak'; ylabel 'Probability'
 
    % plot Fig.6c right panel
    figure(4);
    fill([timex,timex(end:-1:1)],[mean(neuralActivity_odor2)-std(neuralActivity_odor2)/sqrt(100),...
        mean(neuralActivity_odor2(:,duration:-1:1))+std(neuralActivity_odor2(:,duration:-1:1))/sqrt(100)],color_Shade(1,:),'lineStyle','none'); hold on;
    fill([timex,timex(end:-1:1)],[mean(neuralActivity_odor2_jitter)-std(neuralActivity_odor2_jitter)/sqrt(100),...
        mean(neuralActivity_odor2_jitter(:,duration:-1:1))+std(neuralActivity_odor2_jitter(:,duration:-1:1))/sqrt(100)],color_Shade(2,:),'lineStyle','none');
    plot(timex,mean(neuralActivity_odor2),'color',color_PSTH(1,:),'linewidth',2);
    plot(timex,mean(neuralActivity_odor2_jitter),'color',color_PSTH(2,:),'linewidth',2); %ylim([0 .25]);
    legend('Odor 2','Jittered Odor 2');
    title 'PSTH'; xlabel 'Time (s)'
%   ----------------------------------------------------------------------------------------------  

%% Trajectory plot 
    [coeff,score, latent] = pca([neuralActivity_odor1,neuralActivity_odor2]'); % PCA analysis
    variance = latent./sum(latent)*100; % variance captured by components
    % plot Fig.6d
    figure(5); hold on; grid on;
    plot3(score(1:duration,1),score(1:duration,2),score(1:duration,3),'color',[0.8, 0, 0],'LineWidth',3);
    plot3(score(duration+1:2*duration,1),score(duration+1:2*duration,2),score(duration+1:2*duration,3),'color',[0 0 0.8],'LineWidth',3);
    projy3 = neuralActivity_odor1_jitter'*coeff(:,1:3);
    projs3 = neuralActivity_odor2_jitter'*coeff(:,1:3);
    plot3(projy3(:,1),projy3(:,2),projy3(:,3),'color',[1 0 .8],'LineWidth',3);
    plot3(projs3(:,1),projs3(:,2),projs3(:,3),'color',[0 .8 1],'LineWidth',3);
    legend('Odor 1','Odor 2','Jittered Odor 1','Jittered Odor 2')
    xlabel (['PC 1 (',num2str(variance(1)),'%)']);
    ylabel (['PC 2 (',num2str(variance(2)),'%)']);
    zlabel (['PC 3 (',num2str(variance(3)),'%)']);
    title 'Simulated Odor Trajectory'
end
%   --------------------------------------------------------------------------------------------------

%% generate jitter level
function jitter = generateJitter(peakIndex, neuron_number)
% generate jitter index for different neural responses
% input  :  peakIndex, time-to-peak distributions
%           neuron_number,  number of neurons 
% output :  jitter index for each neuron
    jitter = [];
    clear randSamp
    for i = 1:neuron_number,
        randSamp = randi([1,length(peakIndex)],1); % randomly sample time-to-peak value from distribution
        jitter = [jitter, peakIndex(randSamp)];
    end
end
%   --------------------------------------------------------------------------------------------------

%% add jitter to neural activity
function [neuralActivity, peakIndex] = addJitter(spikeActivity,jitter,neuron_number)
    neuralActivity = []; % neural activity with peak index constrained by real data
    for ii = 1:neuron_number
        neuralActivity (ii,:) = [zeros(1,jitter(ii)),spikeActivity(ii,1:end-jitter(ii))];
    end
    [~,peakIndex] = max(neuralActivity,[],2); % peak index
    peakIndex(peakIndex==1)=[]; % exclude non-responsive cell
end
%   --------------------------------------------------------------------------------------------------

%% red blue colormap
function clrmp = redblue(n) 
    rch = [(0:n-1)/n,  ones(1,n)];% red channel
    gch = [(0:n-1)/n, (n-1:-1:0)/n]; % green channel
    bch = [ones(1,n),(n-1:-1:0)/n]; % blue channel
    clrmp = [rch(:), gch(:), bch(:)]; % color map
end
%   --------------------------------------------------------------------------------------------------