AT2 – Neuromodeling: Problem set #4 NETWORKS

PROBLEM 2: Circuit with mutual inhibition

We now consider a circuit of two neurons that are coupled by mutual inhibition. For the firing rates of these two neurons, $x_1$ and $x_2$, we assume the differential equations:

$$ \begin{cases} \dot{x}_1(t) = -x_1(t) + f(w x_2(t) + I)\\ \dot{x}_2(t) = -x_2(t) + f(w x_1(t) + I)\\ \end{cases} $$

where $f$ is defined as before and the inhibitory synaptic weights are given by $w = -0.1$. The external inputs are assumed to be excitatory, $I = 5$.

In [1]:
# /!\ I'm using Python 3

import plotly
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go

import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import stats
%matplotlib inline

plotly.offline.init_notebook_mode(connected=True)
from IPython.core.display import display, HTML, Markdown
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
display(HTML(
    '<script>'
        'var waitForPlotly = setInterval( function() {'
            'if( typeof(window.Plotly) !== "undefined" ){'
                'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
                'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
                'clearInterval(waitForPlotly);'
            '}}, 250 );'
    '</script>'
))

# Colorscales
def colorscale_list(cmap, number_colors, return_rgb_only=False):
    cm = plt.get_cmap(cmap)
    colors = [np.array(cm(i/number_colors)) for i in range(1, number_colors+1)]
    rgb_colors_plotly = []
    rgb_colors_only = []
    for i, c in enumerate(colors):
        col = 'rgb{}'.format(tuple(255*c[:-1]))
        rgb_colors_only.append(col)
        rgb_colors_plotly.append([i/number_colors, col])
        rgb_colors_plotly.append([(i+1)/number_colors, col])
    return rgb_colors_only if return_rgb_only else rgb_colors_plotly

import plotly.figure_factory as ff

def formatted(f): 
    return format(f, '.2f').rstrip('0').rstrip('.')

(a) Plot the null-isoclines ("nullclines") of the system, i.e. the line for which $\dot{x}_1(t) = 0$ and the line for which $\dot{x}_2(t) = 0$. What do the crossing points of these lines indicate?

In [2]:
w = -.1
I = 5

def f(s):
    return 50*(1+np.tanh(s))

f = np.vectorize(f)
In [3]:
x_min, x_max = 0, 126
xs = np.arange(x_min, x_max, 5)

x_1, x_2 = np.meshgrid(xs, xs)

u, v = -x_1+f(w*x_2+I), -x_2+f(w*x_1+I)
In [4]:
# Create quiver figure
fig1 = ff.create_quiver(x_1, x_2, u, v,
                       scale=.11,
                       arrow_scale=.2,
                       line=dict(width=1.1),
                       name='$\\dot{\\vec{x}}$')


layout1 = go.Layout(
    title= '2D flow of the system dynamics (arbitrary units)',
    hovermode= 'closest',
    xaxis= dict(
        title= '$x_1$',
        ticklen= 5,
        zeroline= False,
        range=[xs[0]-5, xs[-1]+5]
    ),
    yaxis=dict(
        title= '$x_2$',
        ticklen= 5,
        zeroline= False,
        range=[xs[0]-5, xs[-1]+5]
    ),
    legend=dict(
        orientation='h',
        y = -.1
    ),
    height = 800
) 

fig1['layout'] = layout1

fixed_points_mask = (np.abs(u) < .1) * (np.abs(v) < .1)
x_fp, y_fp = np.where(fixed_points_mask)

# Fixed points
color_fixed_points = 'ForestGreen'
stable_pts = np.array([True, False, True])

stable = go.Scatter(
    x = x_1[x_fp[stable_pts], y_fp[stable_pts]], 
    y = x_2[x_fp[stable_pts], y_fp[stable_pts]],
    mode='markers',
    marker=dict(size=10, color=color_fixed_points),
    name='stable fixed points'
)

unstable = go.Scatter(
    x = x_1[x_fp[~stable_pts], y_fp[~stable_pts]], 
    y = x_2[x_fp[~stable_pts], y_fp[~stable_pts]],
    mode='markers',
    marker=dict(
        size = 8,
        color = 'transparent',
        line = dict(
            color = color_fixed_points,
            width = 2
        )
    ),
    name='unstable fixed points'
)

# Add points to figure
fig1['data'].extend([stable, unstable])

plotly.offline.iplot(fig1)
In [5]:
x1_nullcline_color = 'DarkViolet'
x2_nullcline_color = 'Goldenrod'

nullcline_x1, nullcline_x2 = [], []

