Figuring out the stimulus of FitzHugh Nagumo equations with ode45

I am trying to write a program to solve the FitzHugh-Nagumo equations for a single cell with these equations below

du/dt =c1u(u −a)(1−u)−c2uv + stim

dv/dt = b (u − v)

where

a=0.13

b=0.013

c1=0.26

c2=0.1

stim is a stimulus current that can be applied for a short time at the beginning of the simulation.

u represents membrane potential and ranges from 0 (rest) to 1 (excited)

v is a recovery variable in the same range.

t is time in milliseconds.

I'm supposed to use MATLAB's ode45() function to integrate the system of differential equations. My program's input should include the duration of the simulation; initial values for u, v, and t; the strength of the stimulus, and the time for which it is applied (typically a few ms) and its output should include vectors for t, u and v.

function [t,u,v] = fhn0D(tf,u0,v0,t0,stim,tstim)

% Solve the FitzHugh-Nagumo equations for a single cell (i.e., without

% spatial coupling)

% du/dt = (c1*u)*(u - a)*(1 - u) ? (c2*u*v) + stim

% dv/dt = b(u - v)

% stim is a stimulus current that can be applied for a short time at

% the beginning of the simulation

% u represents membrane potential and ranges from 0 (rest) to 1 (excited)

% v is a recovery variable in the same range

% t is time in milliseconds

% INPUT: tf duration of simulation

% u0 initial value of u

% v0 initial value of v

% t0 initial value of t

% stim strength of stimulus

% tstim time of applied stimulus

%

% OUTPUT: t time vector (ms)

% u membrane potential vector

% v recovery vector

[t,U] = ode45(@(t,U)fhn(t,U,stim),[t0 tf],[u0; v0]);

% First and second column of U correspond to u and v respectively

u = U(:,1);

v = U(:,2);

% Plot u and v vs. t

plot(t,u,t,v)

title('Solution of FitzHugh-Nagumo Equations for a Single Cell with ODE45');

xlabel('Time t');

ylabel('Solution U');

legend('u','v');

end

function dUdt = fhn(t,U,stim)

% FitzHugh-Nagumo equations defined to be plugged into ode45

%

% INPUT: t time

% U vector that holds u and v

%

% OUTPUT: dUdt vector containing solutions to derivatives of u and v

% U(1) is u and U(2) is v

% Other needed variables

a = 0.13;

b = 0.013;

c1 = 0.26;

c2 = 0.1;

% Preallocate the output vector

dUdt = zeros(2,1);

% Create a two element vector that holds the derivative equations of u and

% v

dUdt = [(c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+ stim; b*(U(1)-U(2))];

end

I've posted the code I've created so far above, but I am confused on how to incorporate the time of stimulus in order to generate an action potential.

[darova]

What is wrong with your earlier question?

[Kate Heinzman]

function [t,u,v] = fhn0D(tf,u0,v0,t0,stim,tstim)

% Solve the FitzHugh-Nagumo equations for a single cell (i.e., without

% spatial coupling)

% du/dt = (c1*u)*(u - a)*(1 - u) ? (c2*u*v) + stim

% dv/dt = b(u - v)

% stim is a stimulus current that can be applied for a short time at

% the beginning of the simulation

% u represents membrane potential and ranges from 0 (rest) to 1 (excited)

% v is a recovery variable in the same range

% t is time in milliseconds

% INPUT: tf duration of simulation

% u0 initial value of u

% v0 initial value of v

% t0 initial value of t

% stim strength of stimulus

% tstim time of applied stimulus

%

% OUTPUT: t time vector (ms)

% u membrane potential vector

% v recovery vector

[t,U] = ode45(@(t,U)fhn(t,U,params),[t0 tf],[u0; v0]);

% First and second column of U correspond to u and v respectively

u = U(:,1);

v = U(:,2);

% Plot u and v vs. t

plot(t,u,t,v)

title('Solution of FitzHugh-Nagumo Equations for a Single Cell with ODE45');

xlabel('Time t');

ylabel('Solution U');

legend('u','v');

end

function dUdt = fhn(t,U,params)

% FitzHugh-Nagumo equations defined to be plugged into ode45

%

% INPUT: t time

% U vector that holds u and v

%

% OUTPUT: dUdt vector containing solutions to derivatives of u and v

% U(1) is u and U(2) is v

% Unpack other needed variables

a = params(1); % 0.13

b = params(2); % 0.013

c1 = params(3); % 0.26

c2 = params(4); % 0.1

stim = params(5);

tstim = params(6);

% Preallocate the output vector

dUdt = zeros(2,1);

% Only want the stimulus applied over a certain time point of the action

% potential not for the entire action potential

% If time is in the range of tstim then stim is incorporated into the

% equation otherwise stim is 0

if t <= tstim

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+ stim;

else

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2));

end

% Create a two element vector that holds the derivative equations of u and

% v

dvdt = b*(U(1)-U(2));

dUdt = [dudt; dvdt];

end

Sorry I understand the stimulus now, but instead of hardcoding the parameters needed into my function i used the vector params and MATLAB is giving me an error saying, Unrecognized function or variable 'params'.

[darova]

Maybe you forgot to assign values to variables in your function fhn0D?

I don't see any value. Where are they?

Answer 1

[Pravin Jagtap]

Hello Kate,

Refer to the following documentation link to understand how to solve system of differential equations. Also refer to the examples given in the link which will give you better insights.

https://www.mathworks.com/help/symbolic/solve-a-system-of-differential-equations.html#buxuui8

Hope this will help you.