# Source code for neurodynex3.phase_plane_analysis.fitzhugh_nagumo

"""
This file implements functions to simulate and analyze
Fitzhugh-Nagumo type differential equations with Brian2.

Relevant book chapters:

- http://neuronaldynamics.epfl.ch/online/Ch4.html
- http://neuronaldynamics.epfl.ch/online/Ch4.S3.html.

"""

# This file is part of the exercise code repository accompanying
# the book: Neuronal Dynamics (see http://neuronaldynamics.epfl.ch)
# located at http://github.com/EPFL-LCN/neuronaldynamics-exercises.

# This free software: you can redistribute it and/or modify it under
# Free Software Foundation. You should have received a copy of the
# GNU General Public License along with the repository. If not,

# Should you reuse and publish the code for your own purposes,
# please cite the book or point to the webpage http://neuronaldynamics.epfl.ch.

# Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski.
# Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition.
# Cambridge University Press, 2014.

import brian2 as b2
import matplotlib.pyplot as plt
import numpy as np

[docs]def get_trajectory(v0=0., w0=0., I_ext=0., eps=0.1, a=2.0, tend=500.):
r"""Solves the following system of FitzHugh Nagumo equations
for given initial conditions:

.. math::

\frac{dv}{dt} = \frac{1}{1ms} v (1-v^2) - w + I \\
\frac{dw}{dt} = eps (v + 0.5 (a - w))

Args:
v0: Intial condition for v [mV]
w0: Intial condition for w [mV]
I: Constant input [mV]
eps: Inverse time constant of the recovery variable w [1/ms]
a: Offset of the w-nullcline [mV]
tend: Simulation time [ms]

Returns:
tuple: (t, v, w) tuple for solutions
"""

eqs = """
I_e : amp
dv/dt = 1/ms * ( v * (1 - (v**2) / (mV**2) ) - w + I_e * Mohm ) : volt
dw/dt = eps/ms * (v + 0.5 * (a * mV - w)) : volt
"""

neuron = b2.NeuronGroup(1, eqs, method="euler")

# state initialization
neuron.v = v0 * b2.mV
neuron.w = w0 * b2.mV

# set input current
neuron.I_e = I_ext * b2.nA

# record states
rec = b2.StateMonitor(neuron, ["v", "w"], record=True)

# run the simulation
b2.run(tend * b2.ms)

return (rec.t / b2.ms, rec.v / b2.mV, rec.w / b2.mV)

[docs]def plot_flow(I_ext=0., eps=0.1, a=2.0):
"""Plots the phase plane of the Fitzhugh-Nagumo model
for given model parameters.

Args:
I: Constant input [mV]
eps: Inverse time constant of the recovery variable w [1/ms]
a: Offset of the w-nullcline [mV]
"""

# define the interval spanned by voltage v and recovery variable w
# to produce the phase plane
vv = np.arange(-2.5, 2.5, 0.2)
ww = np.arange(-2.5, 5.5, 0.2)
(VV, WW) = np.meshgrid(vv, ww)

# Compute derivative of v and w according to FHN equations
# and velocity as vector norm
dV = VV * (1. - (VV**2)) - WW + I_ext
dW = eps * (VV + 0.5 * (a - WW))
vel = np.sqrt(dV**2 + dW**2)

# Use quiver function to plot the phase plane
plt.quiver(VV, WW, dV, dW, vel)

[docs]def get_fixed_point(I_ext=0., eps=0.1, a=2.0):
"""Computes the fixed point of the FitzHugh Nagumo model
as a function of the input current I.

We solve the 3rd order poylnomial equation:
v**3 + V + a - I0 = 0

Args:
I: Constant input [mV]
eps: Inverse time constant of the recovery variable w [1/ms]
a: Offset of the w-nullcline [mV]

Returns:
tuple: (v_fp, w_fp) fixed point of the equations
"""

# Use poly1d function from numpy to compute the
# roots of 3rd order polynomial
P = np.poly1d([1, 0, 1, (a - I_ext)], variable="x")

# take only the real root
v_fp = np.real(P.r[np.isreal(P.r)])
w_fp = 2. * v_fp + a

return (v_fp, w_fp)