%gks and gna parameter sweep



% hodgkin huxley single compartment model with parameters described in Goff
% and Goldberg 2019
%this is a good but simplistic starting point for modelling the effects of Nav1.1 and
%KCNQ5 on VIP interneuron function. Future methods will encorporate more
%detailed modelling of the many conductances that account for VIP-IN
%excitability


clear all
close all
% global variables

global Ek
global El
global ENa
global gl
global gkDR
global gks
global gna
global tau_s



El = -70;%leak reversal mV
Ek = -77;%K reversal mV
ENa = 50;%Na reversal mV

tau_s = 300; %ms
gl = 3*10^-5*10^6;
gkDR = 0.005*10^6;


cm = 1000;%nF/cm^2
Vo = -70;%starting membrane potential
A = 0.0001;%area of membrane cm^2 corresponds to cylinder of l=w=25uM


%plot_voltage_dep_of_conductances_and_rates() % use this function to plot
%voltage dependence of rate constants etc.




length_sw= 700;%length of simulation in milliseconds
dur = 600;%current injection duration mseconds
st_pt = 50;%when current begins to be injected mseconds
si = .025; %sampling interval milliseconds
num_sweeps = 3;%number of simulations
delta_I = 0.1;%change in current for each simulation in nano Amps

matrix_exc = injection(length_sw,si,st_pt,dur,num_sweeps,delta_I);%excitation matrix


T = zeros(length_sw/si,num_sweeps); %timescale matrix at the correct sampling rate for plotting purposes only
for i = 1:num_sweeps 
    T(:,i) = linspace(0,length_sw,length_sw/si);
end

g_ks_vals = [0.0001:0.0001: 0.002].*10^6;
g_na_vals = linspace(0.03,0.07,length(g_ks_vals)).*10^6;%controls the bounds of the parameter sweep, set to single values to test the model at a particular point

n = cell(length(g_ks_vals),length(g_na_vals));

matrix_res = n; m = n; h = n; s = n; Im = n;
for j = 1:length(g_ks_vals)
    gks = g_ks_vals(j);%(0.0001 - 0.002)
    
    for i= 1:length(g_ks_vals)
    gna = g_na_vals(i); %conductances in uS/cm^2 (0.35-0.65)
    
    
    [matrix_res{j,i},n{j,i},m{j,i},h{j,i},s{j,i},Im{j,i}] = response(matrix_exc,length_sw,si,num_sweeps,Vo,cm,A);%neuron response matrix, with n,m,h and Im
    end
end

%save matrix_res and matrix_exc in a single .mat file for analysis with analysis_param_sweep.m

%% functions


%model equations as described in the paper
function [alpha] = am(V)
alpha = (0.1*(V+38))/(1 - exp(-(V + 38)/10));
end
function [beta] = bm(V)
beta = 4*exp(-(V+65)/17);
%(-0.062*(V+38))/(1 - exp((V + 38)/5));
end

function [alpha] = ah(V)
alpha = (0.05* exp(-0.05*(V + 55)));
end
function [beta] = bh(V)
beta = (1)/(exp(-0.1*(V+35))+1);
end

function [alpha] = an(V)
alpha = (0.01*(V+55))/(1 - exp(-0.1*(V + 55)));
end
function [beta] = bn(V)
beta = (0.125 * exp(-0.0125*(V + 65)));
end

function [sinf] = s_inf(V)
sinf = 1/(1+exp(0.2*(-35 - V)));%was -53.4
end


function [time_const] = tau(alpha,beta)
time_const = 1/(alpha + beta);
end

function [x_inf] = steady_state(alpha,beta)
x_inf = alpha/(alpha+beta);
end



function [del_x] = del(x_inf,tau,x)
del_x = (x_inf - x)/tau;
end

function [I] = im(V,n,m,h,s)
global Ek
global El
global ENa
global gl
global gkDR
global gks
global gna

I = (gl*(El-V) + gkDR*(n^4)*(Ek-V) + gna*(m^3)*h*(ENa-V) + gks*s*(Ek-V));%/ms
%I = (gl*(V-El) + gk*(n^4)*(V-Ek) + gna*(m^4)*(V-ENa))*0.001;
end

% build excitation curves, ie define the current to be injected

function [matrix_exc] = injection(length,si,st_pt,dur,num_sweeps,delta_I)
matrix_exc = zeros(length/si,num_sweeps);
for i = 1:num_sweeps
    pA = i*delta_I;%this determines the current injected at each sweep, starting at -5delta_I and incrementing by delta_I each sweep
    for j = 1:length/si
        if j >= st_pt/si && j < (st_pt+dur)/si 
            matrix_exc(j,i) = pA;%only inject current during a defined period of 600ms
        else 
            matrix_exc(j,i) = 0;%otherwise no current injected
        end
    end
