""" File: fhnneurons.py Author: Martin Sauer Email: sauer@math.tu-berlin.de Date: May 20, 2014 Simulation of single (stochastic) FHN neurons and networks with disorder dv/dt = kappa v (v-a) (1-v) - w + I dw/dt = epsilon (v - gamma w + b) where I is the applied current or input signal, e.g. white noise. """ import numpy as np from numpy import random as rnd from matplotlib import pyplot as plt class FitzHughNagumoNeuron: """ Class for the FitzHugh-Nagumo neuron. The init_state is [kappa, a, epsilon, gamma, b, sigma, v0, w0] where kappa : scaling of nonlinearity, a : second root of nonlinearity, epsilon : speed of the recovery variable, gammma : inverse slope of the recovery nullcline, b : shift of the recovery nullcline sigma : the noise intensity, v0, w0 : initial conditions. """ def __init__(self, init_state = [4, 0.5, 0.08, 0.8, -0.29794, 0.05, -1, -2], ): self.params = np.asarray(init_state[0:6], dtype='float') self.state = np.asarray(init_state[6:8], dtype='float') self.time_elapsed = 0 def timestep(self, dt, I, v_noise, w_noise): [kappa, a, epsilon, gamma, b, sigma] = self.params [v, w] = self.state self.state[0] = (v + dt * ( kappa * v * (1 - v) * (v - a) - w + I) + sigma * np.sqrt(dt) * v_noise) self.state[1] = (w + dt * epsilon * (v - gamma * w + b) + sigma * np.sqrt(dt * epsilon) * w_noise) self.time_elapsed += dt class FitzHughNagumoNetwork: """ Class for a network of FitzHugh-Nagumo neurons. The init_state is [N, v0, w0, Q] where N : number of neurons, v0 : vector of initial values for v, w0 : vector of initial values for w, Q : connection matrix, sigma : common noise intensity. """ def __init__(self, init_state = [2, -1 * np.ones(2), -2 * np.ones(2), np.ones(2) - 2 * np.eye(2), 0.] ): self.size = init_state[0] self.neurons = [] for i in xrange(self.size): self.neurons.append(FitzHughNagumoNeuron([1, 0.25, 0.01, 0.95 * 3 / (1 - 0.25 + 0.25**2), 0, init_state[4], init_state[1][i], init_state[2][i]])) self.coupling = np.asarray(init_state[3], dtype='float') self.time_elapsed = 0 def timestep(self, dt, I, v_noise, w_noise): """time evolution of the network""" interaction = [] for i in xrange(self.size): interaction.append(self.neurons[i].state[0]) interaction = np.dot(self.coupling, interaction) + I for i in xrange(self.size): self.neurons[i].timestep(dt, interaction[i], v_noise[i], w_noise[i]) def position(self): state = [] for i in xrange(self.size): state.append(self.neurons[i].state) return state def v_position(self): state = [] for i in xrange(self.size): state.append(self.neurons[i].state[0]) return state """ Setting for global parameters """ dt = 0.1 T = 14000 time = np.linspace(0, T * dt, T) I = 0.06 theta_thr = (1 - 0.25 + 0.25**2) / 6 theta = 0.01 * theta_thr sigma = 0.0 fig = plt.figure() ax = plt.axes(xlim=(0,T * dt), ylim=(-1,2)) neurons = FitzHughNagumoNetwork([2, [-1,1], [-0.25, -0.5], theta * (np.ones(2) - 2 * np.eye(2)), sigma]) v_data = np.zeros([2,T]) for t in xrange(T): neurons.timestep(dt, I, rnd.normal() * np.ones(2), rnd.normal() * np.ones(2)) v_data[:,t] = neurons.v_position() lines = [ax.plot(time, v_data[i], lw=2) for i in xrange(neurons.size)] plt.show()