% ReConv algorithm for reconstructing evoked LFP in cerebellar granular layer
% Plotting the in vitro evoked extracellular field potentials(EP) 
% EP is generated by delaying and convolving many possible responses 
% composed of possible combinations out of a single cell
% Uses function: deljitter - try 'help deljitter' for more info

%  Uses traces from multicompartmental GrC model (see http://senselab.med.yale.edu/ModelDb/showmodel.asp?model=116835)
%  Last updated 11-June-2011
%  Model developer: Shyam Diwakar M.
%  Developed at Amrita School of Biotechnology (India) and at Prof. Egidio D'Angelo's Lab at Univ of Pavia (Italy)
%  Amrita School of Biotechnology, Amritapuri
%  Clappana P.O., Kollam, 690 525, Kerala, India.
%  http://research.amrita.edu/compneuro
%  Email:shyam@amrita.edu

%  Model published as [Diwakar et al., 2011, manuscript accepted, PLoS ONE]
% Shyam Diwakar, Paola Lombardo, Sergio Solinas, Giovanni Naldi, Egidio D'Angelo. "Local field potential modeling predicts dense activation in cerebellar granule cells clusters under LTP and LTD control", PLoS ONE, 2011.


% Summing areas of dendrites,soma, hillock and 3 compartments of axon 
load t4mfnorm.asc
load t3mfnorm.asc
load t2mfnorm.asc
load t1mfnorm.asc

load c4mfnorm.asc
load c3mfnorm.asc
load c2mfnorm.asc
load c1mfnorm.asc



seed=123456789;

% initialize weights - comparitive index
w1=0.15; % w1 percent of cells fire 1 synapse
w2=0.35; % w2 perc of cells fire 2 synapses
w3=0.35; % w3 perc of cells fire 3 synapses
w4=0.15; % w4 perc of cells fire 4 synapses

% initialize weights - for the fields
wn1=w1;
wn2=w2;
wn3=w3;%1;
wn4=w4;%1;

nX1mf0in=deljitter(t1mfnorm(:,2),1300);
nX2mf0in=deljitter(t2mfnorm(:,2),1300);
nX3mf0if=deljitter(t3mfnorm(:,2),1);
nX4mf0in=deljitter(t4mfnorm(:,2),1);
Xnew=[];

nC1mf0in=deljitter(c1mfnorm(:,2),2600);
nC2mf0in=deljitter(c2mfnorm(:,2),2600);
nC3mf0if=deljitter(c3mfnorm(:,2),1);
nC4mf0in=deljitter(c4mfnorm(:,2),1);
X1new=[];

% Plot cases
Xnew(:,1)=deljitter(t4mfnorm(:,1),1);
Xnew(:,2)=((w1*nX1mf0in)+(w2*nX2mf0in)+(w3*nX3mf0if)+(w4*nX4mf0in));
X1new(:,1)=deljitter(c4mfnorm(:,1),1);
X1new(:,2)=((w1*nC1mf0in)+(w2*nC2mf0in)+(w3*nC3mf0if)+(w4*nC4mf0in));
figure
%plot(Xnew(:,1),Xnew(:,2));
%hold on;
%plot(X1new(:,1),X1new(:,2));
%title('Weighted measure-SomaHillDend regions');



% new Vectors for fields
noX(:,1) = Xnew(:,1);
noX(:,2) = Xnew(:,2);
noC(:,1) = X1new(:,1);
noC(:,2) = X1new(:,2);
seed=3e7;
% Going in for the field calculations
ytst = zeros(size(Xnew));
Ynew(:,1) = deljitter(ytst(:,1),1);
Ynew(:,2) = deljitter(ytst(:,2),1);

%Make copies for averaging
Anew(:,1) = deljitter(ytst(:,1),1);
Anew(:,2) = deljitter(ytst(:,2),1);
Bnew(:,1) = deljitter(ytst(:,1),1);
Bnew(:,2) = deljitter(ytst(:,2),1);
Cnew(:,1) = deljitter(ytst(:,1),1);
Cnew(:,2) = deljitter(ytst(:,2),1);
Dnew(:,1) = deljitter(ytst(:,1),1);
Dnew(:,2) = deljitter(ytst(:,2),1);
Enew(:,1) = deljitter(ytst(:,1),1);
Enew(:,2) = deljitter(ytst(:,2),1);
Fnew(:,1) = deljitter(ytst(:,1),1);
Fnew(:,2) = deljitter(ytst(:,2),1);
Gnew(:,1) = deljitter(ytst(:,1),1);
Gnew(:,2) = deljitter(ytst(:,2),1);
Hnew(:,1) = deljitter(ytst(:,1),1);
Hnew(:,2) = deljitter(ytst(:,2),1);
Inew(:,1) = deljitter(ytst(:,1),1);
Inew(:,2) = deljitter(ytst(:,2),1);

savscram=[];
savscrami1=[];

%For the separate T,C
%Tw(:,1) = deljitter(ytst(:,1),340);
%Cw(:,1) = deljitter(ytst(:,1),300);

for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Ynew(:,[2])=Ynew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;	
end
Ynew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Anew(:,[2])=Anew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Anew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Bnew(:,[2])=Bnew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Bnew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Cnew(:,[2])=Cnew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Cnew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Dnew(:,[2])=Dnew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Dnew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Enew(:,[2])=Enew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Enew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Fnew(:,[2])=Fnew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Fnew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Gnew(:,[2])=Gnew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Gnew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Hnew(:,[2])=Hnew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Hnew(:,[1])=deljitter(Xnew(:,[1]),1,0);

clear scram,scrami1,savscram,savscrami1,y_i,y1_i;
for i=1:110
    scram=abs(round((normrnd(1300,340))));
    savscram=[scram;savscram];
    scrami1=abs(round((normrnd(2600,340))));
    savscrami1=[scrami1;savscrami1];
    y_i=deljitter(noX(:,[2]),scram);
    y1_i=deljitter(noC(:,[2]),scrami1);
    Inew(:,[2])=Inew(:,[2])+y_i+y1_i;
%    Tw=Tw+y_i;
%    Cw=Cw+y1_i;
end
Inew(:,[1])=deljitter(Xnew(:,[1]),1,0);


Snew(:,[1])=deljitter(Xnew(:,[1]),1,0);
Snew(:,[2])=((Ynew(:,[2]))+(Anew(:,[2]))+(Bnew(:,[2]))+(Cnew(:,[2]))+(Dnew(:,[2]))+(Enew(:,[2]))+(Fnew(:,[2]))+(Gnew(:,[2]))+(Hnew(:,[2]))+(Inew(:,[2])))/10;


%for new field related x-axis (check this for other lengths of simulations - assumed here tstop=200ms with 40 points
disp0=1:(1/40):((length(Ynew)+340)/40);
disp1=disp0([1:17998]);
dispvar=(disp1)/4;

yindex=0:(length(Ynew)-1);
%plot fields
figure
plot(dispvar,Snew(:,[2]));
title('In vivo evoked LFP reconstruction');

% Plot delay jitters - Uncomment lines for plotting spatial and temporal jitters 
%figure
%subplot(2,1,1)
%plot(dispvar,Tw);
%subplot(2,1,2)
%plot(dispvar,Cw);
 
%figure;
%indexscram1=1:110;
%indexscram2=111:220;
%plot(savscram,indexscram1,'b*'); hold on;
%plot(savscrami1,indexscram2,'g');
%title('Jitter value distributions');
%hold off;

%figure;
%[f1 x1]=ksdensity(savscram,'kernel','normal');
%plot(x1,f1,'b');hold on;