%........................Section 1: usage..................................

% See the README.txt file located in the same folder

%..........................................................................

clear all
close all

%........................Section 2: Parameters.............................

% Define the ion channel conductances/area
gASIC= 10; % mS/cm^2
gK = 45;   % mS/cm^2
gNa = 135; % mS/cm^2  
gL = 0.50; % mS/cm^2 

% Define membrane and channels parameters 
Easic=50;       % Reversal potential for ASIC               
Cm = 1;         % uF/cm^2 ; membrane capacitance
% Parameters for the Hodgkin-Huxley (HH) neuronal model
EK = -12-60;    % mV     ; Reversal potential for the K+ current
ENa = 115-60;   % mV     ; Reversal potential for the Na+ current
EL = 10.613-60; % mV     ; Reversal potential for the leak current
% The offset by 60 mV leads to a resting potential of -60 mV

% Settings for the pH Ramp
% Initial and final pH values
pHi=7.4; 
pHf=6.0;

targetRamp=3000; % ms
SimRamp = targetRamp * 1.1527 - 0.044; % Formula for 10-90 fall time
% Because of the sigmoid-shaped function, a value of e. g. 1153 ms will reproduces the experimental 1000 ms pH ramp 

BeforeRamp=2000; % Time in ms at initial pH before application of the ramp


SimDur= BeforeRamp + SimRamp + 5000; % Simulation duration; extend by 5000 ms after the end of the ramp ramp
% The time after (default 5000 for SimRamp 1153 and lambda 2) must be larger than 
% the ratio SimRamp/lambda
% The parametre lambda is defined section 3

% Set up further simulation parameters
dt = 0.01; % tine step in ms
n_steps= round(SimDur)/dt; % number of steps
% Generate tinebase (time axis)
timebase=dt.*(0:n_steps-1);
% If desired, generate a current stimulation waveform (to stimulate the model electrically) 
IstA=zeros(1,n_steps); % here zero (no electrical stimulation)

%..........................................................................

%...................Section 3: Sigmoid-shaped pH protocol..................
% lambda is a parameter shaping the ramp
lambda = 2;
% Note: The time added after the end of the ramp (default 5000 for SimRamp=1153 and lambda=2) 
% must be larger than the ratio SimRamp/lambda
ParDur=SimRamp/lambda; 
ddProt=zeros(1,n_steps);
% Define second derivative
ddProt(timebase>=BeforeRamp & timebase<BeforeRamp+ParDur)=1;
ddProt(timebase>=BeforeRamp+SimRamp & timebase<BeforeRamp+SimRamp+ParDur)=-1;
% Integrate twice
% The double integration leads to a ramp consisting of piecewise linear and parabolic segments
dProt=cumsum(ddProt);   
Prot=cumsum(dProt);
% Normalize
Prot=Prot-min(Prot);
Prot=Prot/max(Prot);
% Final protocol
pHprotocol = pHi+(pHf-pHi)*Prot;   
%..........................................................................

%........................Section 4: Storing arrays.........................
% Time
t_arr=zeros(1, n_steps);
% Gating variables of the ASIC
a_arr=zeros(1,n_steps);
s_arr=zeros(1,n_steps);
% Gates of the HH neuronal model
m_arr=zeros(1, n_steps);
h_arr=zeros(1, n_steps);
n_arr=zeros(1, n_steps);
% Membrane potential
V_arr=zeros(1, n_steps);
% Currents
I_arr=zeros(1, n_steps); % total current
INa_arr=zeros(1, n_steps); % Na+ current
IK_arr=zeros(1, n_steps); % K+ current
IL_arr=zeros(1, n_steps); % Leak current
Iasic_arr=zeros(1, n_steps); % ASIC current
%External, defined a priori
Istim_arr=zeros(1, n_steps); % electrical stimulation current
pH_arr=zeros(1, n_steps); % pH
%..........................................................................

