(Wong & Wang, 2006) Decision making rate model

import brainpy as bp

bp.math.set_dt(dt=0.0005)

A rate model for decision-making

Consider two excitatory neural assemblies, populations \(1\) and \(2\), that compete with each other through a shared pool of inhibitory neurons.

Let \(r_1\) and \(r_2\) be their respective population-firing rates, and the total synaptic input current \(I_i\) and the resulting firing rate \(r_i\) of the neural population \(i\) obey the following input-output relationship (\(F - I\) curve):

\[r_i = F(I_i) = \frac{aI_i - b}{1-\exp(-d(aI_i-b))}\]

which captures the current-frequency function of a leaky integrate-and-fire neuron. The parameter values are \(a\) = 270 Hz/nA, \(b\) = 108 Hz, \(d\) = 0.154 sec.

Assume that the ‘synaptic drive variables’ \(S_1\) and \(S_2\) obey

\begin{align} \frac{dS_1}{dt} &= F(I_1)\,\gamma(1-S_1)-S_1/\tau_s\\ \frac{dS_2}{dt} &= F(I_2)\,\gamma(1-S_2)-S_2/\tau_s \end{align}

where \(\gamma\) = 0.641. The net current into each population is given by

\begin{align} I_1 &= J_E S_1 + J_I S_2 + I_{b1} + J_{ext}\mu_1\\ I_2 &= J_E S_2 + J_I S_1 +I_{b2} +J_{ext}\mu_2. \end{align}

The synaptic time constant is \(\tau_s\) = 100 ms (NMDA time consant). The synaptic coupling strengths are \(J_E\) = 0.2609 nA, \(J_I\) = -0.0497 nA, and \(J_{ext}\) = 0.00052 nA. Stimulus-selective inputs to populations 1 and 2 are governed by unitless parameters \(\mu_1\) and \(\mu_2\), respectively. \(I_b\) is the background input which has a mean \(I_0\) and a noise component described by an Ornstein-Uhlenbeck process:

\begin{align} dI_{b1} &= (I_0 - I_{b1}) \frac{dt}{\tau_0} + \sigma dW \\ dI_{b2} &= (I_0 - I_{b2}) \frac{dt}{\tau_0} + \sigma dW \\ \end{align}

where \(I_0=0.3255\) nA, filter time constant \(\tau_0=2\) ms, and noise amplitude \(\sigma=0.02\) nA. \(dW\) is a Wiener process and note that when numerially integrating that with step size \(\frac{dt}{\tau_0}\) then \(\Delta W \sim \mathcal{N}(0, \frac{dt}{\tau_0})\), a normal distribution with mean 0 and variance \(\frac{dt}{\tau_0}\)

