% 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.file:///home/seva/work/ReConv/data/Invivo/weighavg.m

%  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.


seed=123456789;

% Summing areas of dendrites,soma, hillock and 3 compartments of axon 
load e4mfnorm.asc
load e3mfnorm.asc
load e2mfnorm.asc
load e1mfnorm.asc

% No Nmda situ
load e4mfnonmda.asc
load e3mfnonmda.asc
load e2mfnonmda.asc
load e1mfnonmda.asc

% Summing areas of dendrites,soma, hillock and 3 compartments of axon
% -Inhibitory cases
load e4i4.asc
load e3i3.asc
load e2i2.asc
load e1i1.asc

% initialize weights - comparitive index
w1=0.2150;% w1 percent of cells fire 1 synapse
w2=0.3310;% w2 perc of cells fire 2 synapses
w3=0.2535;% w3 perc of cells fire 3 synapses
w4=0.2004;% w4 perc of cells fire 4 synapses



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

nX1mf0in=deljitter(e1mfnorm(:,[2]),300);
nX2mf0in=deljitter(e2mfnorm(:,[2]),300);
nX3mf0if=deljitter(e3mfnorm(:,[2]),1);
nX4mf0in=deljitter(e4mfnorm(:,[2]),1);
Xnew=[];

nei1=deljitter(e1i1(:,[2]),300);
nei2=deljitter(e2i2(:,[2]),300);
nei3=deljitter(e3i3(:,[2]),1);
nei4=deljitter(e4i4(:,[2]),1);

% Plot cases
Xnew(:,[1])=deljitter(e4mfnorm(:,[1]),2);
Xnew(:,[2])=((w1*nX1mf0in)+(w2*nX2mf0in)+(w3*nX3mf0if)+(w4*nX4mf0in));
figure
plot(Xnew(:,[1]),Xnew(:,[2]));


Vnew(:,[1])=deljitter(e4i4(:,[1]),2);
Vnew(:,[2])=((w1*nei1)+(w2*nei2)+(w3*nei3)+(w4*nei4));
hold on;


nN1mf0in=deljitter(e1mfnonmda(:,[2]),300);
nN2mf0in=deljitter(e2mfnonmda(:,[2]),300);
nN3mf0if=deljitter(e3mfnonmda(:,[2]),1);
nN4mf0in=deljitter(e4mfnonmda(:,[2]),1);


Nmdanew(:,[1])=deljitter(e4mfnonmda(:,[1]),2);
Nmdanew(:,[2])=((w1*nN1mf0in)+(w2*nN2mf0in)+(w3*nN3mf0if)+(w4*nN4mf0in)); %
title('Weighted measure(ctrl,nonmda,inhib)-SomaHillDend regions');



% new Vectors for fields
noX(:,[1]) = Xnew(:,[1]);
noX(:,[2]) = Xnew(:,[2]);
noV(:,[2]) = Vnew(:,[2])
noV(:,[1]) = Vnew(:,[1]);
noNM(:,[2]) = Nmdanew(:,[2]); 
noNM(:,[1]) = Nmdanew(:,[1]); 


% Going in for the field calculations
ytst = zeros(size(Xnew));
Ynew(:,1) = deljitter(ytst(:,1),1);
Ynew(:,2) = deljitter(ytst(:,2),1);
savscram=[];

for i=1:700
    scram=round((normrnd(300,81)));%(mod(k1,105))*
    savscram=[scram;savscram];
    y_i=deljitter(noX(:,[2]),scram);
    Ynew(:,[2])=Ynew(:,[2])+y_i;
end
Ynew(:,[1])=deljitter(Xnew(:,[1]),1,0);

seed=123456789;
% Going in for the non-NMDA field calculations
ztst = zeros(size(Xnew));
savscram1=[];
Znew(:,1) = deljitter(ztst(:,1),1);
Znew(:,2) = deljitter(ztst(:,2),1);
for i=1:700
    scram1=round((normrnd(300,81))); %
    savscram1=[scram1;savscram1];
    z_i=deljitter(noNM(:,[2]),scram1);
    Znew(:,[2])=Znew(:,[2])+z_i;
end
Znew(:,[1])=deljitter(Xnew(:,[1]),1,0);

seed=1234556789;
% Going in for the inhibitory included field calculations
vtst = zeros(size(Vnew));
savscram2=[];
Unew(:,1) = deljitter(vtst(:,1),1);
Unew(:,2) = deljitter(vtst(:,2),1);
for i=1:700
    scram2=round((normrnd(300,81))); %
    savscram2=[scram2;savscram2];
%    disp('scram2 value is');
%    disp(i);
%    disp(scram2);
    u_i=deljitter(noV(:,[2]),scram2);
    Unew(:,[2])=Unew(:,[2])+u_i;
end
Unew(:,[1])=deljitter(noV(:,[1]),1,0);

seed=123456789;

%for new field related x-axis
disp=1:(1/40):((length(Ynew))/40);
disp1=disp([1:4496]);
dispvar=(disp1)/4;

%plot fields
figure
plot(Ynew([5:4500],[2]));
plot(dispvar,Ynew([5:4500],[2]),'b');
hold on;
plot(dispvar,Znew([5:4500],[2]),'r');
hold on;
plot(dispvar,Unew([5:4500],[2]),'g');
legend('Ctrl- Without inhibition', 'no NMDA', 'with inhibition');
title('Extracellular field reconstruction-ctrl(B), Inhib-on(G) & no-nmda(R)');
hold off

% Plot delay jitters - Uncomment lines for plotting spatial and temporal jitters 
%figure;
%indexscram=1:700;
%plot(savscram,indexscram,'b*'); hold on;
%plot(savscram1,indexscram,'r*');file:///home/seva/work/ReConv/data/Invitro/weighsum.m

%plot(savscram2,indexscram,'g+');
%title('Jitter value distributions no-Inhib,ctrl,no-nmda cases');
%hold off;

%figure;
%[f1 x1]=ksdensity(savscram,'kernel','normal');
%[f2 x2]=ksdensity(savscram1,'kernel','normal');
%[f3 x3]=ksdensity(savscram2,'kernel','normal');
%plot(x1,f1,'b');hold on;
%plot(x2,f2,'r');
%figure; hold on;
%plot(x3,f3,'g');
%hold off;