%.......................Section 5: Initial conditions......................
t = 0; % Time in ms
V = -60; % Membrane potential in mV 
% Gates of the HH neuronal model
m=0.0529; % steady state at at V=-60 mV
h=0.5961; % steady state at at V=-60 mV
n=0.3177; % steady state at at V=-60 mV
% Gating variables of ASIC
[a_inf,s_inf,tau_a,tau_s]=f_HH_ASIC(pHprotocol(1));
a=a_inf;
s=s_inf;
%..........................................................................

%.......................Section 6: Main simulation loop....................
for i = 1:n_steps
    
    % Compute currents of the HH neuronal model in uA/cm^2
    IK = gK.*n.*n.*n.*n.*(V-EK);
    INa = gNa.*m.*m.*m.*h.*(V-ENa);
    IL = gL.*(V-EL);
    
    %Compute Iasic in uA/cm^2
    Iasic=gASIC.*a.*s.*(V-Easic);
    
    % get extrinsically imposed conditions (pH, electrical stimulation current)
    Istim=IstA(i); % get a priori defined Istim
    pH=pHprotocol(i); % get a priori defined pH    
    
    % Compute total current
    Itot = (INa + IK + IL) + Istim + Iasic;
    
    % Compute gating functions for the HH neuronal model (functions of V)
    
    % m gate
    alpham =  (2.5-0.1.*(V+60))./(exp(2.5-0.1.*(V+60))-1); 
    betam =4.*exp(-(V+60)./18);                             
    % h gate
    alphah= 0.007.*exp(-(V+60)./20);                         
    betah = 0.1./(exp(3-0.1.*(V+60))+1);                      
    % n gate
    alphan=(0.1-0.01.*(V+60))./(exp(1-0.1.*(V+60))-1);      
    betan = 0.125.*exp(-(V+60)./80);                        
    
    % Compute gating functions for the HH neuronal model (functions of pH)
    % by calling the function "f_HH_ASIC"
    [a_inf,s_inf,tau_a,tau_s]=f_HH_ASIC(pH);
    
    % Numerical integration of the model variables
    
    % Compute all derivatives
    % derivatives of the HH neuronal model gates
    dndt = alphan.*(1-n)-betan.*n;
    dmdt = alpham.*(1-m)-betam.*m;
    dhdt = alphah.*(1-h)-betah.*h; 
    % derivatives of the ASIC model gates
    dadt=(a_inf-a)/tau_a;
    dsdt=(s_inf-s)/tau_s;
    % derivative of membrane potential
    dVdt = -1./Cm.*Itot;
        
    % Integration using the forward Euler method
    % Time and membrane potential
    t = t + dt;
    V = V + dt.*dVdt;
    % HH neuronal model gates
    n = n + dt.*dndt;
    m = m + dt.*dmdt;
    h = h + dt.*dhdt;
    % ASIC gates
    a = a + dt.*dadt;
    s = s + dt.*dsdt;
        
    % Storage into corresponding arrays
    
    % Time and membrane potential
    t_arr(i)=t;
    V_arr(i)=V;
    % Extrinsic conditions
    Istim_arr(i)=Istim;
    pH_arr(i)=pH;
    % Gating variables of the HH neuronal model
    m_arr(i)=m;
    h_arr(i)=h;
    n_arr(i)=n;
    % Gating variables of the ASIC current
    a_arr(i)=a;
    s_arr(i)=s;    
    % Currents
    INa_arr(i)=INa;
    IK_arr(i)=IK;
    IL_arr(i)=IL;
    Iasic_arr(i)=Iasic;    
    I_arr(i)=Itot; % Total current
        
end

%..........................................................................

%...............Section 7: Downsample the output...........................

% Downsampling factor: reduces the nunmber of data points by this factor. 
% This is appropriate or even necessary to decrease the memory/requirements 
% and the figure file size

