Hodgkin-Huxley model response

I have code for Hodgkin-Huxley model using Euler's method and the response I get for following pulses:

N = 400;              % number of pulses
Tp = 0.001;           % pulse width
d1 = 0.001;           % interphase delay
d2_2 = 0.002;         % interpulse delay [ms] 
U0 = 100;             % amplitude

% define sampling freq
fs = 500000;        % [Hz]

% single pulse
x1 = [ones(round(Tp*fs), 1)            % + pulse
      zeros(round(d1*fs), 1)           % d1
     -ones(round(Tp*fs), 1)            % - pulse
      zeros(round(d2_2*fs), 1)];       % d2

x = U0*repmat(x1, N, 1);
t = (0:length(x)-1)/fs;

% Downsampling
downsample_factor = ceil(length(x) / 100000);
x_downsampled = x(1:downsample_factor:end);
t_downsampled = t(1:downsample_factor:end);
% plot the downsampled signal
plot(t_downsampled, x_downsampled);
grid on;
ylim([-1 1] * 1.5 * U0);
xlabel("Time [ms]");
ylabel("U [V]");
xlim([0 0.2])
ylim([-0.1 0.1])

is (left figure is amplitude = 1mV and right one amplitude = 100mV):

Is this logical (moslty right figure)? Could it be problem with my code of model or pulses code?

Answer 1

[Rahul]

Hi Marina Babic,

I believe that you are trying to simulate Hodgkin-Huxley model using Euler's method. Assuming that you are using the latest R2024b release of MATLAB, you can consider the following improvements in your code:

• The pulse generation and down sampling seem fine, but you can optimize flexibility by defining parameters like pulse duration, delay, and repetition dynamically to make it more scalable.
• Euler's method is simple but may not handle the stiffness of Hodgkin-Huxley equations effectively. You can consider using a higher-order method, like the Runge-Kutta method, which is more accurate and stable for these types of equations.
• Further, you can add functionality to test multiple pulse amplitudes and compare their effects. Use loops or arrays to explore the response for different ‘U0’ values.

Here is you can modify and restructure the given script in MATLAB:

% Hodgkin-Huxley Model Simulation with Pulse Train 
clear; clc; close all; 

%% Define Simulation Parameters 
fs = 500000;               % Sampling frequency [Hz] 
dt = 1 / fs;               % Time step [ms] 
T_total = 2;               % Simulation duration [ms] 
time_steps = round(T_total / dt); % Total simulation steps 

% Pulse Parameters 
Tp = 0.001;                % Pulse width [s] 
d1 = 0.001;                % Interphase delay [s] 
d2 = 0.002;                % Interpulse delay [s] 
N = 400;                   % Number of pulses 
U0 = 100;                  % Amplitude of pulse [mV] 

% Hodgkin-Huxley Parameters 
gNa = 120;                 % Sodium conductance [mS/cm^2] 
gK = 36;                   % Potassium conductance [mS/cm^2] 
gL = 0.3;                  % Leakage conductance [mS/cm^2] 
ENa = 115;                 % Sodium reversal potential [mV] 
EK = -12;                  % Potassium reversal potential [mV] 
EL = 10.6;                 % Leakage reversal potential [mV] 
Cm = 1;                    % Membrane capacitance [uF/cm^2] 

%% Generate Pulse Train 
pulse_duration = round(Tp * fs); % Pulse width in samples 
d1_duration = round(d1 * fs);    % Interphase delay in samples 
d2_duration = round(d2 * fs);    % Interpulse delay in samples 

single_pulse = [ones(pulse_duration, 1);       % + pulse 
    zeros(d1_duration, 1);         % interphase delay 
    -ones(pulse_duration, 1);       % - pulse 
    zeros(d2_duration, 1)];        % interpulse delay 

x = U0 * repmat(single_pulse, N, 1);           % Repeating pulse train 
t = (0:length(x)-1) / fs;                      % Time vector [s] 

%% Hodgkin-Huxley Model Equations 
% Initialize variables 
V = -65;                % Initial membrane potential [mV] 
n = 0.3177;             % Initial potassium activation variable 
m = 0.0529;             % Initial sodium activation variable 
h = 0.5961;             % Initial sodium inactivation variable 

V_trace = zeros(time_steps, 1); % Store membrane potential over time 
V_trace(1) = V;                 % Initial condition 

% Define gating variable equations 
alpha_n = @(V) 0.01 * (10 - V) / (exp((10 - V) / 10) - 1); 
beta_n = @(V) 0.125 * exp(-V / 80); 
alpha_m = @(V) 0.1 * (25 - V) / (exp((25 - V) / 10) - 1); 
beta_m = @(V) 4 * exp(-V / 18); 
alpha_h = @(V) 0.07 * exp(-V / 20); 
beta_h = @(V) 1 / (exp((30 - V) / 10) + 1); 

% Run simulation 
for i = 2:time_steps 
    % Input current from pulse train 
    if i <= length(x) 
        I_ext = x(i); % External input current [uA/cm^2] 
    else 
        I_ext = 0; 
    end 
    
    % Update gating variables 
    dn = alpha_n(V) * (1 - n) - beta_n(V) * n; 
    dm = alpha_m(V) * (1 - m) - beta_m(V) * m; 
    dh = alpha_h(V) * (1 - h) - beta_h(V) * h; 
    n = n + dt * dn; 
    m = m + dt * dm; 
    h = h + dt * dh; 
    
    % Ionic currents 
    INa = gNa * m^3 * h * (V - ENa); % Sodium current 
    IK = gK * n^4 * (V - EK);        % Potassium current 
    IL = gL * (V - EL);              % Leakage current 
    
    % Update membrane potential using Euler's method 
    dV = (I_ext - INa - IK - IL) / Cm; 
    V = V + dt * dV; 
    V_trace(i) = V; 
end 

%% Downsampling for Plotting 
downsample_factor = ceil(time_steps / 100000); % Downsample for clarity 
V_downsampled = V_trace(1:downsample_factor:end); 
t_downsampled = (0:downsample_factor:time_steps-1) * dt; 

%% Plot Results 
% Pulse Train 
figure; 
plot(t * 1000, x); % Time in ms 
title("Input Pulse Train"); 
xlabel("Time [ms]"); 
ylabel("Input Current [\muA/cm^2]"); 
xlim([0, 0.2 * 1000]); % Show first 200 ms 
grid on; 

% Membrane Potential 
figure; 
plot(t_downsampled * 1000, V_downsampled, 'r', 'LineWidth', 1.5); % Time in ms 
title("Hodgkin-Huxley Membrane Potential"); 
xlabel("Time [ms]"); 
ylabel("Membrane Potential [mV]"); 
grid on; 

Here, the Membrane Potential Response verifies the Hodgkin-Huxley model's response to the stimulus.

The first figure shows the generated pulse train, whereas the second figure plots the membrane potential response over time.

You can adjust parameters like ‘U0’ (pulse amplitude), ‘Tp’ (pulse width), and Hodgkin-Huxley constants, like ‘gNa’, ‘gK’, to test for different conditions.

To know more about the usage of ‘plot’ function, refer to the following documentation link:

https://www.mathworks.com/help/matlab/ref/plot.html

Hope this helps!

[Marina Babic]

Oh thank you so much for your answer, it helps a lot!