For the decision-making paradigm, the input rates \(\mu_1\) and \(\mu_2\) are determined by the stimulus coherence \(c'\) which ranges between 0 (0%) and 1 (100%):

\begin{align} \mu_1 &=\mu_0(1+c')\\ \mu_2 &=\mu_0(1-c') \end{align}

References:

• Wong K-F and Wang X-J (2006). A recurrent network mechanism for time integration in perceptual decisions. J. Neurosci 26, 1314-1328.

Implementation

gamma = 0.641  # Saturation factor for gating variable
tau = 0.06  # Synaptic time constant [sec]
tau0 = 0.002  # Noise time constant [sec]
a = 270.
b = 108.
d = 0.154

I0 = 0.3255  # background current [nA]
JE = 0.3725  # self-coupling strength [nA]
JI = -0.1137  # cross-coupling strength [nA]
JAext = 0.00117  # Stimulus input strength [nA]
sigma = 1.02  # nA


# I0 = 0.3255 # background current [nA]
# JE = 0.2609  # self-coupling strength [nA]
# JI = -0.0497  # cross-coupling strength [nA]
# JAext = 5.2e-4  # Stimulus input strength [nA]
# sigma = 0.02 # nA

class DecisionMaking(bp.NeuGroup):
  @staticmethod
  @bp.odeint
  def integral(s1, s2, t, mu0, coh, Ib1, Ib2):
    I1 = JE * s1 + JI * s2 + Ib1 + JAext * mu0 * (1. + coh)
    r1 = (a * I1 - b) / (1. - bp.math.exp(-d * (a * I1 - b)))
    ds1dt = - s1 / tau + (1. - s1) * gamma * r1

    I2 = JE * s2 + JI * s1 + Ib2 + JAext * mu0 * (1. - coh)
    r2 = (a * I2 - b) / (1. - bp.math.exp(-d * (a * I2 - b)))
    ds2dt = - s2 / tau + (1. - s2) * gamma * r2

    return ds1dt, ds2dt

  @staticmethod
  @bp.odeint
  def int_noise(Ib1, Ib2, t, sigma):
    dIb1 = (I0 - Ib1) / tau0 + sigma / tau0 * bp.math.random.randn(*Ib1.shape)
    dIb2 = (I0 - Ib2) / tau0 + sigma / tau0 * bp.math.random.randn(*Ib2.shape)
    return dIb1, dIb2

  def __init__(self, size, **kwargs):
    super(DecisionMaking, self).__init__(size=size, **kwargs)

    # parameters
    self.mu0 = 20.  # Stimulus firing rate [spikes/sec]
    self.coh = 0.2  # # Stimulus coherence [%]
    self.sigma = 0.02  # nA

    # variables
    self.s1 = bp.math.zeros(self.num)
    self.s2 = bp.math.zeros(self.num)
    self.Ib1 = bp.math.zeros(self.num)
    self.Ib2 = bp.math.zeros(self.num)

  def update(self, _t, _dt):
    self.s1, self.s2 = self.integral(self.s1, self.s2, _t,
                                     self.mu0, self.coh, self.Ib1, self.Ib2)
    self.Ib1, self.Ib2 = self.int_noise(self.Ib1, self.Ib2, _t, self.sigma)

Analysis

dm = DecisionMaking(1)

\(\mu_0 = 0, c'=0.\)

phase = bp.symbolic.PhasePlane(
  dm,
  target_vars=dict(s2=[0., 1.], s1=[0., 1.]),
  pars_update=dict(mu0=0., coh=0., Ib1=0.3297, Ib2=0.3297),
  numerical_resolution=0.0005,
  options={'escape_sympy_solver': True})

phase.plot_nullcline()
phase.plot_fixed_point()
_ = phase.plot_vector_field(show=True)

\(\mu_0 = 30, c'=0.\)

phase = bp.symbolic.PhasePlane(
  dm,
  target_vars=dict(s2=[0., 1.], s1=[0., 1.]),
  pars_update=dict(mu0=30., coh=0., Ib1=0.3297, Ib2=0.3297),
  numerical_resolution=0.0003,
  options={'escape_sympy_solver': True})

phase.plot_nullcline()
phase.plot_fixed_point()
_ = phase.plot_vector_field(show=True)

\(\mu_0 = 30, c'=0.5\)

phase = bp.symbolic.PhasePlane(
  dm,
  target_vars=dict(s2=[0., 1.], s1=[0., 1.]),
  pars_update=dict(mu0=30., coh=0.5, Ib1=0.3297, Ib2=0.3297),
  numerical_resolution=0.0003,
  options={'escape_sympy_solver': True})

phase.plot_nullcline()
phase.plot_fixed_point()
_ = phase.plot_vector_field(show=True)

\(\mu_0 = 30, c'=-0.5\)

phase = bp.symbolic.PhasePlane(
  dm,
  target_vars=dict(s2=[0., 1.], s1=[0., 1.]),
  pars_update=dict(mu0=30., coh=-0.5, Ib1=0.3297, Ib2=0.3297),
  numerical_resolution=0.0003,
  options={'escape_sympy_solver': True})

phase.plot_nullcline()
phase.plot_fixed_point()
_ = phase.plot_vector_field(show=True)

\(\mu_0 = 30, c'=1.\)

phase = bp.symbolic.PhasePlane(
  dm,
  target_vars=dict(s2=[0., 1.], s1=[0., 1.]),
  pars_update=dict(mu0=30., coh=1., Ib1=0.3297, Ib2=0.3297),
  numerical_resolution=0.0003,
  options={'escape_sympy_solver': True})

phase.plot_nullcline()
phase.plot_fixed_point()
_ = phase.plot_vector_field(show=True)