AT2 – Neuromodeling: Problem set #4 NETWORKS

PROBLEM 1: Neuron with Autapse

Consider a neuron whose output feeds back onto itself via a synapse (such a synapse is often called an "autapse"). We assume that the neuron's firing rate $x$, i.e. its average number of action potentials per second, is given by the differential equation:

$$\dot{x}(t) = - x(t) + f\big(wx(t) + I\big)$$

where

  • $w = 0.04$ is the strength of the synaptic connection
  • and $I = -2$ is some constant, external (and inhibitory) background input

The input-output (or activation) function of the neuron is given by a sigmoidal function which we assume to be

$$f(s) ≝ 50 \, (1 + \tanh(s))$$

where $s$ is the total input to the neuron.

In [36]:
# /!\ 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

(a) Plot the neuron's activation function $f(s)$ against $s$. Choose a suitable range of inputs.

In [37]:
w = .04
I = -2

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

f = np.vectorize(f)
In [38]:
layout1 = go.Layout(
    title= '$\\text{Plot of the neuron\'s activation function } f$',
    hovermode= 'closest',
    xaxis= dict(
        title= '$\\text{Input } s$',
        ticklen= 5,
        zeroline= False,
    ),
    yaxis=dict(
        title= '$f(s)$',
        ticklen= 5
    ),
) 


# Plotting the evolution

xs = np.linspace(-3, 3, 50)

trace_f = go.Scatter(
    x = xs, 
    y = f(xs),
    mode = 'lines',
    line = dict(
        width = 2,
        dash = 'solid'
    ),
    hoverlabel = dict(
        namelength = -1
    ),
    name='f(s)'
)

plotly.offline.iplot(go.Figure(data=[trace_f], layout=layout1))

(b) Plot the derivatives $\dot{x}$ as a function of the neuron's firing rate $x$. What do the zero-crossings of this plot indicate?

In [39]:
def stable_unstable_alternating_array(length, stableFirst=True):
    arr = np.full((length,), False)
    arr[::2] = stableFirst
    return arr


decrease_color, increase_color = 'MediumTurquoise', 'Salmon'

def behavior_arrows(xs, stable=True, decrease_color=decrease_color, increase_color=increase_color,\
                    arrow_size=9, padding=1.2):
    
    annotations = []
    
    for x in xs:
        annotations.extend([
            dict(
                x = x-(padding if stable else arrow_size),
                y = 0,
                ax = x-(arrow_size if stable else padding),
                ay = 0,
                xref = 'x',
                axref = 'x',
                yref = 'y',
                ayref = 'y',
                showarrow = True,
                arrowhead = 2,
                arrowsize = 1,
                arrowwidth = 3,
                arrowcolor = increase_color if stable else decrease_color,
                opacity = 0.8
            ),
            dict(
                x = x+(padding if stable else arrow_size),
                y = 0,
                ax = x+(arrow_size if stable else padding),
                ay = 0,
                xref = 'x',
                axref = 'x',
                yref = 'y',
                ayref = 'y',
                showarrow = True,
                arrowhead = 2,
                arrowsize = 1,
                arrowwidth = 3,
                arrowcolor = decrease_color if stable else increase_color,
                opacity = 0.8
            )
        ])
    
    return annotations

color_fixed_points = 'DarkOliveGreen'
In [40]:
# Plotting the evolution

xs = np.linspace(0, 150, 301)

trace_der = go.Scatter(
    x = xs, 
    y = -xs + f(w*xs+I),
    mode = 'lines',
    line = dict(
        width = 2,
        dash = 'solid'
    ),
    hoverlabel = dict(
        namelength = -1
    ),
    showlegend=False
)

fixed_points_mask = np.abs(-xs + f(w*xs+I)) < .2
# print(-xs[fixed_points_mask] + f(w*xs[fixed_points_mask]+I))

stable_unstable_mask = stable_unstable_alternating_array(fixed_points_mask.sum())

trace_stable = go.Scatter(
    x = (xs[fixed_points_mask])[stable_unstable_mask], 
    y = [0]*len(stable_unstable_mask==True),
    mode = 'markers',
    marker = dict(
      color = color_fixed_points,
      size = 12
    ),
    name='stable fixed points'
)

trace_unstable = go.Scatter(
    x = (xs[fixed_points_mask])[~stable_unstable_mask], 
    y = [0]*len(stable_unstable_mask==False),
    mode = 'markers',
    marker = dict(
        size = 10,
        color = 'transparent',
        line = dict(
            color = color_fixed_points,
            width = 2
        )
    ),
    name='unstable fixed points'
)


