Supplemental Materials for ”Mean-field Approximations for Coupled
Populations of Generalized Linear Model Spiking Neurons with Markov
Refractoriness” by Taro Toyoizumi, Kamiar Rahnama Rad, and LiamPaninski, Neural Computation, Volume 21, Issue 5, pp. 1203-1243.
Color figures:
BA
54

1.5

I

h3
2

0.01

0.02

0.03

0.04

00.05

0.01

0.02

0.03

0.04

0.05

00.4
simulation
noH
exact
mf

40

0.3

E[h]

firing intensity [Hz]

1160

20
0020.5

1003

20.2
0.1

0.01

0.02

0.03

time [s]0.04

000.050.01

0.02

0.03

0.04

0.05

time [s]Figure 2: Direct solution of the master equation in the one-dimensional case. (A) Top: A single
neuron is simulated with a pulse input I(t) shown. Bottom: Comparison of the firing rates
computed via: direct Monte Carlo sampling (simulation); the inhomogeneous Poisson case, withno history dependence, i.e., w = 0 (noH); the master equation of equation. 3.2 (exact); and themean-field approximation methods in the limit of Nc → ∞ introduced in section 4. the mean-fieldapproximation gives good approximations of the output firing intensity. (B) Top: The probability
distribution P (h, t); color axis indicates the height of the probability density. The self-history
term h(t) was the convolution of the output spike train and a single exponential filter with time
constant τ = 10 ms and weight J = −1; thus the neuron was inhibited for a brief period following
each spike. Note the discontinuous behavior visible at stimulus onset, where probability mass
jumps significantly towards h = 1, then decays back down towards zero, according to equation
3.1. Bottom: The mean of h calculated from the probability distribution P (h, t). This E[h] also
increases following the input pulse, which has the effect of slowing down the firing rate: since theweight J is negative, the history term is inhibitory here.

1

A

B
55

3000

simulation
noH
mf

50

2500

4540

firing intensity [Hz]firing intensity [Hz]simulation
noH
mf

35
30
25
20
15

2000

1500

1000

500

10
50

0.01

0.02

0.03

0.04

000.05time [s]

0.01

0.02

0.03

0.04

0.05

time [s]Figure 3: Comparison of the direct Monte Carlo simulation of a single neuron (solid),
the firing intensity without recurrent input H (dotted), and the mean-field approximation(dashed). (A) Low step input (baseline I = 2, peak I = 4) is applied. The self-inhibition
term is set to J = −5 and τ = 2 ms. (B) High step input (baseline I = 2, peak I = 8)
is applied. The self inhibition terms is set to J = −5 and τ = 2 ms. Because of the large
non-Gaussian fluctuations of input caused by the strong self-inhibition, the mean-fieldapproximation loses accuracy here.

2

B
15

15101010

1000.05

time [s]0.05

0−0.05

00.05

10050

0

50−0.05

0.05

10050

0−0.050

ν 22 ( ·| 0)ν 21 ( ·| 0)1000−0.05

ν 21 ( ·| 0)0−0.0550

time [s]

00.05

00.05

100

50

0−0.050.05

50−0.05ν 22 ( ·| 0)5

ν 12 ( ·| 0)15
ν 11 ( ·| 0)15
ν 12 ( ·| 0)ν 11 ( ·| 0)A0time [s]0.05

500−0.05time [s]Figure 4: Comparison of cross-correlation functions of two neurons calculated by direct
simulation (solid), calculated from the mean-field approximation (see equations 4.9 and4.10; dashed), and the mean-field approximation with finite size correction (see equations 4.11 and 4.12; dot-dashed). (A) Weak self-interaction: parameters are set toJ11 = −1, J22 = −1, J21 = 1, J12 = −0.5, τ11 = 10ms, τ22 = 10ms, τ21 = 10ms, τ12 = 20ms,and I1 = 2, I2 = 4. The result with finite size correction is obtained after 5 iterationsof the self-consistent equations. The calculation in the mean-field limit captures only
the steady state firing intensity but not any temporal cross-correlation function. Crosscorrelation functions are well approximated with the finite size correction terms. (B)Strong self-interaction: parameters are set to J11 = −5, J22 = −5, J21 = 1, J12 = −0.5,
τ11 = 2ms, τ22 = 2ms, τ21 = 10ms, τ12 = 20ms, and I1 = 2, I2 = 4. Because of the large
non-Gaussian fluctuation of the input caused by the strong self-inhibition term, the meanfield approximations does not work very well. In this case the finite size solution converges
rather slowly. Hence, the result is after 10 iterations of the self-consistent equations.

3AB
3C
3

3

20

1018
2161514I (0)

120

108

1

1I (0)

120

10I (0)

26−1

4−220

8−1

6
4−2−1
5−2

2−3−3−2

J−1−3−30−2−1

J−3−30−2−1J0Figure 5: Comparison of (A) a simulated firing rates of a single neuron and (B) theself-consistent solutions of equation 4.13 for different input levels, I (0) , and self-inhibition
strengths, J. (C) The percentage errors of the mean-field approximation are shown withrespect to the simulated firing intensity. We set τ = 10 ms.
AB501−2230
20

