%function [spikes,Vrec]=runhh(Ne,M,T,dt,I_mu,I_dev)
%
%  Simulates the activity in a given network using the HH-type model with
%  four ionic and three synaptic currents (Golomb et al. 2006) and synaptic
%  depression (Golomb & Amitai 1997). Euler-Maruyama method is used for
%  integration.
%
%  Input:
%    Ne    - number of excitatory cells
%    M     - connectivity matrix with synaptic weights
%    T     - simulation time (ms)
%    dt    - simulation time step (ms)
%    I_mu  - constant current injection to each cell (uA/cm2)
%    I_dev - deviation of the Brownian current injection to each cell (uA/cm2)
%    
%  Output:
%    spikes - Nx2 matrix of the spiking data, where first column shows the
%             time of spike and second column the index of the neuron that spiked
%    Vrec   - the membrane potential time course of the first excitatory and
%             first inhibitory neuron (do not use if only excitatory neurons)
%
%  Tuomo Mäki-Marttunen
%  Last modified 8.1.2013

function [spikes,Vrec]=runhh(Ne,M,T,dt,I_mu,I_dev)

if nargin < 1 || isempty(Ne)
    Ne = 1; %number of excitatory cells
end

if nargin < 2 || isempty(M)
    Ni = 0; %number of inhibitory cells
    N = Ne + Ni;
    M = rand(N,N) < 0.5; %the connectivity matrix is picked by random
    M(~~eye(N)) = 0;     %no autapses
else
    N = size(M,1);
    Ni = N - Ne;
end

if nargin < 3 || isempty(T)
    T = 1000; %ms
end

if nargin < 4 || isempty(dt)
    dt = 0.002; %ms
end

if nargin < 5 || isempty(I_mu)
    I_mu = 0; %mS/cm^2
end

if nargin < 6 || isempty(I_dev)
    I_dev = 1; %mS/cm^2
end

spike_thr = -30; %mV; membrane potential peaks that exceed this threshold are considered spikes
N_spikes_max = 500000; %max. number of spikes stored in memory

%The parameters for membrane potential dynamics taken from Golomb et al.
%2006, Journal of Neurophysiology.
C = 1; %microF/cm^2
C_I = 1; %microF/cm^2
g_Na = 35; %mS/cm^2
g_NaP = 0.2;
V_Na = 55; % mV
theta_m = -30; %mV
sigma_m = 9.5; %mV
theta_h = -45; %mV
sigma_h = -7; %mV
theta_th = -40.5; %mV
sigma_th = -6; %mV

theta_p = -47; %mV
sigma_p = 3; %mV

g_Kdr = 6; %mS/cm^2
%Note: in Golomb et al. 2006 g_Kdr is said to be 3mS/cm^2, but in the
%corresponding modelDB simulator it is set 6mS/cm^2. 6mS/cm^2 seems to 
%give the right spike shape.
V_K = -90; %mV
theta_n = -33; %mV
sigma_n = 10; %mV
theta_tn = -27; %mV
sigma_tn = -15; %mV

g_K_slow = 1.8; %mS/cm^2
theta_z = -39; %mV
sigma_z = 5; %mV
tau_z = 75; %ms

g_L = 0.05; %mS/cm^2
V_L = -70; %mV

theta_s = -20; %mV
sigma_s = 2; %mV
k_fP = 1; %(ms)^-1
g_AMPA = 0.08;
g_EI_AMPA = 0.2;
tau_AMPA = 5; %ms
V_Glu = 0; %mV

k_xN = 1; %(ms)^-1
tau_hat_NMDA = 14.3; %ms
k_fN = 1; %(ms)^-1
g_NMDA = 0.07;
g_EI_NMDA = 0.05;
tau_NMDA = 100; %ms
sigma_NMDA = 10; %mV
MgCon = 0.7; %mM (0.7 mM used as a typical Mg-concentration in cultures)
theta_NMDA = 10.5*log(MgCon/38.3); %mV, MgCon in mM
k_t = 1; %(ms)^-1
k_v = 0.001; %(ms)^-1

