{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Dynamical systems with Python: Wilson and Cowan's model\n",
"\n",
"Wilson and Cowan, in the seventies, developed a dynamical systems approach to the study of the large-scale behaviour of neuronal population. Their approach doesn't mind the behaviour of single neurons, but works at the level of population firing rates (roughly, the number of neurons that fire in the unit time) for two subpopulation: the inhibitory neurons and the excitatory neurons.\n",
"\n",
"If you're interested in Wilson and Cowan's model, keep reading. If you came here to see some Python code, skip to the [numerical solution](#numerical).\n",
"\n",
"## The theory\n",
"\n",
"### Assumptions \n",
"\n",
"We are concerned only with the behaviour of **populations** rather than individual cells.\n",
"\n",
"> The cells comprising such populations are assumed to be in close spatial proximity, and their interconnections are assumed to be random, yet dense enough so that it is very probable that there will be at least one path (either direct or via interneurons) connecting any two cells within the population. Under these conditions we may neglect spatial interactions and deal simply with the temporal dynamics of the aggregate. Consistent with this approach, we have chosen as the relevant variable the proportion of cells in the population which become active per unit time. [...]\n",
"\n",
"> There is one final and crucial assumption upon which this study rests: *all nervous processes of any complexity are dependent upon the interaction of excitatory and inhibitory cells*.\n",
"\n",
"Spatial interactions are discussed in a subsequent paper by the same authors.\n",
"\n",
"\n",
"### Getting to a system of equations\n",
"\n",
"In the rest of the post, we will call $E(t)$ and $I(t)$ the istantaneous firing rates, at time $t$, of the excitatory and inhibitory populations respectively. The point $(E=0, I=0)$ should not be interpreted as the state where the activity is completely dead, but rather as the resting state of the network, as we will require it to be a stable fixed point.\n",
"\n",
"To construct a meaningful model, we slowly incorporate the relevant biological assumptions. When does a neuron fire? In the simplest models, two conditions need to be fulfilled:\n",
"1. the neuron should *not* be in its \"refractory period\", that is, it can't fire again just after having fired;\n",
"2. it needs to have received sufficient input in a short interval of time.\n",
"\n",
"##### Non-refractoriness\n",
"The fraction of excitatory neurons that fired between $t_1$ and $t_2$ is\n",
"$$ \\int_{t_1}^{t_2} E(t')dt'. $$ \n",
"Therefore, if $r$ is the length of the refractory period, the fraction of neurons that satisfy condition 1 at time $t$ is: \n",
"$$ 1 - \\int_{t-r}^{t} E(t')dt' \\quad\\quad (1)$$ \n",
"and similarly for the inhibitory subpopulation.\n",
"\n",
"##### Sufficient excitation\n",
"To see if condition 2 is fulfilled, we need the total input to the subpopulation, which is $c_1 E(t) - c_2 I(t) + P(t)$, i.e. a weighed contribution from the excitatory population, a corresponding negative contribution from the inhibitory neurons, and an external input $P$. Then, the neuron responds non-linearly to this input. We call $S_e$ the response function, also called the *input-frequency characteristic* of the excitatory neurons, and correspondingly $S_i$ for the inhibitory ones. \n",
"Now, it's not only the instantaneous behaviour that counts: a spike can still help eliciting a new spike in a downstream neuron even a few milliseconds later. So the probability of being excited at time $t$ is proportional to\n",
"$$ S_e \\left( \\int_{-\\infty}^t \\alpha(t-t') [c_1 E(t') - c_2 I(t') + P(t')]dt' \\right) \\quad\\quad (2)$$\n",
"\n",
"\n",
"##### Coarse graining\n",
"Both in (1) and (2) we can get rid of the integrals and multiply the stimulus by a constant describing the length of the time influence instead, if we are interested in the coarse grained temporal behaviour of the activity, i.e. we focus on variations at a timescale slightly longer than the refractory period $r$ and the \"characteristic length\", which we call $k$, of the function $\\alpha$. So, finally, we can say the activity at time $d + dt$ depends on the simultaneous fulfillment of conditions (1) and (2):\n",
"$$ E(t + dt) = (1- rE(t))\\, S_e(kc_1 E(t) - kc_2 I(t) + kP(t)) $$\n",
"\n",
"\n",
"### Wilson and Cowan's model, a final version\n",
"\n",
"After all these approximations and assumptions, by turning the equation above in differential form and appropriately rescaling $S_e$, we reach a system of coupled, nonlinear, differential equation for the firing rates of the excitatory and inhibitory populations, which constitute the Wilson-Cowan model.\n",
"\n",
"$$ \\tau_e \\frac{dE}{dt} = -E + (k_e - r_e E) \\, S_e(c_1 E - c_2 I + P)$$\n",
"$$ \\tau_i \\frac{dI}{dt} = -I + (k_i - r_i I) \\, S_i(c_3 E - c_4 I + Q),$$\n",
"where:\n",
"* $\\tau_e$ and $\\tau_i$ are time constants;\n",
"* $k_e$ and $k_i$ are adimensional constants\n",
"* $r_e$ and $r_i$ are constants describing the length of the refractory periods;\n",
"* $S_e$ and $S_i$ are sigmoid functions expressing the nonlinearity of the interactions;\n",
"* $c_{1,2,3,4}$ are parameters representing the strength of the excitatory to excitatory, inhibitory to excitatory, excitatory to inhibitory and inhibitory to inhibitory interactions;\n",
"* $P$ and $Q$ are external inputs to the two populations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical solutions\n",
"\n",
"Now, we want to numerically solve the [equations above](#numerical), regardless of their meaning. We will choose some values of their parameters, and some initial conditions, and plot the behaviour of the solution in the state space and in time."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# for fast array manipulation\n",
"import numpy as np\n",
"# for plotting\n",
"import matplotlib.pyplot as plt\n",
"# for numerical ODE integration\n",
"from scipy.integrate import odeint\n",
"# for nonlinear equations\n",
"from scipy.optimize import fsolve\n",
"# to display plots in-line\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parameter definitions\n",
"To represent the nonlinear behaviour of neurons, we use a sigmoid function, dependent on two parameters $a$ and $\\theta$:\n",
"\n",
"$$ S(x) = (1+e^{-a(x-\\theta)})^{-1} $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def sigmoid(x, a, thr):\n",
" return 1 / (1 + np.exp(-a * (x - thr)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### For a stable limit cycle and a stable fixed point"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# couplings\n",
"c1 = 16\n",
"c2 = 12\n",
"c3 = 15\n",
"c4 = 3\n",
"\n",
"# refractory periods\n",
"rE = 1\n",
"rI = 1\n",
"\n",
"# external inputs\n",
"P = 1.\n",
"Q = 1\n",
"\n",
"# nonlinear functions\n",
"def Se(x):\n",
" aE = 1.3\n",
" thrE = 4\n",
" return sigmoid(x, thrE, aE) - sigmoid(0, thrE, aE)\n",
"\n",
"def Si(x):\n",
" aI = 2\n",
" thrI = 3.7\n",
" return sigmoid(x, thrI, aI) - sigmoid(0, thrI, aI)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### More preparation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# this function returns the right hand side of the Wilson-Cowan equation\n",
"# (both, in a 2-vector)\n",
"def WilsonCowan(y, t):\n",
" E = y[0]\n",
" I = y[1]\n",
" y1 = -E + (1 - rE * E) * Se(c1 * E - c2 * I + P)\n",
" y2 = -I + (1 - rI * I) * Si(c3 * E - c4 * I + Q)\n",
" return [y1, y2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# minimum and maximum E and I values we want displayed in the graph\n",
"minval = -.1\n",
"maxval = .6\n",
"resolution = 50\n",
"# State variables\n",
"x1 = np.linspace(minval, maxval, resolution)\n",
"x2 = np.linspace(minval, maxval, resolution)\n",
"# Create a grid for evaluation of the vector field\n",
"x1, x2 = np.meshgrid(x1, x2)\n",
"# Evaluate the slopes\n",
"X1, X2 = WilsonCowan([x1, x2], 0)\n",
"# Compute the magnitude vector\n",
"M = np.hypot(X1, X2)\n",
"# Normalize the slopes vectors (for the field plot)\n",
"#X1, X2 = X1/M, X2/M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving and plotting\n",
"#### Numerical search of stationary points\n",
"\n",
"This is a tricky part. To view the stationary points, which are an important property of a dynamical system, we need to study the points where the derivatives of $E$ and $I$ are zero, i.e., where the function called WilsonCowan is zero. Since it's highly nonlinear, we have to do that numerically.\n",
"\n",
"Numerical root finding basically works by taking an initial guess, then following the shape of the function until we reach a zero. If we want *all* of the zeros, we need to try many initial guesses, and we are not guaranteed to succeed.\n",
"\n",
"In practice, it seems fast enough to use all the points of the grid I'll use for plotting, and compute a zero for each of them. But this may not always work."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"fixed_p = []\n",
"y1 = x1.ravel()\n",
"y2 = x2.ravel()\n",
"for i in range(resolution**2):\n",
" # find a zero\n",
" sol, infodict, ier, mesg = fsolve(WilsonCowan, [y1[i], y2[i]], args=(0), full_output=1)\n",
" if ier == 1: # I exclude the cases where fsolve didn't converge\n",
" fixed_p.append(sol)\n",
"\n",
"fixed_p = np.array(fixed_p).T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Numerical ODE integration\n",
"\n",
"Here is where we actually integrate the dynamical system in time. **Try changing the values of E0 and I0** to obtain different paths in the phase space."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# simulation duration and step size\n",
"time = np.linspace(0, 100, 2000)\n",
"\n",
"# starting point, hopefully inside the basin of attraction of our attractor\n",
"E0, I0 = 0.39, 0.49 # try changing this\n",
"\n",
"# find the solution with scint.odeint\n",
"odesol = odeint(WilsonCowan, [E0, I0], time)\n",
"# separate the two solutions\n",
"exc_timeseries, inh_timeseries = odesol.T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Plotting"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [],
"source": [
"# plotting the vector field in the state space (E, I)\n",
"plt.figure(figsize=(10, 10))\n",
"plt.quiver(x2, x1, X2, X1, pivot='mid', alpha=.5)\n",
"plt.xlim([minval, maxval])\n",
"plt.ylim([minval, maxval])\n",
"plt.xlabel(r'$I$', fontsize=16) # yes, you can use Latex code!\n",
"plt.ylabel(r'$E$', fontsize=16)\n",
"plt.grid()\n",
"\n",
"# plot the solution in the state space\n",
"plt.plot(inh_timeseries, exc_timeseries, '.-');\n",
"\n",
"# plot the starting point\n",
"plt.scatter(I0, E0, marker='*', s=300, label=\"Starting point\")\n",
"plt.legend(loc=\"upper left\")\n",
"\n",
"# plot the fixed points we identified\n",
"plt.scatter(fixed_p[1], fixed_p[0], marker='o', s=50, label=\"Stationary points\")\n",
"\n",
"# plot the solution in time\n",
"plt.figure(figsize=(10.3,3))\n",
"plt.ylabel(r'$E, I$')\n",
"plt.xlabel(r'$t$')\n",
"plt.plot(time, exc_timeseries, '.-', label=\"excitatory\");\n",
"plt.plot(time, inh_timeseries, '.-', label=\"inhibitory\");\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see the system exhibits:\n",
" - a **stable** stationary point at (0.5, 0.5)\n",
" - an **unstable** stationary point at (0.5, 0.4)\n",
" - a **limit cycle** about the stationary point at (0.1, 0.05)\n",
" \n",
"The different nature of the stationary point, and the consequent influence they have on the behaviour of the solution in their neighbourhood, is due to the eigenvalues of the system when linearised near that point. Refer to a textbook for that!"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.5"
}
},
"nbformat": 4,
"nbformat_minor": 0
}