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:
Hope this helps!
