When using ODE45, can I specify a variable to assume two different values during the timespan?

Question

Samuele Bolotta 2021-3-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/765356-when-using-ode45-can-i-specify-a-variable-to-assume-two-different-values-during-the-timespan

编辑： Samuele Bolotta 2021-3-9

Functions.zip

I have tried in different ways to see what happens to voltage V and gating conductances m, n and h when, at time step x, current I switched from 0 to 0.1, and then at time step x + n it gets back to 0. However, it looks like regardless of what I do, ODE45 still assumes I is either 0 or 0.1 for the whole time span.

This is the function. In this case the current I is 0.1 for the whole time span. Is there a way to make it 0 everywhere except for 10 ms in the middle, for example?

function ODE_Hodgkin_Huxley (varargin)
    t=0:0.1:25; %Time Array ms    
    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);
    [r,c] = size(time);
     
    I = ones (r,c) ./ 10; %Current
    
    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'
    
end
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
    I = 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);
    
    dydt = [((1/Cm)*(I-(INa+IK+Il))); % Normal case
        alphan(V)*(1-n)-betan(V)*n;
        alpham(V)*(1-m)-betam(V)*m;
        alphah(V)*(1-h)-betah(V)*h];
end

I attach the other functions here.

Thank you!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Steven Lord 2021-3-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/765356-when-using-ode45-can-i-specify-a-variable-to-assume-two-different-values-during-the-timespan#answer_641196

Solve the system with V = 0 up until the time when it should change. If you're not sure of the exact time when it should change, use an events function to determine when it should change. See the ballode example for how to use events functions with the ODE solvers.

Use the final result from the first call to the ODE solver to generate initial conditions for the second call to the ODE solver, which solves the ODE with V = 0.1 over the interval when it has that value.

Use the final result from the second call to the ODE solver to generate initial conditions for the third call to the ODE solver, which solves the ODE with V = 0 over the remaining interval.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Samuele Bolotta 2021-3-9

编辑：Samuele Bolotta 2021-3-9

Thank you!

What are the advantages of using this approach rather than the "ODE with Time-Dependent Terms" example that is mentioned in the ode45 documentation?

Best,

Samuele

请先登录，再进行评论。

Answer 2

Jan 2021-3-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/765356-when-using-ode45-can-i-specify-a-variable-to-assume-two-different-values-during-the-timespan#answer_641206

在 MATLAB Online 中打开

ODE45 is designe to integrate smooth functions only. To change a parameter you have to integrate in chunks:

tSwitch1 = 10.0
tSwtich2 = 10.1;
t0 = 0;
...
I = 0.1;
[time1, V1] = ode45(@(t,y), @ODEMAT(t, y, I), [t0, tSwitch1], y0);
I  = 0.0;
y0 = V1(end, :);
[time2, V2] = ode45(@(t,y), @ODEMAT(t, y, I), [tSwitch1, tSwitch2], y0);
I  = 0.1;
y0 = V2(end, :);
[time3, V3] = ode45(@(t,y), @ODEMAT(t, y, I), [tSwitch1, tEnd], y0);
time = cat(1, time1, time2, time3);
V    = cat(1, V1, V2, V3);
...
 
function [dydt] = ODEMAT(t, y, I)
% Omit the definition of I inside this function.
...
end