nullcline_x1.append(go.Contour(
    x = xs,
    y = xs,
    z = u,
    contours=dict(
        start = 0,
        end = 0,
    ),
    colorscale = [[0, 'transparent'], [1, 'transparent']],
    line=dict(
        color = x1_nullcline_color,
        width = 3,
        smoothing=1.3,
        dash = 'longdash'
    ),
    name = '$x_1\\text{-nullcline}$',
    showscale=False,
    showlegend=True,
))


nullcline_x2.append(go.Contour(
    x = xs,
    y = xs,
    z = v,
    contours=dict(
        start = 0,
        end = 0,
    ),
    colorscale = [[0, 'transparent'], [1, 'transparent']],
    line=dict(
        color = x2_nullcline_color,
        width = 3,
        smoothing=1.3,
        dash = 'longdash'
    ),
    name = '$x_2\\text{-nullcline}$',
    showscale=False,
    showlegend=True,
))

# Workaround to add legends, because legends aren't supported for contour traces in plotly yet
nullcline_x1.append(go.Scatter(
    x = [0],
    y = [0],
    mode = 'lines',
    line = dict(
        color = x1_nullcline_color,
        width = 3,
        dash = 'dash'
    ),
    name = '$x_1\\text{-nullcline}$',
    showlegend=True
))

nullcline_x2.append(go.Scatter(
    x = [0],
    y = [0],
    mode = 'lines',
    line = dict(
        color = x2_nullcline_color,
        width = 3,
        dash = 'dash'
    ),
    name = '$x_2\\text{-nullcline}$',
    showlegend=True
))

flow_nullclines = fig1['data']+nullcline_x1+nullcline_x2

fig1['layout']['title'] = 'Phase portait of the system dynamics, with nullclines (arbitrary units)'

plotly.offline.iplot(go.Figure(data=flow_nullclines, layout=fig1['layout']))

(b) Simulate the system (time step $dt = 0.1$)! Choose an initial condition for $x_1(0)$ and $x_2(0)$ and plot the evolution of the firing rates of the two neurons into the same plot as the nullclines from (a). What happens for different initial conditions?

