solving 2D FitzHugh-Nagumo system

function [U,V,t] = fhn1d(tend,dx,stimMask,stimStr,stimDur,U0,V0,t0)

% function [U,V,t] = fhn1d(tend,dx,stimMask,stimStr,stimDur,U0,V0,t0)

% simulates propagation of a FitxHugh-Nagumo cable

%

% INPUT: tend time to integrate to

% dx spatial time step in nm

% stimMask length of this vector specifies how many nodes are in

% the cable and specifies which nodes will receive a

% stimulus current

% stimStr strength of current stimulus

% stimDur duration of current stimulus

% U0 vector the same length of stimMask that contains the

% initial condition for u on each node

% V0 vector the same length of stimMask that contains the

% initial condition for v on each node

% t0 time at which to begin integrating

%

% OUTPUT: U contains solution where it has a column for each node on

% the cable and a row for each timestep

% V contains solution where it has a column for each node on

% the cable and a row for each timestep

% t vector with the time of each timestep

% Set argument defaults

% The first three arguments are required (others are optional)

if nargin < 8 || isempty(t0)

t0 = 0;

end

if nargin < 7 || isempty(V0)

V0 = zeros(size(stimMask));

end

if nargin < 6 || isempty(U0)

U0 = zeros(size(stimMask));

end

if nargin < 5 || isempty(stimDur)

stimDur = 0;

end

if nargin < 4 || isempty(stimStr)

stimStr = 0;

end

% Set membrane parameters and cable conductivity

a = 0.13;

b = 0.013;

c1 = 0.26;

c2 = 0.1;

D = 5;

% stimMask is a column vector

stimMask = stimMask(:);

% Number of nodes

nx = length(stimMask);

% Stack up initial values for u and v into a long column vector

% All u values are first and then the v values

UV0 = [U0(:);V0(:)];

% Set time at which stimulus turns off

stimOff = stimDur + t0;

% Make a handle to an anonymous function that will evaluate the

% right-hand-sides of all the ODEs in the system

fhnHandle = @(t,UV) fhn1(t,UV,stimStr,stimOff,stimMask,dx,a,b,c1,c2,D);

% Call ode 45 to do the integration

% On output, t is a column vector of times

% UV has one column for each ODE in the system and one row per timestep

% fhn1 is written so that the first nx columns in UV contain u values going

% down the cable and the next nx columns contain the v values

[t,UV] = ode45(fhnHandle,[t0, tend],UV0);

% Unpack UV into seperate U and V arrays

U = UV(:,1:nx);

V = UV(:,nx+1:end);

plotFiber(U,V,t)

end

function dUVdt = fhn1(t,UV,stimStr,stimOff,stimMask,dx,a,b,c1,c2,D)

% function dUVdt = fhn1(t,UV,stimStr,stimOff,stimMask,dx,a,b,c1,c2,D)

% evaluates the right-hand-sides of a system of ODEs that stimulate

% propagation on an FHN cable; called by ode45()

%

% INPUT: t current time

% UV column vector with the current values of the dependent

% variables

% stimStr strength of current stimulus

% stimOff time at which the stimulus switches off

% stimMask vector that specifies which nodes receive the stimulus

% dx spatial time step in nm

% a,b,c1,c2,D model parameters

%

% OUTPUT: dUVdt column vector containing solutions

% Get number of nodes

nx = length(stimMask);

% Unpack u and v

u = UV(1:nx);

v = UV(nx+1:end);

% Compute the diffusive term for each node, d2u/dx2, and store in column

% vector L

% Preallocate L

L = zeros(nx,1);

% Special term for first node that imposes no-flux boundary

L(1) = D*(u(2)-u(1))/(dx*dx);

% Diffusive term for all interior nodes

for i = 2:nx-1

L(i) = D*(u(i-1) + u(i+1) - 2*u(i))/(dx*dx);

end

% No-flux boundary condition for last node

L(nx) = D*(u(nx-1)-u(nx))/(dx*dx);

% Make a vector that contains the current stimulus for each node

if t < stimOff

IStim = stimMask*stimStr;

else

IStim = stimMask*0;

end

% Get fhn current for each node

I = c1.*u.*(u-a).*(1-u)-c2*u.*v;

% Add all the terms together to give dUdt on each node

dUdt = L + I + IStim;

% Compute dVdt

dVdt = b*(u-v);

% Stack all of our dUdt and dVdt values into one column vector

dUVdt = [dUdt;dVdt];

end

How do I take my code for the solving the fitzhugh-nagumo system in one dimension and solve it two dimensions according to the equations below?

The membrane parameters (a,b,c1,c2) are the same. D11 and D22 are conductivity constants in the x and y directions respectively; both can be set to five. How do you specify the number of nodes in the sheet in both directions? The spatial steps dx and dy should be input arguments to the function. How do you do the no flux boundary conditions for d2u/dx2 on the left and right edges of the sheet and the top and bottom edges for d2u/dy2? Then once the diffusive terms are figured out for this system, how do you pass them into ode45? Overall, I want to know how to convert my 1d version into what needs to be solved for the 2d version.

[Kate Heinzman]

The program should output a time vector and separate 3D arrays for u and v. The first dimension of the u and v arrays should be time and the second and third should be the two spatial dimensions.

[darova]

Why are you producing the same questions? I have answered you. Didn't my answer help you?

[Kate Heinzman]

The answer helped for the 0-dimension version, but not for the 2D version I am asking about in this question. This question is not the same as my others.

[darova]

I already have answered you about your difference scheme. Did you see that?

LINK

[Kate Heinzman]

Yes, but this solution did not help me. Can you please read over this current problem that I am asking about? I am not concerned with the difference scheme; I am confused on how to pass the data to ode45 such that it returns a 3D array where the first layer is time and second and third layers are the two spatial dimensions. What I am asking is how to take the correct 1 dimensional code I already have posted above and update it to a 2 dimensional version.

Answer 1

[darova]

• How do you do the no flux boundary conditions for d2u/dx2 on the left and right edges of the sheet and the top and bottom edges for d2u/dy2?

It's question to you. What are boundary conditions (temperature or heat flux) ?

• Then once the diffusive terms are figured out for this system, how do you pass them into ode45? Overall, I want to know how to convert my 1d version into what needs to be solved for the 2d version.

Once you have "2D version" you can't use ode45

I attach similar problem solved using difference scheme. See if it helps

[Kate Heinzman]

I've been tasked with using ode45 to solve it; I can't use the difference scheme. Boundary conditions are the same as in the above code except they have to be applied for d2u/dx2 on the left and right columns of the sheet, d2u/dy2 on the top and bottomo row of sheet, and both d2u/dx2 and d2u/dy2 on the corners since this is where they overlap.

[darova]

ode45 (Ordinary Differencial Equation) is for ordinary diff equations

You have partial differential equations (function depends on 2 or more variables)

You can't solve your equations using ode45