layout3 = go.Layout(
    title = '$\\text{Plot of the derivative } \dot{x} \\text{ as a function of the neuron\'s firing rate } x '+\
    '\\text{ (arbitrary units)}$',
    hovermode= 'closest',
    xaxis= dict(
        title= '$\\text{Firing rate } x$',
    ),
    yaxis=dict(
        title= '$\\text{Derivative } \dot{x}$',
    ),
    annotations=behavior_arrows((xs[fixed_points_mask])[stable_unstable_mask])+\
    behavior_arrows((xs[fixed_points_mask])[~stable_unstable_mask], stable=False),
    legend=dict(
        orientation='h'
    )
) 


plotly.offline.iplot(go.Figure(data=[trace_der, trace_stable, trace_unstable], layout=layout3))

(c) Simulate the system for different initial conditions, $x(0) = 49, x(0) = 50$, and $x(0) = 51$, using a time step $dt = 0.1$. Chose a reasonable value of $T$. What are the outcomes of these simulations? Why?

We solve the differential equation numerically with resort to the Euler method:

$$ \begin{align*} \quad & \frac{Δx}{ Δt} = - x(t) + f\big(wx(t) + I\big)\\ ⟺ \quad & x(t+Δt) - x(t) = - x(t) Δt + f\big(wx(t) + I\big) Δt \\ ⟺ \quad & x(t+Δt) = (1 - Δt) x(t) + f\big(wx(t) + I\big) Δt \\ \end{align*} $$