01

102

00.6

00.05

0.1

00.02

0.04

00.4
0G 12 [Hz]0.2
0−2−2znorm. filter

0.8

40z

simulation (I0=2)
theory (I0=2)
simulation (I0=4)
theory (I0=4)

G 11 [Hz]1.2

1005−42−0.2
−0.01

00.01

0.02

0.03

0.04

0.05

−600.05

time [s]time [s]

0.1

00

0.02

0.04

time [s]Figure 6: Comparison of the analytically derived and simulated input-output filters. (A)A simulated normalized linear input-output filter for I (0) = 2 (thin-dotted) and I (0) = 4(thick-dotted) are compared with analytical results at I (0) = 2 (thin-solid) and I (0) =4 (thick-solid) of equation 4.18. Because the strength of the self-inhibition component
changes with the baseline input through f ′ (I (0) + µ(0) ), the input-filter is sharpened forhigh baseline input. Other parameters are set to J = 1, τ = 10 ms, K(t) = e−t/τI Θ(t)with τI = 20 ms and σξ = 1 Hz. All filter amplitudes are normalized by their maximum
amplitudes. (B) Left panels show the comparison of analytically derived input-output
filters G11 (t) and G12 (t) of equation 4.20 (solid lines) and simulated input-output filters(dotted lines) for I (0) = 2 (thin lines) and I (0) = 4 (thick lines). Right panels show neuron
1’s spatio-temporal input-output filter G̃1 (t, z) for I (0) = 2 (top) and I (0) = 4 (bottom).
See text for parameter values.

4A

B
55

3000
simulation
noH
mf

50

firing intensity [Hz]45

firing intensity [Hz]simulation
noH
mf

2500

40
35
30
25
20

2000

1500

1000

15

500

10
5
0

0.01

0.02

0.03

0.04

000.05time [s]

C

0.02

0.03

time [s]

0.04

0.05

D
140

150
simulation
noH
mf

120

simulation
noH
mf

100

firing intensity [Hz]firing intensity [Hz]

0.01

8060
40

10050

20
000.01

0.02

0.03

0.04

000.05

0.01

0.02

0.03

0.04

0.05

time [s]time [s]Figure 7: Comparison of the mean firing intensity obtained from direct simulation of a single
neuron (solid), without recurrent input H and no Markov refractory effect (dotted), and themean-field equation in the limit of Nc → ∞ (dashed). In all the cases, the markov refractory
effect with τr = 2ms/(M −1) was included. (A) Weak step input (baseline I = 2, peak I = 4) andrefractoriness M = 3; J = −1, τ = 10 ms. Although the model includes strong Markov refractory
effect, the mean-field approximation is not abolished. (B) Strong step input (baseline I = 2,peak I = 8) and refractoriness M = 10; J = −1, τ = 10ms. The mean-field approximation also
approximates very well the transient oscillation caused by the the strong refractory effect and steep
current. (C) The mean-field approximation works for colored noise input: τI dI(t)dt = −I(t) + ξ(t)′′with E[ξ(t)] = 4, Cov[ξ(t), ξ(t )] = 0.1δ(t − t ), and τI = 20 ms; M = 3, J = −1, τ = 10ms. The
initial error is due to the preset initial values deviate from the true ones. (D) Sinusoidal input
input: I(t) = 2 + sin(2π · 10Hz · t); M = 3, J = −1, τ = 10ms. The mean-field results give good
approximations.

5B

10

100

0.05

0−0.0500time [s]

0.05

50

0−0.0550−0.05

0.05

100

50

151050−0.05

ν 22 ( ·| 0)ν 21 ( ·| 0)100

1500.05

time [s]0−0.0515

1550−0.05

0.05

520

1001020
ν 21 ( ·| 0)0−0.0520

ν 22 ( ·| 0)520ν 12 ( ·| 0)15

ν 11 ( ·| 0)15
ν 12 ( ·| 0)ν 11 ( ·| 0)A0time [s]

0.05

00.05

00.05

10
50−0.05

time [s]Figure 8: Comparison of spike cross-correlation functions of two neurons obtained by direct
numerical simulation (solid), mean-field approximation in the limit of Nc → ∞ (dashed),and with finite size correction terms (dot-dashed). While the approximation in the limitof Nc → ∞ only captures the strong Markov refractory effect, the weaker refractory effectcaused by the self-interaction terms is well captured with the finite size correction. The finite size correction provides a good approximation despite the strong refractoriness becauseof an explicit evaluation of the Markov refractory effect. Mean-field equations are iteratively solved over 5 iterations. (A) Two neurons are directly coupled: parameters are set toJ11 = −1, J22 = −1, J12 = −0.5, J21 = 1, τ11 = 10ms, τ22 = 10ms, τ12 = 20ms, τ21 = 10ms,
M = 3, I1 = 2, I2 = 4. (B) Two neurons are not directly connected but receive common
input from a third neuron: parameters are set to J13 = J23 = 2, τ13 = τ23 = 10ms, M = 3,
I1 = I2 = I3 = 2.

6