Example Data Project: the Train Illusion

Example Data Project: the Train Illusion#

By Neuromatch Academy

Content creators: Marius ‘t Hart, Megan Peters, Paul Schrater, Jean Laurens, Gunnar Blohm

Production editor: Spiros Chavlis

Disclaimer: this is a “toy” data neuroscience project used to demonstrate the 10 step procedure of how-to-model. It is not meant to be state of the art research.

Setup#

Install dependencies#

# Imports

# for matrices and plotting
import numpy as np
import matplotlib.pyplot as plt

# for random distributions
from scipy.stats import norm, poisson

# for logistic regression
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

Plot Functions#

Generate Data#

Show code cell source Hide code cell source

# @title Generate Data

def generateSpikeTrains(seed=37):

  gain  = 2
  neurons  = 50
  movements  = [0, 1, 2]
  repetitions = 800

  np.random.seed(seed)

  # set up the basic parameters:
  dt = 1/100
  start, stop = -1.5, 1.5
  t = np.arange(start, stop+dt, dt)  # a time interval
  Velocity_sigma = 0.5  # std dev of the velocity profile
  Velocity_Profile = norm.pdf(t ,0, Velocity_sigma)/norm.pdf(0, 0, Velocity_sigma)  # The Gaussian velocity profile, normalized to a peak of 1

  # set up the neuron properties:
  Gains = np.random.rand(neurons) * gain  # random sensitivity between 0 and `gain`
  FRs   = (np.random.rand(neurons) * 60 ) - 10  # random base firing rate between -10 and 50

  # output matrix will have this shape:
  target_shape = [len(movements), repetitions, neurons, len(Velocity_Profile)]

  # build matrix for spikes, first, they depend on the velocity profile:
  Spikes = np.repeat(Velocity_Profile.reshape([1, 1, 1, len(Velocity_Profile)]),
                     len(movements)*repetitions*neurons,axis=2).reshape(target_shape)

  # multiplied by gains:
  S_gains = np.repeat(np.repeat(Gains.reshape([1, 1, neurons]),
                                len(movements)*repetitions, axis=1).reshape(target_shape[:3]),
                      len(Velocity_Profile)).reshape(target_shape)
  Spikes = Spikes * S_gains

  # and multiplied by the movement:
  S_moves = np.repeat(np.array(movements).reshape([len(movements), 1, 1, 1]),
                      repetitions*neurons*len(Velocity_Profile),
                      axis=3).reshape(target_shape)
  Spikes = Spikes * S_moves

  # on top of a baseline firing rate:
  S_FR = np.repeat(np.repeat(FRs.reshape([1, 1, neurons]),
                             len(movements)*repetitions, axis=1).reshape(target_shape[:3]),
                   len(Velocity_Profile)).reshape(target_shape)
  Spikes = Spikes + S_FR

  # can not run the poisson random number generator on input lower than 0:
  Spikes = np.where(Spikes < 0, 0, Spikes)

  # so far, these were expected firing rates per second, correct for dt:
  Spikes = poisson.rvs(Spikes * dt)

  return Spikes


def subsetPerception(spikes, seed=0):

  movements = [0, 1, 2]
  split = 400
  subset = 40
  hwin = 3

  [num_movements, repetitions, neurons, timepoints] = np.shape(spikes)

  decision = np.zeros([num_movements, repetitions])

  # ground truth for logistic regression:
  y_train = np.repeat([0, 1, 1], split)
  y_test  = np.repeat([0, 1, 1], repetitions - split)

  m_train = np.repeat(movements, split)
  m_test  = np.repeat(movements, split)

  # reproduce the time points:
  dt = 1/100
  start, stop = -1.5, 1.5
  t = np.arange(start, stop+dt, dt)

  w_idx = list((abs(t) < (hwin*dt)).nonzero()[0])
  w_0 = min(w_idx)
  w_1 = max(w_idx) + 1  # python...

  # get the total spike counts from stationary and movement trials:
  spikes_stat = np.sum(spikes[0, :, :, :], axis=2)
  spikes_move = np.sum(spikes[1:, :, :, :], axis=3)

  train_spikes_stat = spikes_stat[:split, :]
  train_spikes_move = spikes_move[:, :split, :].reshape([-1, neurons])

  test_spikes_stat = spikes_stat[split:, :]
  test_spikes_move = spikes_move[:, split:, :].reshape([-1, neurons])

  # data to use to predict y:
  x_train = np.concatenate((train_spikes_stat, train_spikes_move))
  x_test = np.concatenate((test_spikes_stat, test_spikes_move))

  # this line creates a logistics regression model object, and immediately fits it:
  population_model = LogisticRegression(solver='liblinear',
                                        random_state=seed).fit(x_train, y_train)

  # solver, one of: 'liblinear', 'newton-cg', 'lbfgs', 'sag', and 'saga'
  # some of those require certain other options
  #print(population_model.coef_)       # slope
  #print(population_model.intercept_)  # intercept

  ground_truth = np.array(population_model.predict(x_test))
  ground_truth = ground_truth.reshape([3, -1])

  output = {}
  output['perception'] = ground_truth
  output['spikes'] = spikes[:, split:, :subset, :]

  return output


def getData():

  spikes = generateSpikeTrains()
  dataset = subsetPerception(spikes=spikes)

  return dataset


dataset = getData()
perception = dataset['perception']
spikes = dataset['spikes']

Phenomenon#

Part of Steps 1-2

The train illusion occurs when sitting on a train and viewing another train outside the window. Suddenly, the other train seems to move, i.e. you experience visual motion of the other train relative to your train. But which train is actually moving?

Often people mix this up. In particular, they think their own train might be moving when it’s the other train that moves; or vice versa. The illusion is usually resolved once you gain vision of the surroundings that lets you disambiguate the relative motion; or if you experience strong vibrations indicating that it is indeed your own train that is in motion.

We would like to understand more about why this illusion occurs. We know from the literature that in the brain, the inner ear organs (vestibular organs) sense self-motion acceleration. However, these acceleration signals have low signal-to-noise ratio, especially for small accelerations (like in a train). Thus, judging self-motion from vestibular signals can be tricly, which is believed to result in the train illusion.