In [6]:
delta_t = 0.1
T = 10
In [7]:
def evolution(x1_0, x2_0, delta_t=delta_t, T=T):
    x = np.zeros((len(np.linspace(0, T, int(T//delta_t))), 2))
    x[0] = x1_0, x2_0
    
    for i in range(1, x.shape[0]):
        x[i] = (1-delta_t)*x[i-1] + delta_t*f(w*x[i-1,::-1] + I)
    return x
In [8]:
x1_0, x2_0 = 1, 1
x = evolution(x1_0, x2_0) 

trace_updated = []
colors = colorscale_list('Reds', number_colors=int(1.5*len(x)), return_rgb_only=True)
legend_step = len(x)//5

for i, updated_pos in enumerate(x):
    trace_updated.append(go.Scatter(
        x = [x[i][0]],
        y = [x[i][1]],
        mode = 'markers',
        hovertext = 'time: {}'.format(formatted(i*delta_t)),
        hoveron = "points",
        marker = dict(
            size = min(7+i//legend_step, 18),
            color = colors[i],
            line = dict(
                width = 2,
                color = 'Maroon'
            )
        ),
        name = 'At time {}'.format(formatted(i*delta_t)),
        showlegend = (i % legend_step == 0)
    ))

trace_identity = go.Scatter(
    x = xs,
    y = xs,
    mode = 'line',
    line = dict(
        width = 2,
        color = 'black'
    ),
    opacity = 0.5,
    name = 'Identity line'
)

fig1['layout']['title'] = '$\\text{Simulation of the system with initial condition } (x_1(0), x_2(0)) = '+\
'({}, {})'.format(x1_0, x2_0)+' \\text{ (arbitrary units)}$'

fig = dict(data=flow_nullclines + trace_updated + [trace_identity], layout=fig1['layout'])
plotly.offline.iplot(fig)
In [9]:
x1_0, x2_0 = 1, 2
x = evolution(x1_0, x2_0) 

trace_updated = []
colors = colorscale_list('Reds', number_colors=int(1.5*len(x)), return_rgb_only=True)
legend_step = len(x)//5

for i, updated_pos in enumerate(x):
    trace_updated.append(go.Scatter(
        x = [x[i][0]],
        y = [x[i][1]],
        mode = 'markers',
        hovertext = 'time: {}'.format(formatted(i*delta_t)),
        hoveron = "points",
        marker = dict(
            size = min(7+i//legend_step, 18),
            color = colors[i],
            line = dict(
                width = 2,
                color = 'Maroon'
            )
        ),
        name = 'At time {}'.format(formatted(i*delta_t)),
        showlegend = (i % legend_step == 0)
    ))

trace_identity = go.Scatter(
    x = xs,
    y = xs,
    mode = 'line',
    line = dict(
        width = 2,
        color = 'black'
    ),
    opacity = 0.5,
    name = 'Identity line'
)


fig1['layout']['title'] = '$\\text{Simulation of the system with initial condition } (x_1(0), x_2(0)) = '+\
'({}, {})'.format(x1_0, x2_0)+' \\text{ (arbitrary units)}$'

fig = dict(data=flow_nullclines + trace_updated + [trace_identity], layout=fig1['layout'])
plotly.offline.iplot(fig)
In [10]:
x1_0, x2_0 = 1, 0
x = evolution(x1_0, x2_0)

trace_updated = []
colors = colorscale_list('Reds', number_colors=int(1.5*len(x)), return_rgb_only=True)
legend_step = len(x)//5

for i, updated_pos in enumerate(x):
    trace_updated.append(go.Scatter(
        x = [x[i][0]],
        y = [x[i][1]],
        mode = 'markers',
        hovertext = 'time: {}'.format(formatted(i*delta_t)),
        hoveron = "points",
        marker = dict(
            size = min(7+i//legend_step, 18),
            color = colors[i],
            line = dict(
                width = 2,
                color = 'Maroon'
            )
        ),
        name = 'At time {}'.format(formatted(i*delta_t)),
        showlegend = (i % legend_step == 0)
    ))

trace_identity = go.Scatter(
    x = xs,
    y = xs,
    mode = 'line',
    line = dict(
        width = 2,
        color = 'black'
    ),
    opacity = 0.5,
    name = 'Identity line'
)

fig1['layout']['title'] = '$\\text{Simulation of the system with initial condition } (x_1(0), x_2(0)) = '+\
'({}, {})'.format(x1_0, x2_0)+' \\text{ (arbitrary units)}$'

fig = dict(data=flow_nullclines + trace_updated + [trace_identity], layout=fig1['layout'])
plotly.offline.iplot(fig)

(c) So far we have treated each neuron separately, i.e. each neuron got its own equation. When we use matrix-vector notation, we can summarize these two equations into one:

$$\dot{\textbf{x}}(t) = - \textbf{x}(t) + f(W \textbf{x}(t) + \textbf{I})$$

where

  • $$\textbf{x}(t) = \begin{pmatrix} x_1(t)\\ x_2(t)\end{pmatrix}$$

  • $$W = \begin{pmatrix} 0 & w\\ w & 0\end{pmatrix}$$

  • $$I = \begin{pmatrix} I_1\\ I_2\end{pmatrix}$$

We here follow the convention in the literature by which vectors are denoted by bold-faced letters, and matrices by capitalized letters. Reprogram your program from (b) using Scipy/MATLAB/Numpy's intrinsic ability to compute with matrices and vectors! (Note that $\textbf{x}$ and $\textbf{I}$ are column-vectors!)

In [11]:
W = np.array([[0, w], [w, 0]])
I_bold = np.array([I, I])
In [12]:
def evolution_matrix(x_0, delta_t=delta_t, T=T):
    x = np.zeros((len(np.linspace(0, T, int(T//delta_t))), 2))
    x[0] = x_0
    
    for i in range(1, x.shape[0]):
        x[i] = (1-delta_t)*x[i-1] + delta_t*f(W.dot(x[i-1]) + I_bold)
    return x
In [14]:
x1_0, x2_0 = 1, 2
x = evolution_matrix([x1_0, x2_0]) 

trace_updated = []
colors = colorscale_list('Reds', number_colors=int(1.5*len(x)), return_rgb_only=True)
legend_step = len(x)//5

for i, updated_pos in enumerate(x):
    trace_updated.append(go.Scatter(
        x = [x[i][0]],
        y = [x[i][1]],
        mode = 'markers',
        hovertext = 'time: {}'.format(formatted(i*delta_t)),
        hoveron = "points",
        marker = dict(
            size = min(7+i//legend_step, 18),
            color = colors[i],
            line = dict(
                width = 2,
                color = 'Maroon'
            )
        ),
        name = 'At time {}'.format(formatted(i*delta_t)),
        showlegend = (i % legend_step == 0)
    ))

trace_identity = go.Scatter(
    x = xs,
    y = xs,
    mode = 'line',
    line = dict(
        width = 2,
        color = 'black'
    ),
    opacity = 0.5,
    name = 'Identity line'
)


fig1['layout']['title'] = '$\\text{Matrix-based simulation of the system with initial condition } (x_1(0), x_2(0)) = '+\
'({}, {})'.format(x1_0, x2_0)+' \\text{ (arbitrary units)}$'

fig = dict(data=flow_nullclines + trace_updated + [trace_identity], layout=fig1['layout'])
plotly.offline.iplot(fig)

(d) (Advanced) Plotting the vector field of derivatives in addition to the nullclines, use the Scipy/MATLAB/matplotlib function quiver to plot the vector field of derivatives!