# Modeling neural excitation with a non-linear differential equation

I think I have a special problem and I am not sure how to search for an answer, so I thought I would try here.

I am working with the so called FitzHugh-Nagumo model which describes very simple neurons by the following equations:

u' = v - 1/3*v^3 - u
v' = u + a


where a is the excitation threshold.

If a is greater than 1, the system is in the so called excited state. So there is one fix point in which all trajectories will lead for t -> ∞. With the NDSolve function there is no problem computing this for every initial condition I want, e.g., the following:

s = NDSolve[{ε*u'[t] == u[t] - 1/3*u[t]^3 - v[t], v'[t] == u[t] + a,
u[0] == 0.5, v[0] == 0.5}, {u, v}, {t, 200}, MaxSteps -> 10^6]


If the system now is excited a little, nothing special happens, it just relaxes back to the fixpoint again. But if the excitation is large, the system can make a loop before it relaxes. And that would be exactly a neural pulse.

So my problem is how to tell Mathematica that it should excite the system again after it reaches its fix point. I already found the EvenLocator function, which gives me the t value, when u or v is close to the fixed point:

Method -> {"EventLocator", "Event" -> u'[t] - 0.00001, "EventAction" -> ?}


with the value in Event I can control how close my system is at a fix point. Now I want the dynamic variables u and v to jump to a new point above the threshold and start to move through the dynamics again. So here are the questions:

• Can I use the EventLocator for what I want and is EventAction the right instrument?
• What possibilities do I have to simulate such excitation (something continuous would be nice)?
• Or am I totally on wrong direction?

Maybe it is important to state that I use Mathematica 6.0 at the moment.

Update:

thanks for your help! Your solution is very interesting for me because I also want to solve the problem in more Dimensions, but in my special case in 0D I have another problem: There is no travelling pulse that excitates element by element. So I need another trick to simulate something like this. But studying your solution I had an idea: because there is no diffusion (of course!) I need something like a virtual excitation from the neighbors to restart the pulse. So maybe it is possible to simulate just two neighbors and have something like very small 1D model. So the more dimensional approach maybe can help me in terms of simulating such a "virtual diffusion". Do you know what I mean?

The point is, that this is only interesting if it would save resources well compared to a full 1D model. Otherwise it probably would not make sense to have a 0D model ...

What do you think?

-
1) The formula u'=v-1/3*v^3-u is not the same as epsilon u'[t] == u[t] - 1/3*u[t]^3 - v[t] 2) can you give example of value for epsilon (+1 -1 ... ?) – andre Jan 28 '13 at 19:54
Related question: mathematica.stackexchange.com/questions/1511/… – Michael E2 Mar 31 '13 at 12:49

ε = .1;
a = 5;
trip = -2 a {1, 1, 1};
dt = .2;
s = NDSolve[{ε*u'[t] == -u[t] - 1/3*v[t]^3 + v[t],
v'[t] == u[t] + a, u[0] == 0.5, v[0] == 0.5,
WhenEvent[Mod[t, dt] == 0,
trip = Join[{u[t]}, Rest@RotateRight@trip];
If[Norm[trip + a] < .1, u[t] -> u[t] + a, u[t] -> u[t]]]
}, {u, v}, {t, 10}, MaxSteps -> 10^6];
Plot[{u[t], v[t]} /. s[[1]], {t, 0, 10}]


-