g_IE_GABAA = 0.05; %mS/cm^2
g_II_GABAA = 0.05; %mS/cm^2; different from (Golomb et al. 2006), where I-I synapses were not considered
V_GABAA = -70; %mV

g_I_Na = 35; %mS/cm^2
V_I_Na = 55; %mV
g_I_Kdr = 9; %mS/cm^2
V_I_Kdr = -90; %mV

g_I_L = 0.1; %mS/cm^2
V_I_L = -65; %mV


k_fA = 1; %(ms)^-1
tau_GABAA = 10; %ms

V = -71*ones(Ne,1); %the excitatory neuron membrane potential is at rest at around -71
h = 0.98+zeros(Ne,1);
n = 0.022+zeros(Ne,1);
z = 0+zeros(Ne,1);
V_I = -64*ones(Ni,1); %the inhibitory neuron membrane potential is at rest at around -64
h_I = 0.78+zeros(Ni,1);
n_I = 0.089+zeros(Ni,1);

%The synaptic gating variables:
s_AMPA = 0+zeros(Ne,1);
s_NMDA = 0+zeros(Ne,1);
x_NMDA = 0+zeros(Ne,1);
s_GABAA = 0+zeros(Ni,1);
%The glutamatergic resource variable (see Golomb & Amitai 1997):
T_Glu = 1+zeros(Ne,1);

