未能使用插值而不是在带有 ODE45 的块中集成

编程问答 2022-05-17

问题描述

我编写了这段代码来在 30ms 内对电压 V 和门控电导 m、n 和 h 进行积分，并注入从 10 毫秒到 20 毫秒的电流。所以基本上电流 I 是一个向量，它在 10ms 之前具有 0 值，在 20ms 之前具有 1 值和 0 值直到 30ms。代码运行良好，显示了我想看到的行为。

我正在尝试使用插值而不是块集成来实现相同的代码。但是，我没有这样做。

这是与插值的集成，它给出了不正确的结果：

myode = @(T,y0) [((1/Cm)*interp1(t,I,T,'linear','extrap')-(INa+IK+Il)); % normal case
    alpham(y0(1,1))*(1-y0(2,1))-betam(y0(1,1))*y0(2,1);
    alphan(y0(1,1))*(1-y0(3,1))-betan(y0(1,1))*y0(3,1);
    alphah(y0(1,1))*(1-y0(4,1))-betah(y0(1,1))*y0(4,1)];

[time,V] = ode45(myode,t,y0,I);

这是正确的，分块：

[time,V] = ode45(@ODEMAT,[t(chunk),t(chunk+1)],y);   

%%% Calls ODEMAT,in which theres' the integration (extensive code posted below):

dydt = [((1/Cm)*(I(chunk)-(INa+IK+Il))); % normal case
    alpham(V)*(1-m)-betam(V)*m;
    alphan(V)*(1-n)-betan(V)*n;
    alphah(V)*(1-h)-betah(V)*h];

这是分块集成的代码：

function Z2_chunks ()
    %% Initial values
    V=-60; % Initial Membrane voltage
    m1=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
    n1=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
    h1=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
    y0=[V;m1;n1;h1];
    
    t = [1:30];
    I = [zeros(1,10),ones(1,zeros(1,10)];
    
    % Plotting purposes (set I(idx) equal to last value of I)
    idx = numel(t);
    I(idx) = 0.1;
    
    chunks = numel(t) - 1;
    
    for chunk = 1:chunks
        
        if chunk == 1
            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
            y=[V;m;n;h];
        else
            y = V(end,:);  % Final position is initial value for next interval
        end
        [time,y);
        
        if chunk == 1
            def_time = time;
            def_v = V;
        else
            def_time = [def_time; time];
            def_v = [def_v; V];
        end       
    end
    
    OD = def_v(:,1);
    ODm = def_v(:,2);
    ODn = def_v(:,3);
    ODh = def_v(:,4);
    time = def_time;
    
    %% Plots
    %% Voltage
    figure
    subplot(3,1,1)
    plot(time,OD);
    legend('ODE45 solver');
    xlabel('Time (ms)');
    ylabel('Voltage (mV)');
    title('Voltage Change for Hodgkin-Huxley Model');
    
    %% Current
    subplot(3,2)
    stairs(t,I)
    ylim([0 5*max(I)])
    legend('Current injected')
    xlabel('Time (ms)')
    ylabel('Ampere')
    title('Current')
    
    %% Gating variables
    subplot(3,3)
    plot(time,[ODm,ODn,ODh]);
    legend('ODm','ODn','ODh');
    xlabel('Time (ms)')
    ylabel('Value')
    title('Gating variables')
        
    function [dydt] = ODEMAT(t,y)
        
        %% Constants
        ENa=55; % mv Na reversal potential
        EK=-72; % mv K reversal potential
        El=-49; % mv Leakage reversal potential
        
        %% Values of conductances
        gbarl=0.003; % mS/cm^2 Leakage conductance
        gbarNa=1.2; % mS/cm^2 Na conductance
        gbarK=0.36; % mS/cm^2 K conductancence
        Cm = 0.01; % Capacitance
        
        % Values set to equal input values
        V = y(1);
        m = y(2);
        n = 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);
    
        
        dydt = [((1/Cm)*(I(chunk)-(INa+IK+Il))); % normal case
            alpham(V)*(1-m)-betam(V)*m;
            alphan(V)*(1-n)-betan(V)*n;
            alphah(V)*(1-h)-betah(V)*h]; 
    end
end

这是与插值集成的代码。如您所见，情节确实不同：

function Z1_interpol () 
    %% Initial values
    V=-60; % Initial Membrane voltage
    m1=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
    n1=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
    h1=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
    y0=[V;m1;n1;h1];
    
    t = [1:30];
    I = [zeros(1,10)];
    
    % Plotting purposes (set I(idx) equal to last value of I)
    idx = numel(t);
    I(idx) = 0.1;   
    
    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
    y=[V;m;n;h];
    
    %% Constants
    ENa=55; % mv Na reversal potential
    EK=-72; % mv K reversal potential
    El=-49; % mv Leakage reversal potential
    
    %% Values of conductances
    gbarl=0.003; % mS/cm^2 Leakage conductance
    gbarNa=1.2; % mS/cm^2 Na conductance
    gbarK=0.36; % mS/cm^2 K conductancence
    Cm = 0.01; % Capacitance
    
    %% Initial values
    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;m;n;h];
    
    gNa = gbarNa*m^3*h;
    gK = gbarK*n^4;
    gL = gbarl;
    
    INa=gNa*(V-ENa);
    IK=gK*(V-EK);
    Il=gL*(V-El);
    
    myode = @(T,'extrap')-(INa+IK+Il)); % normal case
        alpham(y0(1,1);
        alphan(y0(1,1);
        alphah(y0(1,1)];
    
    [time,I);
    
    OD=V(:,1);
    ODm=V(:,2);
    ODn=V(:,3);
    ODh=V(:,4);
    time = time;
    
    
        %% Plots
    %% Voltage
    figure
    subplot(3,I)
    ylim([0 5*max(I)])
    legend('Current injected')
    xlabel('Time (ms)')
    ylabel('Ampere')
    title('Current')
    
    
    %% Gating variables
    subplot(3,'ODh');
    xlabel('Time (ms)')
    ylabel('Value')
    title('Gating variables')       
    
end

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

interpolationmatlabmatlabodeode45