end

end


%this function is what models the neurons response
function [V,n,m,h,s,Im] = response(matrix_exc,length,dt,num_sweeps,Vo,cm,A)
global tau_s

V = zeros(length/dt,num_sweeps);%same size as excitation matrix
V(1,:) = Vo;%starting condition for each sweep
n = zeros(length/dt,num_sweeps);
m = n; h = n; s = n;

Im = zeros(length/dt,num_sweeps);
    n(1,:) = steady_state(an(Vo),bn(Vo));%starting value for rate constants
    m(1,:) = steady_state(am(Vo),bm(Vo));
    h(1,:) = steady_state(ah(Vo),bh(Vo));
    s(1,:) = s_inf(Vo);
for i = 1:num_sweeps

    
    for j = 1:length/dt-1
        delV = dt*(((matrix_exc(j,i)+0.05*randn(1))/A + im(V(j,i),n(j,i),m(j,i),h(j,i),s(j,i)))/cm);%A is the area per unit for injected current
        deln = dt*del(steady_state(an(V(j,i)),bn(V(j,i))),tau(an(V(j,i)),bn(V(j,i))),n(j,i));
        delm = dt*del(steady_state(am(V(j,i)),bm(V(j,i))),tau(am(V(j,i)),bm(V(j,i))),m(j,i));
        delh = dt*del(steady_state(ah(V(j,i)),bh(V(j,i))),tau(an(V(j,i)),bh(V(j,i))),h(j,i));
        dels = dt*(s_inf(V(j,i)) - s(j,i))/tau_s;
        
        V(j+1,i) = V(j,i) + delV;%eulers method for the differential equation with added noise and inhibitory current
   
        n(j+1,i) = n(j,i) + deln; 
        m(j+1,i) = m(j,i) + delm;
        h(j+1,i) = h(j,i) + delh;
        s(j+1,i) = s(j,i) + dels;
        Im(j,i) = im(V(j,i),n(j,i),m(j,i),h(j,i),s(j,i));
    end
end
end

function [] = plot_voltage_dep_of_conductances_and_rates()
V = linspace(-100,20);
alpha_N = zeros(size(V));
for i = 1 : size(V,2)
    alpha_N(1,i) = an(V(1,i));
end
beta_N = zeros(size(V));
for i = 1 : size(V,2)
    beta_N(1,i) = bn(V(1,i));
end


N_inf = zeros(size(V));
for i = 1 : size(V,2)
    N_inf(1,i) = steady_state(an(V(1,i)),bn(V(1,i)));
end

M_inf = zeros(size(V));
for i = 1 : size(V,2)
    M_inf(1,i) = steady_state(am(V(1,i)),bm(V(1,i)));
end

H_inf = zeros(size(V));
for i = 1 : size(V,2)
    H_inf(1,i) = steady_state(ah(V(1,i)),bh(V(1,i)));
end

tau_n = zeros(size(V));
for i = 1 : size(V,2)
    tau_n(1,i) = tau(an(V(1,i)),bn(V(1,i)));
end

tau_m = zeros(size(V));
for i = 1 : size(V,2)
    tau_m(1,i) = tau(am(V(1,i)),bm(V(1,i)));
end

tau_h = zeros(size(V));
for i = 1 : size(V,2)
    tau_h(1,i) = tau(ah(V(1,i)),bh(V(1,i)));
end

figure
subplot(1,3,1)
plot(V,alpha_N)
hold on
plot(V,beta_N)
title('rate functions of potassium channel')
xlabel('mV')
ylabel('1/ms')
legend('alpha_n','beta_n')
subplot(1,3,2)
plot(V,N_inf)
title('steady state K channel voltage dependence')
xlabel('mV')
ylabel('ninf')
subplot(1,3,3)
plot(V,tau_n)
title('rate constant')
xlabel('mV')
ylabel('tau(ms)')

figure
hold on
plot(V,N_inf)
plot(V,M_inf)
plot(V,H_inf)
title('steady state activation/inactivation')
xlabel('mV')
legend('Ninf','Minf','Hinf')

figure
hold on
plot(V,tau_n)
plot(V,tau_m)
plot(V,tau_h)
title('voltage dependent time constants')
xlabel('mV')
ylabel('tau(ms)')
legend('tau_n','tau_m','tau_h')

end%use this to monitor changes to the model parameters