spikes = zeros(N_spikes_max,2);
N_spikes = 0;
spike_check_interval = 10; %the interval (#time steps) to check spikes
Vrec_save_interval = 1;    %the interval (#spike check steps) to save the membrane potential data

if nargout > 1 %save membrane potential of first cells only if asked for
    Vrec = zeros(2,ceil(T/dt/spike_check_interval/Vrec_save_interval));
end

Vold = [V; V_I]*ones(1,2); %always save the past two membrane potential values to find a spike

for ti = 1:T/dt
    I_app = I_mu + I_dev/sqrt(dt)*randn(N,1); 
    % NOTE: _divided_ by the sqrt of dt. This way, as dV = ... + dt*I_app/C, we have the coefficient of dW proportional to sqrt(dt).
    % Hence, the variance of the random increment is proportional to dt, as required by the Brownian motion framework.

    I_Na = g_Na./(1+exp(-(V-theta_m)/sigma_m)).^3.*h.*(V-V_Na);
    I_NaP = g_NaP./(1+exp(-(V-theta_p)/sigma_p)).*(V-V_Na);
    I_Kdr = g_Kdr*n.^4.*(V-V_K);
    I_K_slow = g_K_slow*z.*(V-V_K);
    I_L = g_L*(V-V_L);
    
    alpha_m = 0.5*(V_I+35)./(1-exp(-(V_I+35)/10));
    beta_m = 20*exp(-(V_I+60)/18);
    I_I_Na = g_I_Na.*(alpha_m./(alpha_m + beta_m)).^3.*h_I.*(V_I-V_I_Na);
    I_I_Kdr = g_I_Kdr*n_I.^4.*(V_I-V_I_Kdr);
    I_I_L = g_I_L*(V_I-V_I_L);

    I_AMPA = g_AMPA*(V-V_Glu).*(M(1:Ne,1:Ne)'*s_AMPA); % AMPA current to excitatory cells
    I_EI_AMPA = g_EI_AMPA*(V_I-V_Glu).*(M(1:Ne,Ne+1:N)'*s_AMPA); % AMPA current to inhibitory cells

    I_NMDA = g_NMDA./(1+exp(-(V-theta_NMDA)/sigma_NMDA)).*(V-V_Glu).*(M(1:Ne,1:Ne)'*s_NMDA); % NMDA current to excitatory cells
    I_EI_NMDA = g_EI_NMDA./(1+exp(-(V_I-theta_NMDA)/sigma_NMDA)).*(V_I-V_Glu).*(M(1:Ne,Ne+1:N)'*s_NMDA); % NMDA current to inhibitory cells

    I_IE_GABAA = g_IE_GABAA*(V-V_GABAA).*(M(Ne+1:N,1:Ne)'*s_GABAA); % GABA current to excitatory cells
    I_II_GABAA = g_II_GABAA*(V_I-V_GABAA).*(M(Ne+1:N,Ne+1:N)'*s_GABAA); % GABA current to inhibitory cells
    
    dh = (1./(1+exp(-(V-theta_h)/sigma_h)) - h)./(0.1+0.75./(1+exp(-(V-theta_th)/sigma_th)));
    dn = (1./(1+exp(-(V-theta_n)/sigma_n)) - n)./(0.1+0.5./(1+exp(-(V-theta_tn)/sigma_tn)));
    dz = (1./(1+exp(-(V-theta_z)/sigma_z)) - z)/tau_z;
    
    alpha_h =  0.35*exp(-(V_I+58)/20);
    beta_h = 5./(1+ exp(-(V_I + 28)/10));
    dh_I = (alpha_h.*(1 - h_I) - beta_h.*h_I);
    alpha_n =  0.05*(V_I+34)./(1-exp(-(V_I+34)/10));
    beta_n = 0.625*exp(-(V_I + 44)/80);
    dn_I = (alpha_n.*(1-n_I) - beta_n.*n_I);
    
    ds_AMPA = (k_fP./(1+exp(-(V-theta_s)/sigma_s)).*T_Glu.*(1-s_AMPA) - s_AMPA/tau_AMPA);
    dx_NMDA = (k_xN./(1+exp(-(V-theta_s)/sigma_s)).*(1-x_NMDA) - (1-1./(1+exp(-(V-theta_s)/sigma_s))).*x_NMDA/tau_hat_NMDA);
    ds_NMDA = (k_fN.*x_NMDA.*T_Glu.*(1-s_NMDA) - s_NMDA/tau_NMDA);
    ds_GABAA = (k_fA*1./(1+exp(-(V_I-theta_s)/sigma_s)).*(1-s_GABAA)-s_GABAA/tau_GABAA);
    dT_Glu = -k_t./(1+exp(-(V-theta_s)/sigma_s)).*T_Glu + k_v*(1-T_Glu);
    
    dV = 1/C*(-I_Na-I_NaP-I_Kdr-I_K_slow-I_L-I_AMPA-I_NMDA-I_IE_GABAA+I_app(1:Ne));
    dV_I = 1/C_I*(-I_I_Na-I_I_Kdr-I_I_L-I_EI_AMPA-I_EI_NMDA-I_II_GABAA+I_app(Ne+1:N));
    
    % Euler-Maruyama:
    V = V+dt*dV;
    h = h+dt*dh;
    n = n+dt*dn;
    z = z+dt*dz;
    V_I = V_I+dt*dV_I;
    h_I = h_I+dt*dh_I;
    n_I = n_I+dt*dn_I;
    s_AMPA = s_AMPA+dt*ds_AMPA;
    s_NMDA = s_NMDA+dt*ds_NMDA;
    x_NMDA = x_NMDA+dt*dx_NMDA;
    s_GABAA = s_GABAA+dt*ds_GABAA;
    T_Glu = T_Glu + dt*dT_Glu;
    
    %Check the spiking and save the membrane potential values if needed:
    if mod(ti,spike_check_interval) == 0
        %There was a spike iff there was a local maximum of membrane
        %potential and the maximum value is above spike_thr.
        spiking = [V; V_I] > spike_thr & Vold(:,2) < Vold(:,1) & Vold(:,1) >= [V;V_I];
        if any(spiking) && N_spikes < N_spikes_max
            spiking = find(spiking);
            spvec = [(ti-1)*dt*ones(length(spiking),1), spiking];
            spikes(N_spikes+1:min(N_spikes+length(spiking),N_spikes_max),:) = spvec(1:min(length(spiking),N_spikes_max-N_spikes),:);
            N_spikes = N_spikes + length(spiking);
        end      
        if nargout > 1 && mod(ti/spike_check_interval,Vrec_save_interval)==0
            Vrec(:,ti/spike_check_interval) = [V(1); V_I(1)];
        end
	    Vold = [[V;V_I] Vold(:,1)];
    end
    
end

end