ds_factor=20; % for a value of 20 and dt=0.01 ms, one point is shown efery 0.2 ms
n_downsampled_points=floor(n_steps/ds_factor);
dsi=1+ds_factor.*(0:n_downsampled_points-1); % indices of downsampled values

% Downsample
t_dsi=t_arr(dsi);
Istim_dsi=Istim_arr(dsi);
pH_dsi=pH_arr(dsi);
V_dsi=V_arr(dsi);
m_dsi=m_arr(dsi);
h_dsi=h_arr(dsi);
n_dsi=n_arr(dsi);
a_dsi=a_arr(dsi);
s_dsi=s_arr(dsi);
INa_dsi=INa_arr(dsi);
IK_dsi=IK_arr(dsi);
IL_dsi=IL_arr(dsi);
Iasic_dsi=Iasic_arr(dsi);

%..........................................................................

%..............Section 8: Graphical output in one figure ..................

% Get screen resolution and setup the figure accordingly
Pix_SS = get(0,'screensize');
xscr=Pix_SS(3);
yscr=Pix_SS(4);
figure;
set(gcf,'Position',[xscr*0.45 yscr*0.03 xscr*0.5 yscr*0.9])

% Subplot for the pH ramp
ax(1)=subplot(511);
plot(t_dsi, pH_dsi,'k','LineWidth',2)
title('pH');
ylabel('pH')
ylim([5 7.5]);
xlim([0 SimDur])

% Subplot for membrane potential
ax(2)=subplot(512);
plot(t_dsi, V_dsi,'k','LineWidth',2) 
hold on
plot([0 SimDur],[EK EK],':','Color',[0 0.6 1]) 
plot([0 SimDur],[ENa ENa],':','Color',[1 0.5 0]) 
plot([0 SimDur],[0 0],':k')
hold off
title('Membrane potential');
ylabel('mV')
ylim([EK-5 ENa+5]);
xlim([0 SimDur])

% Subplot for the gating variables of the HH neuronal model (gate open probabilities)
ax(3)=subplot(513);
plot(t_dsi, m_dsi,'r','LineWidth',2)
hold on
plot(t_dsi, h_dsi,'g','LineWidth',2)
plot(t_dsi, n_dsi,'Color',[0 0.6 1],'LineWidth',2)
plot([0 SimDur],[0 0],':k')
plot([0 SimDur],[1 1],':k')
hold off
title('Gates')
ylabel('Open probability')
ylim([-0.05 1.05]);
legend('m','h','n')
xlim([0 SimDur])

% Subplot for the gating variables of the ASIC model (gate open probabilities)
ax(4)=subplot(514);
plot(t_dsi, a_dsi,'r','LineWidth',2)
hold on
plot(t_dsi, s_dsi,'Color',[0 0.6 1],'LineWidth',2)
plot(t_dsi, a_dsi.*s_dsi,'g','LineWidth',2)
plot([0 SimDur],[0 0],':k')
plot([0 SimDur],[1 1],':k')
hold off
title('ASIC probabilities')
ylabel('Probability')
ylim([-0.05 1.05]);
legend('a','s','a*s')
xlim([0 SimDur])

% Subplot for the ion currents
ax(5)=subplot(515);
plot(t_dsi, INa_dsi,'Color',[1 0.5 0],'LineWidth',2)
hold on
plot(t_dsi, IK_dsi,'Color',[0 0.6 1],'LineWidth',2)
plot(t_dsi, IL_dsi,'b','LineWidth',2)
plot(t_dsi, Iasic_dsi,'m','LineWidth',2)
plot([0 SimDur],[0 0],':k')
hold off
title('Ion currents')
ylabel('\muA/cm^2')
xlabel('Time (ms)')
legend('I_{Na}','I_K','I_L','I_{ASIC}')
xlim([0 SimDur])

% Link all the x-axes (time)
linkaxes(ax,'x');
xlim([0 SimDur])

%..........................................................................
