A complete mean-field theory for dynamics of binary recurrent neural networks

0.47 MB·

by  Farzad Farkhooi

A complete mean-field theory for dynamics of
binary recurrent neural networks
7
1
0
Farzad Farkhooia,b,1 and Wilhelm Stannata,b
2

n
aInstitut fu¨r Mathematik, Technische Universit¨at Berlin, D-10623
a
J Berlin, Germany.

5 bBernstein Center for Computational Neuroscience, D-10115
2
Berlin, Germany.

] [email protected]
h
p
-
o January 26, 2017
i
b
.
s
c Abstract
i
s Recurrent networks of binary units provide a prototype for complex
y interactionsinneuronalnetworks. Inthetheoreticalanalysisoftheircol-
h lectivedynamics,itisoftenassumedthatunitinteractionsinthenetwork
p
are homogeneous. Here, we derive a unified theory that encompasses
[
the macroscopic dynamics of recurrent networks with arbitrary connec-

1 tivity architectures. We identify conditions to the connectivity patterns
v for which the population averages, according to mean-field theory, con-
8 verge to a deterministic limit as the size of the system increases. The
2 mathematical analysis, using the martingale structure of the network’s
1 Markovian dynamics, further quantifies the aggregate fluctuation of the
7
population activity in a finite-size system. Our theory allows, without
0
additionaleffort,theinvestigationofsystemsthatdonotexhibitapoint-
.
1 wiseconvergencetothemean-fieldlimit. Inparticular,weuncoveranovel
0 dynamicstateinwhichasystemwithinhomogeneousmicro-structuresin
7 the network connectivity, along with an asynchronous behavior, also ex-
1
hibitsstochasticsynchronization. Thisphenomenonresemblesnon-trivial
:
v fluctuationsintheaveragepopulationactivity,whichcansurviveirrespec-
i tive of system size. The concise mathematical analysis that is presented
X
hereenablesformulationofmean-fieldequationforthesecasesintheform
r ofastochasticordinarydifferentialequation. Ourresultsindicatethees-
a
sential elements for implementation of connectivity motifs in the general
theory of recurrent networks.
1
1 Introduction
The fundamental characteristics of collective spiking dynamics in neural sys-
tems remain elusive. For instance, it is a challenge to consistently explain
experimental data that exhibit wide modes of complex macroscopic dynam-
ics, ranging from a variety of synchronous to asynchronous states, in cortical
circuitry [13, 14, 16, 10, 11]. The prevailing theoretical models, so-called bal-
anced networks, exhibit self-generated fluctuations [18, 12, 6] that successfully
explain the irregular asynchronous state that has been observed in the cortical
spiking activity. The description of the asynchronous state is closely related
to the formulation of a mean-field theory (MFT) that replaces the statistics of
individual units with population averages; this is similar to the statistical me-
chanics of spin glasses [5] in the thermodynamic limit. In the analysis of the
balanced state, the number of units in the network, N, is assumed to be large
(i.e. N → ∞) and the average number of direct interactions per unit, K, is
assumed to be substantial and homogeneous (i.e. K →∞). To date, no theory
has been developed for finite-size corrections with respect to realistic values of
N and K. Additionally, a growing body of experimental evidence is emerging
[15, 7, 9] that contradicts the assumptions that the network architecture fol-
lows a homogeneous and uniformly random Bernoulli directed graph (directed
Erd˝os-R´enyinetwork)thatiscommonlyusedintheclassicalanalysisofcortical
models[18,12,6]. Itisnotyetunderstoodhowthecombinationoftheserealistic
connectivitymicro-structuresandthefinite-sizeofbiologicalnetworksinfluence
the modes of the network dynamics. Therefore, it is important to determine
whether the existence of these inhomogeneous structures in the architecture of
connectivity is able to alter the MFT predictions for cortical networks behav-
ior. Understanding the various spiking dynamic states that can emerge within
these realistic cortical networks has important implications for elucidating the
temporal components of neuronal coding. Here, we develop a mathematical
theory to determine the collective dynamics of a finite-sized neuronal popula-
tion. We first develop a general MFT for binary units that interact within an
arbitrary network architecture. To simplify our computations, we demonstrate
our general method for connectivity matrices that have a fixed row sum, i.e.,
a fixed number K of inputs to units. Interestingly, we identify the conditions
on the connectivity structure that are necessary to guarantee the convergence
of the average population activity to a deterministic limit as N → ∞. We
show that whenever these conditions are satisfied, the classical results can be
recovered in the thermodynamic limit. Using the martingale structure of the
Markoviandynamicsofthenetwork,wecomputethefinitesizeeffectandderive
acorrespondingdynamicalequationintermsofanOrnstein-Uhlenbeckprocess.
Furthermore, our approach reveals a novel dynamical state in a network with
inhomogeneous coupling, in which the fluctuations of the average population
activity survive irrespective of the network size. In this polychronous state [8],
precise time-locked population events are able to emerge, whereas the overall
network activity remains asynchronous.
2
2 Mean-field theory (MFT) for binary recurrent
network
Consider a network that is described by an adjacency matrix J=(J ) of N bi-
ij
naryunits, whosetheircurrent statesaredenotedas n(t):=(n (t),...,n (t)).
1 N
The vector n(t) is a time-continuous Markov chain on {0,1}N with Q-matrix,
where Q(n,m)=0 if and only if ||n−m||≥2 and
(cid:40)
f (n) if m−n=e
Q(n,m)= i i
1−f (n) if m−n=−e
i i
where e(i) = δ denotes the i-th unit vector. In order to comply with the
j ij
centralization property of Q-matrices, it follows that:
N
(cid:88) (cid:0) (cid:1)
Q(n,n)=− n (t) 1−2f (n(t)) +f (n(t)).
i i i
i=1
The analytical gain function f (n(t)) defines the state transition rate of a unit
i
i, given the state of the network, n(t), and it is assumed to take values in the
range [0,1]. Typically, this function is written as f (u (t)), where u (t), which
i i i
represents the input to the unit i with the scaling parameter 0 < γ ≤ 1, is
written as
N
u (t):=J¯K−γ(cid:88)J n (t)+Kγµ (1)
i i ij j i 0,i
j=1
where J¯, K , and µ are the coupling strength, the number of input units
i 0,i
and the external drive to the i-th unit, respectively. For the convenience of
the current presentation, we consider here networks with K = K, f = f and
i i
µ =µ for all i. We will provide below (in eqns. 10 and 11) conditions on the
0,i 0
network structure, J, that imply the convergence of the averaged population
activity in the network towards a deterministic limit
N
! 1 (cid:88)
m(t)= lim n (t). (2)
N→∞N i
i=1
Here, m(t) is known as the mean-field limit and has the following temporal
dynamics:
d
m(t)=−m(t)+F(m(t)) (3)
dt
for some a priori unknown function F. In order to determine F, we use the
following semimartingale decomposition, that specifies the difference between
n¯(t) := 1 (cid:80)N n (t) (i.e., the average population activity of a finite-size net-
N i=1 i
3
work) and the mean-field limit m(t) of the system,
(cid:90) t
(cid:0) (cid:1)
n¯(t)−m(t)=(n¯(0)−m(0))− ds n¯(s)−m(s)
0
(cid:90) t (cid:0) 1 (cid:88)N (cid:1)
+ ds f(u (t))−F(m(s)
N i
0 i=1
+M(t), (4)
whereM(t)issomesquareintegrablemartingalethat, accordingtothegeneral
theory of Markov processes [17], satisfies
E(cid:2)M(t)2(cid:3)= 1 (cid:90) tds −Q(n(s),n(s))≤ t . (5)
N2 N
0
Note that E[M(t)]=0 and, in general, M(t) specify finite-size fluctuations in
the average population activity. Provided that m(t) exists (refer to eqns. 10
and 11 for a justification of this ansatz), we can construct the function F by
expanding 1 (cid:80)N f(u (t)) as N →∞ in eqn. 4 around
N i=1 i
µ (t):=K1−γJ¯m(t)+Kγµ . (6)
1 0
Using the lemma that is described later in the Materials and Methods section,
we obtain the following series expansion
F(m(t))=f(µ (t))+(cid:88)∞ f(r)(µ1(t))µ (t). (7)
1 r! r
r=2
where µ represents the average input to a unit in the network at time t. The
1
higherordercoefficientscanbecomputedbyexpandingµ :=lim 1 (cid:80)N [(u −
r N→∞ N i=1 i
µ )r]. The second order coefficient is given by:
1
µ (t)=J¯2K1−2γm(t)(1−m(t)) (8)
2
and the subsequent coefficients are given by:
r r−q
µ (t)=J¯rK−rγ(cid:88)a m(t)q(cid:88)b m(t)s (9)
r q s
q=0 s=0
where,
(cid:18) (cid:19)
r
a := (−1)qKq
q q
and
b :=S(r−q,s)(K)
s s
Here, S is a Stirling number of the second kind and (.) denotes the falling
s
factorial. In the binomial expansion of µ (t) given in eqn. 9, the summation
r
4
over j is performed using the ansatz that m(t) exists; thereafter, summation
overiintheaverageoperatorlim 1 (cid:80)N [.]isapplied. Inordertoprovide
N→∞ N i=1
the necessary conditions for the existence of a deterministic limit, m(t), the
summation order must be changed. Therefore, the first condition for m(t) and
µ (t) to exist is
1
N (cid:34) N (cid:35)
1 (cid:88) (cid:88)
lim J −K =0. (10)
N→∞N ij
j=1 i=1
Eqn. 10 implies that the coefficient in front of f(cid:48) in the series expansion that
leads to eqn. 7 vanishes in the thermodynamic limit (refer to the Materials
and Methods section for a detailed description). The second condition for the
pointwiseconvergenceofanaveragedpopulationactivitytothemean-fieldlimit
in eqn. 2 is given by
N (cid:34) N (cid:35)
1 (cid:88) (cid:88)
lim J J −K(K−1) =0. (11)
N→∞N ij1 ij2
j1(cid:54)=j2 i=1
Thisconditionspecifiesthat,asN →∞,themeancovarianceofcolumnsinthe
connectivity matrix J must decay to zero. The higher order condition can be
similarlydeterminedinordertoachieveapointwiseconvergenceoftheaveraged
population activity to its mean-field limit, as described in eqn. 3 (refer to the
Materials and Methods section for calculation details).
It is often of interest to study network dynamics when the number of inputs
into units is large (e.g. K →∞). In order to study this classical case, we must
investigate the asymptotic behavior of µ in eqn. 9 in the order of K; it can be
r
observed that the odd coefficients are given by
µ ∼O(K1−(2k+1)γ)
2k+1
and the even coefficients are given by
µ ∼O(K1−2kγ)+(2k)!!µk,
2k 2
for k ∈ N. Hence, it is apparent that the scaling parameter γ plays a critical
role in the large K limit. In classical balanced-state theory [18], the scaling
parameter γ is generally assumed to take the value 0.5; in this case, µ ∼O(1)
2
and the mean-field coefficients of eqn. 3 converge as K → ∞, towards the
central moments of a Gaussian distribution function. As a result, the related
power series that is given by eqn. 7 can be reformulated in terms of a simple
Gaussian integral so that, in this special case, eqn. 3 reduces to
d (cid:90) f(x) (x−µ (t))2
m(t)=−m(t)+ dx exp(− 1 ). (12)
(cid:112)
dt 2πµ2(t) 2µ2(t)
In the above analysis, we first take N →∞ to arrive at the mean-field of eqn.3
and then we consider K → ∞ in order to recover eqn. 12. This derivation
5
generalizes the proposed limit (K → ∞) in eqn. A.7 of van Vreeswijk and
Sompolinsky’s seminal work [18]. It is noteworthy that the analysis here shows
that the finite K correction to eqn. 12 is relatively small for a homogeneous
network; this is consistent with the approximation that is made in section 4.2
of [18].
Thesemimartingaledecompositionthatisgivenineqn.4providesadditional
informationontheexpansionofnetworksthathavefinitesizeN. Usingeqn. 5,
wecandeterminethefluctuationsmagnitudeoftheaveragepopulationactivity
in finite networks in the mean-square sense as
E(cid:0)M(t)2(cid:1)= 1 (cid:90) ds(cid:88)N (cid:0)n (s)(1−2f(u (s)))+f(u (s))(cid:1).
N2 i i i
i=1
and, by expanding 1 (cid:80)N f(u (t)) at µ (t), we arrive at
N i=1 i 1
E(cid:0)M(t)2(cid:1)= 1 (cid:90) tds(cid:0)m(s)(1−2(g(µ )+R))+g(µ )+R)(cid:1)
N 1 1
0
where g(µ ) := f(µ (t))+f(cid:48)(cid:48)(µ )µ /2 and R := (cid:80)∞ f(r)(µ1)µ denotes the
1 1 1 2 r=3 r! r
remaindertermsoftheexpansion. Theaveragepopulationactivitydynamicsof
a finite size network can be described approximately in terms of the following
Ornstein-Uhlenbeck process
(cid:0) (cid:1) σ(t)
dn¯(t)≈ −m(t)+F(m(t)) dt+ √ dB (13)
t
N
whereσ2(t):=m(t)(1−2g(µ (t))+g(µ (t))andB isBrownianmotion. Inthe
1 1 ·
approximation of eqn. 13, we ignore the contribution of remainder terms (e.g.,
R) to σ(t). It now becomes clear, that the fluctuating term of eqn. 13 vanishes
in the thermodynamic limit.
3 Convergence to MFT in networks of binary
neurons
In order to illustrate our complete MFT, we consider two scenarios that are
relevant to the theoretical analysis of neural systems. The first scenario is that
of a population network with a constant external input µ > 0. When J¯< 0
0
this system is the classical balanced network for which the external input is
canceled by internal fluctuations. We choose a widely-used gain function in
neural networks theory which it is given by
1+Erf(αx)
f(x):= . (14)
2
The parameter α describes the intrinsic noise intensity of the individual units
and therefore must be positive. When α → ∞, this transfer function approxi-
mates to the well-studied Heaviside step function [18]. Using the transfer func-
tion given by eqn. 15 (with α = 5) and a directed fixed-in-degree Erd˝os-R´enyi
6
y 0.40
t
vi
ti 0.004
ac 0.35 S
M

p. R
po 0.30 0.000
. 1.0 0.3
e
v
a 0.25

.
uil
q 0.20
e
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3
coupling strength J¯
Figure 1: Convergence of the average population activity to the steady-state
MFT predictions. The red line indicates the predictions of complete MFT
in eqn. 3 and is computed as described in the Materials and Methods section.
The blue line is the predictions of mean-field eqn. 12, assuming only Gaussian
fluctuations. TheblackdotsrepresentsimulationsofanetworkwithN =1000
units averaged over 20 independent trails (error bars are smaller than symbol
size). The inset is the root mean square of error between simulations and
the complete theory (upward red triangles) and the Gaussian approximate
theory (downward blue triangles). The simulations were performed using a
Doob-Gillespie algorithm for T =5×105 steps with the gain function given by
eqn. 15 . The averaged activity was estimated in the last 5×103 steps across
all trials. Parameters: γ =0.5, α=5, K =10 and µ =0.1.
0
7
network (with K =10), we compute the complete steady-state mean-field limit
of eqn. 3, using the approximation described in the Materials and Methods
section (Fig.1, red line). We compare the complete mean-field theory (Fig.1,
red line) with the mean-field prediction that assumes only Gaussian statistics
(i.e., K →∞) in eqn. 12 (Fig.1, blue line). The difference between the predic-
tions becomes apparent as |J¯| increases. Numerical simulations of a finite-size
network (N = 1000) are used to estimate the steady-state population activity
byaveraging20independenttrials(Fig.1,blackdots). Theequilibriumpopula-
tionaverageactivityofsimulatednetworks(Fig.1,blackdots)exhibitsexcellent
agreement with both the complete (Fig.1, red line) and the Gaussian approx-
imation (Fig.1, blue line) of the mean-field limit in the case of weak coupling.
However, in cases where the coupling is strong, the average population equi-
librium activity deviates from the Gaussian approximation (Fig.1, blue line)
and,instead,followsthepredictionsofthecompletemean-fieldlimit(Fig.1,red
line). Therefore, the Gaussian approximation that is given in eqn. 12 is only
reasonableforweakcouplingandrelativelylargevalueofK. Theerrorbetween
steady-statepopulationactivityfromthesimulations(Fig.1,blackdots)andthe
Gaussianapproximation(Fig.1,blueline)increasesas|J¯|becomeslarger(Fig.1,
downward blue triangles in the inset). However, the error between averaged ac-
tivity of the simulated networks and the complete MFT stays almost constant
(Fig.1, upward red triangles in the inset); this error is due to the truncation of
the series expansion and the numerical stability of power-series representation
underlies the function F in eqn. 7 (refer to Materials and Methods section for
more details).
In order to illustrate our complete MFT, we consider two scenarios that are
relevant to the theoretical analysis of neural systems. The first scenario is that
of a population network with a constant external input µ > 0. When J¯< 0
0
this system is the classical balanced network for which the external input is
canceled by internal fluctuations. We choose a widely-used gain function in
neural networks theory which it is given by
1+Erf(αx)
f(x):= . (15)
2
The parameter α describes the intrinsic noise intensity of the individual units
and therefore must be positive. When α → ∞, this transfer function approxi-
mates to the well-studied Heaviside step function [18]. Using the transfer func-
tion given by eqn. 15 (with α = 5) and a directed fixed-in-degree Erd˝os-R´enyi
network (with K =10), we compute the complete steady-state mean-field limit
of eqn. 3, using the approximation described in the Materials and Methods
section (Fig.1, red line). We compare the complete mean-field theory (Fig.1,
red line) with the mean-field prediction that assumes only Gaussian statistics
(i.e., K →∞) in eqn. 12 (Fig.1, blue line). The difference between the predic-
tions becomes apparent as |J¯| increases. Numerical simulations of a finite-size
network (N = 1000) are used to estimate the steady-state population activity
byaveraging20independenttrials(Fig.1,blackdots). Theequilibriumpopula-
tionaverageactivityofsimulatednetworks(Fig.1,blackdots)exhibitsexcellent
8
0.28
)
t
(
¯n 0.26

y
t
vi 0.24
ti
c
a 0.22

.
p
o 0.20
p
. 0.18
e
v
a 0.16
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8
time 1e8
Figure 2: Emergence of stochastic MFT. The green line shows the temporal
evolutionofasimulatednetworkwithN =5000andK =10; thisnetworkdoes
nothaveadeterministicMFTsincetheconditionineqn.10isnotsatisfied. The
out-degreeofasingleunitinthenetworkwassettobeN (i.e.,ρ=1ineqn.16).
Forcomparison,thered lineshowstheaveragepopulationactivityofasimilar
network for which ρ=K/N. The system time in x-axis is calculated according
to a Doob-Gillespie algorithm (the average update interval is dt=1368.6, refer
to the Materials and Methods section for details). Parameters: J¯= −0.7 and
all other parameters are as in Fig.1.
agreement with both the complete (Fig.1, red line) and the Gaussian approx-
imation (Fig.1, blue line) of the mean-field limit in the case of weak coupling.
However, in cases where the coupling is strong, the average population equi-
librium activity deviates from the Gaussian approximation (Fig.1, blue line)
and,instead,followsthepredictionsofthecompletemean-fieldlimit(Fig.1,red
line). Therefore, the Gaussian approximation that is given in eqn. 12 is only
reasonableforweakcouplingandrelativelylargevalueofK. Theerrorbetween
steady-statepopulationactivityfromthesimulations(Fig.1,blackdots)andthe
Gaussianapproximation(Fig.1,blueline)increasesas|J¯|becomeslarger(Fig.1,
downward blue triangles in the inset). However, the error between averaged ac-
tivity of the simulated networks and the complete MFT stays almost constant
(Fig.1, upward red triangles in the inset); this error is due to the truncation of
the series expansion and the numerical stability of power-series representation
underlies the function F in eqn. 7 (refer to Materials and Methods section for
more details).
The hierarchy of conditions on the connectivity matrix J (eqns. 10 and 11)
guarantees the convergence of the average population activity to the prediction
of MFT. In particular, eqn. 10 indicates that, as N → ∞, the average column
sum of the connectivity matrix J should be K. It is straightforward to con-
struct networks that do not obey this rule; such networks lose their pointwise
convergence to a deterministic mean-field limit in eqn. 3. An extreme example
9
of a network of this kind is a network that that has a single unit, n , that
j∗
connects into ρN neurons in the circuits, where 0 < ρ ≤ 1 is the fraction of
units in the network that are post-synaptic for n . Numerical simulations of
j∗
such a network (N = 5000 and ρ = 1) show large-amplitude population activ-
ity fluctuations (Fig.2, green line), in contrast to the smaller fluctuations of a
homogeneous network (Fig.2, red line). Our approach allows the construction
of stochastic correction terms to the mean-field limit by isolating the unit n
j∗
from the network and then taking the limit N → ∞. Therefore, a first-order
correction to the function F of eqn. 7 can be derived as
F (m(t))≈F(m(t))+ρJ¯K−γf(cid:48)(µ (t))n (t). (16)
s 1 j∗
F isastochasticfunctionsincen isabinomialrandomvariableforwhichthe
s j∗
probability of being at state one is m(t); the mean-field equation is thus trans-
formed into an ordinary stochastic differential equation. The correction term
in eqn. 16 indicates that the observed large fluctuations (Fig.2, green line) are
indeed a finite K phenomenon. Therefore, in large networks that have a finite
number of connections between units (e.g. finite K networks), it suffices that
onlyoneunitbreaksthecondition(i.e.,ρ>0)and,asaresult,thedeterministic
MFT collapses. The emergence of large-amplitude population events in Fig.2
has been observed previously as the indication of the “synfire chain” in cortical
networks[2]. Itisnoteworthythatthereiscompellingevidencethatsuchstruc-
tures exist in cortical microcircuits [9] and recent experimental results indicate
that the diverse couplings between single-cell activity and population averages
[11] suggest the existence of this polychronous state [8] in cortical networks.
4 Discussion
In this paper, we considered a simple network of N binary units of which the
recurrent interactions can be in balance with the excitatory feedforward input.
This simplified model captures the essential aspects of cortical asynchronous
state [18], and it allowed us to demonstrate, with extra clarity, the calculation
of a complete MFT for finite-size systems dynamics with arbitrary connectiv-
ity structure. The model neurons are two-state units that update their state
stochastically. Theupdatedstateofasingleunitiseitheractivewithaprobabil-
ity given by f(n(t)) orin a state ofquiescence. Our results show that the MFT
for binary units is fundamentally constructed from the Law of Large Numbers
(LLNs) and the emergence of intrinsic fluctuations does not require the appli-
cation of the central limit theorem. We further demonstrate the dependence
of averaged population activity fluctuations on both the system size, N, and
the average number of connections per unit, K. Furthermore, our formalism
allows the study of a finite-size network in the form of an ordinary stochastic
differential equation (eqns. 13 and 16). Noteworthy, the method presented here
reliesonourabilitytocalculatethefunctionF ofeqn.7. Amoreefficientrepre-
sentationforfiniteK networks, similartothelargeK limitineqn.12, provides
a further means to investigate systems with very strong coupling schemes.
10