MATLAB: Solve ODE with a time dependent parameter

differential equationsMATLABode45

Hi,
I am working through Comp Neruo Book in my spare time.
I have been manually coding solutions for the ODE's as it has provided me personally with more flexibility, but I am trying to learn how the numerical ODE suite works.
If I have a time varing parameter (Iapp) that changes when t is divisible by 0.1 seconds, how would that be modelled. I would also like to return the paramater as an array for when it. I don't see how that would would be implemented with the "Event" .
The equations below are the hodgkin-huxley model and y is
y = [V, m, h, n]
though this detail is not relevant.
The Iapp value that is returned makes no sense to me and it has no effect on the calculation anyway…
What is the correct approach to modelling this in the numerical ode suite? 
clear
clc
tspan = 0:0.0001:0.35;
y0 = [0;0;0;0];
Iapp = 0;
[t,y,Iapp] = ode45(@(t,y) HHmodelD(t,y,Iapp),tspan,y0);
   
plot(t,y(:,1))
title("Membrane Potential")
ylabel("Volts")
xlabel("seconds")
function [dydt,Iapp] = HHmodelD(t,y,Iapp)
    
    GL = 30e-9; GnaK = 12e-6; Gk_Max = 3.6e-6; Ena = 45e-3;
    Ek = -82e-3; EL = -60e-3; Cm = 100e-12; 
    
%-------------------------------------------------------

% Where I am modifying the Iapp Value
    if(mod(t,0.1) == 0)    
        Iapp(end+1) = Iapp(end) + 0.22e-9;
    else
        Iapp(end+1) = Iapp(end);
    end
%-------------------------------------------------------
% Note Iapp(end) in dydt1
    dydt1 = (GL/Cm)*(EL - y(1)) +(Gk_Max/Cm)*(y(2)^3) * y(3)*(Ena - y(1)) + (Gk_Max/Cm)*(y(4)^4) *(Ek - y(1)) + Iapp(end)/Cm;
    dydt2 = alphaM(y(1))*(1-y(2))-betaM(y(1))*y(2);
    dydt3 = alphaH(y(1))*(1-y(3))-betaH(y(1))*y(3);
    dydt4 = alphaN(y(1))*(1-y(4))-betaN(y(1))*y(4);
    
    dydt = [dydt1;dydt2;dydt3;dydt4];
        
            
      function aM = alphaM(V)
            aM = ((10^5)*(-V-0.045))/(exp(100*(-V-0.045)) - 1);
      end
      function aH = alphaH(V)
            aH = 70*exp(50*(-V-0.070));
      end
      function aN = alphaN(V)
            aN = ((10^4)*(-V-0.060))/(exp(100*(-V-0.060)) - 1);
      end
      function bM = betaM(V)
            bM = (4*10^3)*exp((-V - 0.070)/0.018);
      end
      function bH = betaH(V)
            bH = (10^3)/(1+exp(100*(-V-0.040)));
      end
      function bN = betaN(V)
            bN = 125*exp((-V-0.070)/0.08); 
      end
      
end

Best Answer

Using discrete values for ‘Iapp’ is likely not appropriate because of the adaptive algorithms used in the MATLAB ODE solvers.
See ODE with Time-Dependent Terms for an approach that will likely work.
Related Question