%%%% script to illustrate exactly what happens during bimodal "jumps"
%%%% induced by large NMDA conductance fluctuations
%%%% Mark Humphries 23/9/2009

clear all

% found values
load fit_model_results_NEWtuning

rand('state',1); randn('state',1);
% -------------------------------------------------------------------------
% Input parameters
% spike-train parameters: start from 500 events/s
% step to range of spike events/s
N_nmda = 84; alpha_nmda = 0;
N_ampa = 84; alpha_ampa = 0;
N_gaba = 84; alpha_gaba = 0; 

%%% the parameter values for the example histogram/spike-trains (though
%%% this spike-train will be different of course]
mNMDA = 100;
r_nmda = 4; 
r_ampa = 4;
r_gaba = 4;

%% or on a huge scale
% mNMDA = 1;
% N_nmda = 10000;
% r_nmda = 1;

% dopamine levels
D1 = 0;
D2 = 0;

% -------------------------------------------------------------------------
% all PSP parameters in saved file
Egaba = -60; 
Enmda = 0;
Eampa = 0;

% these should stay in the same ratio
PSPampa = Xsyn; %% loaded from file...
PSPnmda = PSPampa / ampa_nmda; PSPgaba = PSPampa ./ ampa_gaba;

% MS neuron parameters in saved file
k = izipars(1); a = izipars(2); b = izipars(3); c = izipars(4); vr = izipars(5); vpeak = izipars(6);

% found MS parameters: X = [C,vt,d]
C = X(1); vt =X(2); d = X(3);

% extra DA model parameters in saved file
KIR = XD1(1);    % KIR modifier 
LCA = XD1(2);    % LCA modifier
vrD1 = vr*(1+D1*KIR);
dD1 = d*(1-D1*LCA);

% D2 - intrinsic
alpha = XD2;
kD2 = k*(1-alpha*D2);

% synaptic
cD1 = Xd1all;
cD2 = Xd2all;

% simulation parameters
T = 1000; % max duration of simulation (milliseconds)
dt = 0.1; % time step 

% init simulation 
t = 0:dt:T;
n = length(t); % number of time points
SynExp_ampa = exp(-dt / ts_ampa);
SynExp_nmda = exp(-dt / ts_nmda);
SynExp_gaba = exp(-dt / ts_gaba);


% scale NMDA conductance
PSPnmda = mNMDA * PSPampa / ampa_nmda;

Ggaba = zeros(1,n); Gampa = zeros(1,n);Gnmda = zeros(1,n);
Igaba = zeros(1,n); Iampa = zeros(1,n); Inmda = zeros(1,n);
BD1all_nmda = zeros(1,n);
vD1all = vr*ones(1,n); uD1all=0*v;

% generate the spike trains
Sampa = spkgen([0:dt:T], N_ampa, r_ampa, alpha_ampa);
Snmda = spkgen([0:dt:T], N_nmda, r_nmda, alpha_nmda);
Sgaba = spkgen([0:dt:T], N_gaba, r_gaba, alpha_gaba);       
S = sum(Sampa + Snmda + Sgaba);   % total spike-events


% do simulation
for i = 1:n-1
    % NMDA gate 
    BD1all_nmda(i+1)  = 1 ./ (1 + (Mg/3.57) * exp(-vD1all(i)*0.062));    % from Moyer et al 
    
    if Snmda(i) >= 1
        i
    end
    % delete events
%     if i == 9709 % 9706
%         Sampa(i) = 0;
%     end
%     if i >= 9650
%         Snmda(i) = 0;
%     end
% 
    Gampa(i+1) = Gampa(i) + (PSPampa .* Sampa(i)./ts_ampa);
    Gampa(i+1) = Gampa(i+1) * SynExp_ampa;

    Gnmda(i+1) = Gnmda(i) + (PSPnmda .* Snmda(i)./ts_nmda);
    Gnmda(i+1) = Gnmda(i+1) * SynExp_nmda;

    Ggaba(i+1) = Ggaba(i) + (PSPgaba .* Sgaba(i)./ ts_gaba); % add the MS PSPs
    Ggaba(i+1) = Ggaba(i+1) * SynExp_gaba;

    % D1 intrinsic + synaptic: 
    
    Iampa(i+1) = (Gampa(i+1) .* (Eampa - vD1all(i)));
    Inmda(i+1) = (1+cD1*D1)*BD1all_nmda(i+1)*(Gnmda(i+1) .* (Enmda - vD1all(i)));
    Igaba(i+1) = (Ggaba(i+1) .* (Egaba - vD1all(i)));
    
    % standard model
    vD1all(i+1) = vD1all(i) + dt*(k*(vD1all(i)-vrD1)*(vD1all(i)-vt)-uD1all(i) + ...
        Iampa(i+1) + Inmda(i+1) + Igaba(i+1))/C;
    uD1all(i+1) = uD1all(i) + dt*a*(b*(vD1all(i)-vrD1)-uD1all(i));
    % spikes?   
    if vD1all(i+1)>=vpeak
        vD1all(i)=vpeak; vD1all(i+1)=c; 
        uD1all(i+1)=uD1all(i+1)+dD1;
    end

end
        

%--------------------------------------------------------------------------
% plot voltage trace 
figure(1); clf
plot(t,vD1all)

spkTs = t(vD1all == vpeak);
%%% pick a first spike to look at... e.g. 975, 1187
tspk = find(spkTs <= 975,1,'last');
ix = spkTs(tspk) / dt;

ix = 9741;

%%% take 100 ms prior to this of v, and all S
seg = 15 / dt; 
vseg = vD1all(ix-seg:ix+1);
Bseg = BD1all_nmda(ix-seg:ix);
Sgaba_seg = Sgaba(ix-seg:ix); gabaix = find(Sgaba_seg >= 1);
Snmda_seg = Snmda(ix-seg:ix); nmdaix = find(Snmda_seg >= 1);
Sampa_seg = Sampa(ix-seg:ix); ampaix = find(Sampa_seg >= 1);
Igaba_seg = Igaba(ix-seg:ix); 
Inmda_seg = Inmda(ix-seg:ix);
Iampa_seg = Iampa(ix-seg:ix);

figure(2); clf
subplot(411), plot(vseg)
subplot(412), plot(Inmda_seg), hold on, plot(nmdaix,ones(numel(nmdaix),1)*300,'r.')
subplot(413), plot(Bseg)

subplot(414), plot(Iampa_seg+Igaba_seg), hold on, plot(ampaix,ones(numel(ampaix),1),'r.'); 
plot(gabaix,ones(numel(gabaix),1),'b.');

% save all 
fname = ['Detailed_single_trial_bimodality_test_D1_' num2str(D1) '.mat'];
save(fname)