# How can I calculate values for this system of differential equations?

Disclaimer: I'm aware this is likely basic math for many on this forum, but I'm lacking in the rudimentary maths needed to even... ask this properly, so I'll just go for it.

For my own curiosity, I'm studying several papers about modelling biological neurons, and am stuck on plotting the FitzHugh-Nagumo Model:

\begin{matrix} \dot{V} & = & V-V^3/3 - W + I \\ \dot{W} & = & 0.08(V+0.7 - 0.8W) \end{matrix}

I want to know how V and W evolve over time, where V represents a neuron's "voltage", and W is a "recovery variable" that is needed to model a neuron's firing pattern (and I is a constant, representing a constant current input into the neuron).

As far as I know, V' is just the change in V over time, so I've tried running the model considering it's evolution over time as V + V', but I think that's too simplistic, and I certainly am not able to generate the plots I'm seeing in the relevant papers.

What do I need to know to get these differential equations into the right form, IE V = ... instead of V' = ...?

Specifically, I'd love to be able to generate:

• phase diagrams
• and the evolution of V (or W) over time
• What software do you have available? Commented Dec 27, 2017 at 19:38
• @dbx anything open source, I'm using Haskell, R, python, perhaps Octave (like matlab). I'm comfortable using most anything if I can understand the math. Commented Dec 27, 2017 at 19:40
• There's not really a simple answer. It will depend on what you mean by "the math" and is beyond the scope of a stack exchange comment. Differential equations is the field. Commented Dec 27, 2017 at 19:50
• @Josh.F: Here are some notes that might prove helpful www4.ncsu.edu/~msolufse/LectureFitzHughNagumo.pdf
– Moo
Commented Dec 28, 2017 at 5:59
• You can plot a rough phase portrait in Wolfram Alpha. Commented Dec 28, 2017 at 15:59

Here's a numerical solution I was able to work up in Octave:

function result = f(u,t)
curr = 1;

V = u(1); W = u(2);

result = [V-V^3/3-W+curr 0.08*(V+0.7-0.8*W)]
end;

function demo()
V_0 = 0.7;
W_0 = 0.7;

f_handle = @f;

t = (0:0.1:100);
u = lsode(f_handle, [V_0 W_0], t);

V = u(:,1);
W = u(:,2);

plot(t, V, 'b', t, W, 'r');
figure;

end;


The workhorse is lsode that, as far as I know, just starts from the initial conditions of V_0, W_0 and evolves the dynamic system forward.