Function for ODE45 solver doesn't work

First function .m file:

function dVm=dVm(t,y)
% Parameters are time t, y for Vm, m, n and h, input for Istim)
% and given Constants
gNa = 120; % mS/cm^2 
gK = 36; % mS/cm^2 
gL = 0.3; % mS/cm^2 
ENa = 115; % mV
EK = -12; % mV
EL = 10.6; % mV 
Cm = 1; % uF/cm^2

% System Functions 
dVm(1,1)= @(Istim_sys)(-(1/Cm)*(gNa*((y(2)^3)*y(4)*(y(1)-ENa))+(gK*((y(3)^4)*(y(1)-EK))+(gL*(y(1)-EL))-Istim_sys)));
dVm(2,1) = (0.01*(0.5-y(1)) / (1-exp((9.5-y(1))/10)))*(1-y(2)) - (0.125*exp(-y(1)/ 80))*y(2); % dn/dt where y(2) = n
dVm(3,1) = (0.1*(25-y(1)) / (1-exp((25-y(1))/10)))*(1-y(3)) - (4*exp(-y(1)/ 18))*y(3); % dm/dt where y(3) = m
dVm(4,1) = (0.07*exp(-y(1) / 20))*(1-y(4)) - (1 / (1+exp((30-y(1))/10)))*y(4); %  dh/dt where y(4) = h
end

2nd function .m file:

function Istim_sys
Ts1 = 0.5;
Ts2 = 0.55;
Tst = 2;
ts = (0.00:0.001:Ts1); % msec
tsts = ((Ts1+0.001):0.001:Ts2);
tststs = ((Ts2+0.001):0.001:Tst);
t_sys = [ts tsts tststs];
Istim_1 = zeros(1,length(ts));
Istim_2 = 200.* ones(1,length(tsts));
Istim_12 = vertcat(Istim_1',Istim_2');
Istim_3 = zeros(1,length(tststs));
Istim_sys = vertcat(Istim_12,Istim_3');
end

My .m that I build, calling function above:

% Hodgkin-Huxley Model solved with ODE45 MATLAB function
tspan = [0:0.001:2];

% Initial Conditions t0 = 0 msec
Vm0 = 0; % mV
alpha_n0 = 0.01*(9.5-Vm0) / (-1+exp((9.5-Vm0)*0.1));
beta_n0 = 0.125*exp(-(Vm0)/80);
nss0 = alpha_n0 / (alpha_n0+beta_n0); % nss(0mV)
alpha_m0 = 0.1*(25-Vm0) / (-1+exp(0.1*(25-Vm0)));
beta_m0 = 4*exp(-(Vm0)/18);
mss0 = alpha_m0 / (alpha_m0+beta_m0); % mss(0mV)
alpha_h0 = 0.07*exp(-(Vm0)/20);
beta_h0 = 1 / (1+exp((30-Vm0)*0.1));
hss0 = alpha_h0 / (alpha_h0+beta_h0); % hss(0mV)
y0 = [Vm0; nss0; mss0; hss0];

% ODE45 plot
[t,y]=ode45(@dVm,tspan,y0);

% HH parameters
Vm = y(:,1);
n = y(:,2);
m = y(:,3);
h = y(:,4);

% a) Plot Voltage (sol of dV/dt) as function of time
plot(t,Vm);
title('Hodgkin-Huxley Model');
xlabel('Time (ms)');
ylabel('Potential Voltage (mV)');

Error:

Conversion to function_handle from double is not possible.

Error in dVm (line 14)
dVm(2,1) = (0.01*(0.5-y(1)) / (1-exp((9.5-y(1))/10)))*(1-y(2)) - (0.125*exp(-y(1)/ 80))*y(2); % dn/dt where y(2) = n

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in HHproject_a_b_c (line 18)
[t,y]=ode45(@dVm,tspan,y0);

Answer 1

[Star Strider]

You have not provided enough information to determine what the problem is. However, considering that ‘Istim_sys’ may be a vector that may be time-dependent, see the documentation section on ODE with Time-Dependnet Terms as a possible solution.

This is purely a guess, however my interpretation of the problem would be compatible with the error your code has thrown.

[Laura Zichella]

I fixed that and updated my code from my question posting, I am now stuck with:

Conversion to function_handle from double is not possible.

Error in dVm (line 14)
dVm(2,1) = (0.01*(0.5-y(1)) / (1-exp((9.5-y(1))/10)))*(1-y(2)) - (0.125*exp(-y(1)/ 80))*y(2); % dn/dt where y(2) = n

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in HHproject_a_b_c (line 18)
[t,y]=ode45(@dVm,tspan,y0);
 

[Star Strider]

Note that:

function Istim_sys

does not take any arguments nor return any outputs.

I suspect that is the problem. Since I have no idea what ‘Istim_sys’ should take as an argument or should return as an output, you will have to solve that. I cannot, given the currently available information.

[Laura Zichella]

How and Where do I return an output? Istim_sys is supposed to be the current that activates action potential in cells. The function is defined as 200 for 0.5<t<0.55 and 0 for other t and this is what I am making the function output, Istim_sys that reflects that, and can be used as an input into function dVm

[Star Strider]

The input argument would appear to be ‘t’. I have no idea what the output would be, since I do not understand the function you wrote. You may need to use interp1 to calculate the output.