Problem using ODE45 - current not recognized

Question

Samuele Bolotta 2021-3-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/764801-problem-using-ode45-current-not-recognized

评论： Samuele Bolotta 2021-3-6

AandB.zip

I am using this function:

%Constants set
Cm=0.01; % Membrane Capcitance
dt=0.01; % Time Step ms
t=0:dt:25; %Time Array ms
I = zeros (145,1); %Applied constant
I(40:50,1) = 0.1; %Applied constant
V=-60; % Initial Membrane voltage
m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
y0=[V;n;m;h];
tspan = [0,max(t)];
%Matlab's ode45 function
[time,V] = ode45(@ODEMAT,tspan,y0);
OD=V(:,1);
ODn=V(:,2);
ODm=V(:,3);
ODh=V(:,4);
%% Plots
figure
subplot(3,1,1)
plot(time,OD);
legend('ODE45 solver');
xlabel('Time (ms)');
ylabel('Voltage (mV)');
title('Voltage Change for Hodgkin-Huxley Model');
subplot(3,1,2)
plot(time,I);
legend('Current injected')
xlabel('Time (ms)')
ylabel('Ampere')
title('Current')
subplot(3,1,3)
plot(time,[ODn,ODm,ODh]);
legend('ODn','ODm','ODh');
xlabel('Time (ms)')
ylabel('Value')
title('Gating variables')

With ODEMAT being:

function [dydt,I] = ODEMAT(t,y)
    
    %Constants
    ENa=55.17; % mv Na reversal potential
    EK=-72.14; % mv K reversal potential
    El=-49.42; % mv Leakage reversal potential
    gbarNa=1.2; % mS/cm^2 Na conductance
    gbarK=0.36; % mS/cm^2 K conductance
    gbarl=0.003; % mS/cm^2 Leakage conductance
    
    I = zeros (145,1); %Applied constant
    I(40:50,1) = 0.1; %Applied constant
    Cm = 0.01; % Capacitance
    
    % Values set to equal input values
    V = y(1);
    n = y(2);
    m = y(3);
    h = y(4);
    
    gNa = gbarNa*m^3*h;
    gK = gbarK*n^4;
    gL = gbarl;
    
    INa=gNa*(V-ENa);
    IK=gK*(V-EK);
    Il=gL*(V-El);
    
    for curr = 1:length(I)
        dydt = [((1/Cm)*(I(curr)-(INa+IK+Il))); % This is voltage
            alphan(V)*(1-n)-betan(V)*n; %This is gating conductance A
            alpham(V)*(1-m)-betam(V)*m; %This is gating conductance B
            alphah(V)*(1-h)-betah(V)*h]; %This is gating conductance C
    end
end

So what I'm trying to do here is see what happens to voltage and gating conductances when current I switched from 0 to 0.1 (at timestep 40 until 50). However, it looks like the function doesn't really take 0.1 as input, and keeps using 0 instead. Where am I getting it wrong? I also attach the six alpha and beta functions.

Thanks,

Samuele