In [88]:
delta_t = 0.1
T = 15
In [89]:
def Euler_x(x_0, delta_t=delta_t, T=T):
    
    x = np.zeros(len(np.linspace(0, T, int(T//delta_t))))
    x[0] = x_0
            
    for i in range(1, len(x)):
            x[i] = (1-delta_t)*x[i-1] + f(w*x[i-1] + I)*delta_t

    return x
In [90]:
layout_xs = go.Layout(
    title= 'Simulation of the system for different initial conditions (arbitrary units)',
    hovermode= 'closest',
    xaxis= dict(
        title= 'Time',
        ticklen= 5,
        zeroline= False,
        gridwidth= 2,
    ),
    yaxis=dict(
        title= '$\\text{Firing rate } x$',
        ticklen= 5,
        gridwidth= 2,
    ),
    showlegend= True
) 


legend_every = 1


# Plotting the evolution
traces_xs = []
values = [49, 50, 51]

for i, x_0 in enumerate(values):
    traces_xs.append(
        go.Scatter(
            x = np.linspace(0, T, int(T//delta_t)), 
            y = Euler_x(x_0),
            mode = 'lines',
            name = 'Initial condition: x(0) = {}'.format(x_0),
            line = dict(
                width = 2,
                dash = 'solid'
            ),
            hoverlabel = dict(
                namelength = -1
            ),
            showlegend = (i % legend_every == 0)
        )
    )

plotly.offline.iplot(go.Figure(data=traces_xs, layout=layout_xs))

(d) Redo the simulation of (c), but add a noise component to the system, so that:

$$\dot{x}(t) = - x(t) + f\big(wx(t) + I\big) + σ η(t)$$

where $η(t)$ is a Gaussian noise.

What happens for a noise value of $σ = 5$? What happens if you increase the noise (up to $σ = 80$)? Why?

We solve the stochastic differential equation numerically:

$$ \begin{align*} ⟺ \quad & x(t+Δt) - x(t) = - x(t) Δt + f\big(wx(t) + I\big) Δt + σ η(t) \sqrt{Δt} \\ ⟺ \quad & x(t+Δt) = (1 - Δt) x(t) + f\big(wx(t) + I\big) Δt + σ η(t) \sqrt{Δt} \end{align*} $$

In [44]:
delta_t = 0.1
T = 10
sigma = 5
In [45]:
def Euler_stochastic_x(x_0, delta_t=delta_t, T=T, sigma=sigma):
    
    x = np.zeros(len(np.linspace(0, T, int(T//delta_t))))
    x[0], sigma_sqrt_t = x_0, sigma*np.sqrt(delta_t)
            
    for i in range(1, len(x)):
            x[i] = (1-delta_t)*x[i-1] + f(w*x[i-1] + I)*delta_t + sigma_sqrt_t*np.random.standard_normal()
    return x
In [92]:
for seed in [1, 11, 21]:
    np.random.seed(seed)

    layout_xs = go.Layout(
        title= 'Simulation of the stochastic system for different initial conditions, with σ = {}'.format(sigma),
        hovermode= 'closest',
        xaxis= dict(
            title= 'Time',
            ticklen= 5,
            zeroline= False,
            gridwidth= 2,
        ),
        yaxis=dict(
            title= '$\\text{Firing rate } x$',
            ticklen= 5,
            gridwidth= 2,
        ),
        showlegend= True
    ) 


    legend_every = 1


    # Plotting the evolution
    traces_xs = []
    values = [49, 50, 51]

    for i, x_0 in enumerate(values):
        traces_xs.append(
            go.Scatter(
                x = np.linspace(0, T, int(T//delta_t)), 
                y = Euler_stochastic_x(x_0),
                mode = 'lines',
                name = 'Initial condition: x(0) = {}'.format(x_0),
                line = dict(
                    width = 2,
                    dash = 'solid'
                ),
                hoverlabel = dict(
                    namelength = -1
                ),
                showlegend = (i % legend_every == 0)
            )
        )


    plotly.offline.iplot(go.Figure(data=traces_xs, layout=layout_xs))
In [93]:
for seed in [2, 12, 22]:
    np.random.seed(seed)

    layout_xs_stoch = go.Layout(
        title= 'Simulation of the system for different initial conditions and noise magnitudes (arbitrary units)',
        hovermode= 'closest',
        xaxis= dict(
            title= 'Time',
            ticklen= 5,
            zeroline= False,
            gridwidth= 2,
        ),
        yaxis=dict(
            title= '$\\text{Firing rate } x$',
            ticklen= 5,
            gridwidth= 2,
        ),
        showlegend= True
    ) 


    legend_every = 1

    # Plotting the evolution
    traces_xs_stoch = []

    sigma_values = np.arange(1, 82, 20)

    colors = []
    colors.extend([
        colorscale_list('Blues', len(sigma_values)+3, return_rgb_only=True),
        colorscale_list('Oranges', len(sigma_values)+3, return_rgb_only=True),
        colorscale_list('Greens', len(sigma_values)+3, return_rgb_only=True)])

    x_values = [49, 50, 51]

    for j, x_0 in enumerate(x_values):
        for i, sig in enumerate(sigma_values):
            traces_xs_stoch.append(
                go.Scatter(
                    x = np.linspace(0, T, int(T//delta_t)), 
                    y = Euler_stochastic_x(x_0),
                    mode = 'lines',
                    name = 'x(0) = {}, σ = {}'.format(x_0, sig),
                    line = dict(
                        width = 2,
                        dash = 'solid',
                        color = colors[j][i+2],
                    ),
                    hoverlabel = dict(
                        namelength = -1
                    ),
                    showlegend = (i % legend_every == 0)
                )
            )


    plotly.offline.iplot(go.Figure(data=traces_xs_stoch, layout=layout_xs_stoch))
In [94]:
for seed in [3, 13, 23]:
    np.random.seed(seed)

    layout_xs_stoch = go.Layout(
        title= 'Simulation of the system for different initial conditions and noise magnitudes (arbitrary units)',
        hovermode= 'closest',
        xaxis= dict(
            title= 'Time',
            ticklen= 5,
            zeroline= False,
            gridwidth= 2,
        ),
        yaxis=dict(
            title= '$\\text{Firing rate } x$',
            ticklen= 5,
            gridwidth= 2,
        ),
        showlegend= True
    ) 


    legend_every = 1

    # Plotting the evolution
    traces_xs_stoch = []

    sigma_values = np.arange(1, 82, 20)

    colors = []
    colors.extend([
        colorscale_list('Blues', len(sigma_values)+3, return_rgb_only=True),
        colorscale_list('Reds', len(sigma_values)+3, return_rgb_only=True)])

    x_values = [49, 51]

    for j, x_0 in enumerate(x_values):
        for i, sig in enumerate(sigma_values):
            traces_xs_stoch.append(
                go.Scatter(
                    x = np.linspace(0, T, int(T//delta_t)), 
                    y = Euler_stochastic_x(x_0),
                    mode = 'lines',
                    name = 'x(0) = {}, σ = {}'.format(x_0, sig),
                    line = dict(
                        width = 2,
                        dash = 'solid',
                        color = colors[j][i+2],
                    ),
                    hoverlabel = dict(
                        namelength = -1
                    ),
                    showlegend = (i % legend_every == 0)
                )
            )


    plotly.offline.iplot(go.Figure(data=traces_xs_stoch, layout=layout_xs_stoch))
In [79]:
np.random.seed(4)

legend_every = 1

# Plotting the evolution

traces_xdot_stoch, trace_fixed_pts = [], []

sigma_values = [0, 5, 10, 80]

colors_curve = colorscale_list('Purples', len(sigma_values)+3, return_rgb_only=True)
colors_fixed_pts = colorscale_list('Greens', len(sigma_values)+3, return_rgb_only=True)

xs = np.linspace(0, 120, 2001)

deterministic_values = -xs + f(w*xs+I)

for i, sig in enumerate(sigma_values):
    
    derivative_values = deterministic_values + sig*np.random.standard_normal(len(xs))
    
    traces_xdot_stoch.append(
        go.Scatter(
            x = xs, 
            y = derivative_values,
            mode = 'lines',
            name = 'Noise magnitude: σ = {}'.format(sig),
            line = dict(
                width = 2,
                dash = 'solid',
                color = colors_curve[i+2],
            ),
            hoverlabel = dict(
                namelength = -1
            ),
            showlegend = (i % legend_every == 0),
            xaxis='x{}'.format(i+3 if i <= 1 else i-1), 
            yaxis='y{}'.format(i+3 if i <= 1 else i-1)
        )
    )
    
    fixed_points_mask = np.abs(derivative_values) < .11

    trace_fixed_pts.append(
        go.Scatter(
            x = xs[fixed_points_mask], 
            y = [0]*len(fixed_points_mask==True),
            mode = 'markers',
            marker = dict(
              color = colors_fixed_pts[i+2],
              size = 12
            ),
            name='fixed points for σ = {}'.format(sig),
            xaxis='x{}'.format(i+3 if i <= 1 else i-1), 
            yaxis='y{}'.format(i+3 if i <= 1 else i-1)
        )
    )


layout4 = go.Layout(
    title = '$\\text{Plot of the stochastic derivative } \dot{x} \\text{ as a function of the neuron\'s firing rate } x '+\
    '\\text{ for different noise magnitudes (arbitrary units)}$',
    xaxis=dict(
        title= '$x$',
        domain=[0, 0.45],
        anchor='y1'
    ),
    yaxis=dict(
        title= '$\\dot{x}$',
        domain=[0, 0.45],
        anchor='x1'
    ),
    xaxis2=dict(
        title= '$x$',
        domain=[0.55, 1]
    ),
    xaxis3=dict(
        title= '$x$',
        domain=[0, 0.45],
        anchor='y3'
    ),
    xaxis4=dict(
        title= '$x$',
        domain=[0.55, 1],
        anchor='y4'
    ),
    yaxis2=dict(
        title= '$\\dot{x}$',
        domain=[0, 0.45],
        anchor='x2'
    ),
    yaxis3=dict(
        title= '$\\dot{x}$',
        domain=[0.55, 1]
    ),
    yaxis4=dict(
        title= '$\\dot{x}$',
        domain=[0.55, 1],
        anchor='x4'
    )
) 

plotly.offline.iplot(go.Figure(data=traces_xdot_stoch+trace_fixed_pts, layout=layout4))
In [95]:
for seed in [6, 16, 26]:
    np.random.seed(seed)

    delta_values = 5
    values = list(np.arange(0, 81, delta_values))

    colors = colorscale_list('Greens', len(values)+3, return_rgb_only=True)

    xs = np.linspace(0, 120, 3001)
    deterministic_values = -xs + f(w*xs+I)

    traces_fp = []
    backgrounds = []

    for i, sig in enumerate(values):

        derivative_values = deterministic_values + sig*np.random.standard_normal(len(xs))
        fixed_points = (np.abs(derivative_values) < .07).astype(float)

        traces_fp.append(
            go.Scatter(
                x = xs, 
                y = (sig+5)*fixed_points-5,
                mode = 'markers',
                marker = dict(
                    symbol = 'square',
                    color = colors[i+2]
                ),
                name = 'σ = {}'.format(sig),
                showlegend = False
            )
        )

        if i%2 == 1:
            backgrounds.append(
                dict(
                    fillcolor='rgb(230, 230, 230)',
                    line=dict(
                        width=0
                    ),
                    opacity=0.45,
                    type='rect',
                    x0=xs[0]-10,
                    x1=xs[-1]+10,
                    y0=sig-delta_values/2.5,
                    y1=sig+delta_values/2.5,
                    layer='below'
                ))

    layout5 = go.Layout(
        title= '$\\text{Evolution of fixed points when there is a noise component for different noise magnitudes } σ '+\
        ' \\text{ (arbitrary units)}$',
        hovermode= 'closest',
        xaxis= dict(
            title= 'Fixed point values',
            ticklen= 5,
            zeroline= False,
            gridwidth= 2,
        ),
        yaxis=dict(
            title= 'White noise magnitude σ',
            range=[values[0]-1.5, values[-1]+0.5],
            tickvals=values,
            zeroline= False
        ),
        legend=dict(x=-.1, y=-0.2),
        shapes=backgrounds
    ) 

    plotly.offline.iplot(go.Figure(data=traces_fp, layout=layout5))