System of ODEs with interdependent variables

I have to write a function that solves a differential equation which relies on three parameters. Each parameter is also described by a differential equation that calls the result of the first equation. I am unsure how to nest these as they all require each other as an input. I have initial values for each of them, but don't know where to go from there. I had listed all of the equations and called the ode45 function for each individually, which obviously didn't work as they rely on each other's outputs.

I've attached the equations below. The "main" ODE solves for V in terms of m, n, and h. Additionally, m, n, and h all vary with V based on unique exponential functions which I have already defined. All other variables are constants, I apologize for the clutter. How can I combine these?

dVdt = gna.*m^3.*h.*(Vna-V)+gk.*n^4.*(Vk-V)+gl.*(Vl-V)+I; %solve for V based on m, n, h
dndt = HH_an(V).*(1-n)-HH_bn(V).*n; %solve for n based on V
dhdt = HH_hm(V).*(1-h)-HH_bh(V).*h; %solve for h based on V
dmdt = HH_am(V).*(1-m)-HH_bm(V).*m; %solve for m based on V

Answer 1

[James Tursa]

I would advise writing a derivative function for this. E.g.,

dydt = myderivative(t,y,gna,gk,gl,Vna,Vk,Vl,I)
V = y(1);
n = y(2);
h = y(3);
m = y(4);
dVdt = gna.*m^3.*h.*(Vna-V)+gk.*n^4.*(Vk-V)+gl.*(Vl-V)+I; %solve for V based on m, n, h
dndt = HH_an(V).*(1-n)-HH_bn(V).*n; %solve for n based on V
dhdt = HH_hm(V).*(1-h)-HH_bh(V).*h; %solve for h based on V
dmdt = HH_am(V).*(1-m)-HH_bm(V).*m; %solve for m based on V
dydt = [dVdt;dndt;dhdt;dmdt];
end

Then use your own solver or ode45( ).

[Anna Jacobsen]

Thanks so much for the advice! I set my function up this way, but the ODE solver (ode45) is not accepting dydt as a valid function for double inputs. Any ideas as to why? My initial conditions for ode45 were [V0 m0 h0 n0].

Edit: Figured it out. Need to call dYdt outside of the function where it is defined in order for the ODE solver to recognize it. However, I think I am supposed to both generate and solve this system of differential equations in the same function. Is it possible to do so? Or do I have to do it separately?

[James Tursa]

Are any of your variables symbolic? Can you show the code you are using?

[Anna Jacobsen]

function dYdt = HH(t,varargin)
    % V = HH(t, V0, I)
    % HH (Hodgkins-Huxley model) returns membrane voltage V given input of time vector t, initial
    % voltage V0, and input current I.
    % V0 default -60mV. I default 0nA. 
    %t must be of the form [t_start t_end].
    
% set cell membrane parameters
cm	=    10;	% capacitance, fF/um2
gna	=  1200;	% maximum sodium conductance, pS/um2
gk	=   360;	% maximum potassium conductance, pS/um2
gl 	=     3;	% leak conductance, pS/um2
Vna	=    50;	% sodium reversal potential, mV
Vk	=   -77;	% potassum reversal potential, mV
Vl	= -54.4;	% leak reversal potential mV

% set requirement for time vector
if ~exist('t')
    error('HH requires at least one input parameter, t');
end

% define default initial voltage and current
if numel(varargin) < 1
    V0 = -60;
else 
    V0 = varargin{1};
end

if numel(varargin) < 2
    I = 0;
else 
    I = varargin{2};
end
   
m0 = HH_am(V0)./(HH_am(V0)+HH_bm(V0));
n0 = HH_an(V0)./(HH_an(V0)+HH_bn(V0));
h0 = HH_ah(V0)./(HH_ah(V0)+HH_bh(V0));

% set ODE solver options    
odeset('AbsTol', 10^-6, 'RelTol', 10^-6, 'MaxStep', 0.1);    

%set initial values of parameters based on starting membrane voltage
V = V0;
m = m0;
h = h0;
n = n0;
Y0 = [V0 m0 h0 n0];

% differential equation describing m parameter over time
    % calls specific alpha & beta functions for m
dmdt = HH_am(V).*(1-m)-HH_bm(V).*m;


% differential equation describing h parameter over time
    % calls specific alpha & beta functions for h
dhdt = HH_ah(V).*(1-h)-HH_bh(V).*h;


% differential equation describing n parameter over time
    % calls specific alpha & beta functions for n
dndt = HH_an(V).*(1-n)-HH_bn(V).*n;


% differential equation describing change in voltage over time
dVdt = (gna.*m^3.*h.*(Vna-V)+gk.*n^4.*(Vk-V)+gl.*(Vl-V)+I)./cm;


% combine all differential equations into one function
dYdt = zeros(1,4);
dYdt = [dVdt; dmdt; dhdt; dndt];
end

[James Tursa]

What are all of the HH_etc( )? Are they functions?

[Anna Jacobsen]

Yes! Sorry for the lack of clarification. They are all functions in terms of V and each parameter (m, n, h) has an "alpha" and a "beta" function.

They all look like this, but with unique coefficients:

function ah = HH_ah(V)

ah = .07 * exp(-(V+65)/20);