Mean Field Theory For Non-Equilibrium Network Reconstruction

0.27 MB·English

by  Yasser Roudi

Mean Field Theory For Non-Equilibrium Network Reconstruction
Yasser Roudi1 and John Hertz1,2
1NORDITA, Stockholm, Sweden
2The Niels Bohr Institute, Copenhagen, Denmark
There has been recent progress on inferring the structure of interactions in complex networks
when they are in stationary states satisfying detailed balance, but little has been done for non-
equilibrium systems. Here we introduce an approach to this problem, considering, as an example,
the question of recovering the interactions in an asymmetrically-coupled, synchronously-updated
Sherrington-Kirkpatrick model. We derive an exact iterative inversion algorithm and develop effi-
1
cientapproximationsbasedondynamicalmean-fieldandThouless-Anderson-Palmerequationsthat
1
expressthe interactions in terms of equal-time and one time step-delayed correlation functions.
0
2
PACSnumbers: 05.10.-a,02.50.Tt,75.10.Nr

n
a
J Introduction.—Findingthe connectivityincomplex changingthetemperaturecanberealizedbychangingthe

5 networks is crucial for understanding how they oper- coupling parameter g and the field strengths. Even for

ate. Gene and multi-electrode microarrayshave recently time-independent fieldandina steadystate,this system
]
n madethetypeofdatarequiredforthispurposeavailable. is not in a Gibbs equilibrium [10]. However, we show
n What is needed now is appropriate theoretical tools for that, like its equilibrium counterpart, the nonequilib-
- analyzing these data and extracting the connectivity. rium inverse problem for this model can be solved using
s
i In much recent work on this subject [1–3], the prob- a gradient descent method and also via systematic ap-
d
lem has been posed as that of inferring the parameters proximateinferencesderivedusingdynamicalversionsof
.
t
a of a stationary Gibbs distribution modeling the system. naive mean-field (nMF) and Thouless-Anderson-Palmer
m While satisfied in many applications, the assumption of (TAP) equations. We show that for both the stationary
- Gibbs equilibrium is unlikely to hold for many biological and nonstationary systems these methods provide effi-
d systems since they are usually driven by time-dependent cient reconstructionof interactions. We alsoanalytically
n
externalfields,theirinteractionsmaynotsatisfydetailed quantify their errors.
o
balance, or they may only be observed while the tran- Exact, nMF and TAP learning rules.— Sup-
c
[ sients dominate the dynamics. Applying the equilibrium pose that we have observed R realizations of duration

approach to such cases usually yields effective interac- L time steps of the process in (1). We denote the ob-
2
v tionsthatdonotbearanobviousrelationshiptothereal served state of the system at time t of realization r by
6 ones [3]. Kinetic and nonequilibrium models provide a sr(t)={sr(t),··· ,sr (t)}. Tofindthecouplingsandex-
1 N
4 much richer platform for studying such systems [4–6]. ternal fields, we maximize the likelihood of the observed
9 Whereas for equilibrium models the development of states under the model (1). This maximization can be
5
systematic mean field inference methods [7] has led to done using an iterative algorithm, analogous to Boltz-
.
9 great practical and conceptual advancements, a mean mann learning for the equilibrium model: starting from
0
field theory for nonequilibrium network reconstruction aninitialsetofcouplingsandfields, oneadjuststhemit-
0 is still lacking. In this paper, we show how a mean field erativelybystepsofsizesδh =η ∂L andδJ =η ∂L ,
1 i h∂hi ij J∂Jij
: theory for inference canalso be developedfor a nonequi- L being the log-likelihood. The learning steps thus are
v
libriumsystem. Weconsiderthisproblemforaparticular
i
X simple nonequilibriummodel: a kinetic Isingmodelwith δhi(t)=ηh(cid:8)hsi(t+1)ir−htanh[θi(t))]ir](cid:9) (2a)
ar randomasymmetricinteractions(JjiindependentofJij), δJij =ηJ(cid:8)hsi(t+1)sj(t)i−htanh[θi(t)]sj(t)i(cid:9) (2b)
inanexternalfieldwhichmaybetime-dependent. Thisis
a discrete-time, synchronouslyupdated model composed where η and η are learning rates. Here and in what
h J
of N spins si =±1 with transition probability followsh···ir,h···irepresentaveragingoverrepeats,and
bothrepeatsandtime,respectively. Anoverline,instead,
exp[s (t+1)θ (t)]
Pr(s(t+1)|s(t))=Y i i (1) will indicate averaging over the spins. One can think of
2cosh(θi(t)) Eq. (2b) as performing a logistic regression to explain
i
one-step separated correlations. This is similar to what
where θi(t) = hi(t) + PjJijsj(t). The couplings Jij is proposed in [13] as an approximation for inferring the
are independent Gaussian variables with variance g2/N. connectivity in an equilibrium Ising model.
This model can be readily applied to time-binned neu- Since performing the steps in this algorithm does not
ral data, where t labels the bins, and s (t) = ±1 repre- requireMonteCarloruns,itisfasterthantheequilibrium
i
sents a spike or no spike by neuron i in bin t [1]. The Boltzmannlearning. However,twofactorsstillmakethis
temperature has been set equal to 1, since any effects of algorithm slow for large systems and/or data sets, war-
2
ranting the development of fast approximations. First, (the one approaching zero as g → 0). This root cannot
(2) is still aniterative algorithmwhichcould take a long exceed 1/3, restricting this technique to weak couplings.
timetoconvergeifnotprovidedwithagoodinitialcondi- For both nMF and TAP reconstruction, the external
tionandlearningrates. Second,ateachsteptheaverages fieldsh canalsobefoundbysolvingtherespectivemag-
i
on the right hand side of (2) should be calculated from netization equations after the J have been obtained,
ij
the data de novo, given the adjusted parameters. just as in the equilibrium problem [7].
Two fast approximations, nMF and TAP learning Performance of the algorithms.— We have veri-
rules,arederivedandstudiedbelow. Toimplementthem fied that the algorithm (2) recovers the couplings of an
in the stationary case, one first uses the data to calcu- asymmetricSKmodelexactlyinthe limitofL→∞, for
late the one-step delayed and equal time correlations, a wide range of coupling strengths g, external fields and
D = hδs (t+1)δs (t)i and C = hδs (t)δs (t)i, where system sizes. The mean square error, ǫ , is in gen-
ij i j ij i j exact
m =hs i and δs =s −m . The approximations are eral proportional to 1/L, and in the weak-coupling limit
i i i i i
a quadratic expansion of log-likelihood yields
JnMF/TAP =AnMF/TAP−1DC−1 (3)
1
ǫ =δJ2 ≡(J −Jtrue)2 = , (7)
whereAnMF =(1−m2)δ ,ATAP =AnMF(1−F )andF exact ij ij ij (1−m2)L
ij i ij ij ij i i i
is the root of the cubic equation (6) below. In the non-
where J (Jtrue) are the inferred (true) couplings.
stationarycase too, similarlearning rules canbe derived ij ij
WefindthatthenMFalgorithmleadstoanerror,ǫ ,
as shown later in the paper. MF
of the form ǫ +ǫ∞ , where ǫ∞ is independent of
Derivation of nMF and TAP inversion.— For exact nMF nMF
Landproportionalto1/N. Thus,fordatasetsoflength
simplicity,weconsiderfirstthestationarycase,forwhich
L≪L∗ =1/ǫ∞ ∝N, nMF does almost as well as the
the sequence index r is superfluous, as averaging over nMF
exact algorithm. Furthermore, the larger the network,
time and repeats would be equivalent. We start with
the better nMF does. The errorsfor the exact and nMF
the maximum likelihood conditions, i.e. δh = δJ = 0
i ij
algorithms vs L are shown in Fig. 1a.
in (2). Using the nMF equations m = tanh(h +
i i
JnMFm ), and writing the s in (2) as m +δs , we
Pj ik k i i i a) b)
expand the tanh in the δs . The first nonzero term gives
i
1E-5 1E-5
hδsi(t+1)δsj(t)i=(1−m2i)Xk JinkMFhδsk(t)δsj(t)i. (4) εnMF gg==00..1146εTAP
1E-7 g=0.12 1E-7
g=0.1
which can be written as (3) for the nMF case. g=0.16
g=0.14
To get the TAP inversion formula, we start instead g=0.12
1E-9 1E-9 g=0.1
by assuming that the mi satisfy the TAP equations 1E4 1E6 L 1E8 1E10 1E4 1E6 L 1E8 1E10
m =tanh[h + JTAPm −m (JTAP)2 (1−m2)],
i i Pk ik k iPk ik k
whichtakeintoaccounttheOnsagerreactionterm. Kap-
FIG. 1. (Color online) Performance of thealgorithms. Exact
pen and Spanjers [8] provedthat the TAP equations, al- andnMF(a). andtheTAP(b)errosareshownvsdatalength
though usually derived for the equilibrium (symmetric- L for g = 0.1 (blue stars), 0.12 (magenta crosses), 0.14 (red
J) SK model, also hold for the asynchronously updated, circles) and 0.16 (black x), all for N = 20 and zero external
asymmetric-Jmodel in a stationary state. We have veri- field. Theoretical predictions are thesolid lines.
fiedthattheyarealsovalidinoursynchronously-updated
model [9]. We again write s = m + δs , expand- For weak coupling, we can calculate the asymptotic
i i i
ing the tanh to third order in powers of PkJiTkAPδsk + nMFerror,ǫ∞nMF, analyticallyasfollows. We presentthe
m (JTAP)2 (1−m2). Keeping terms up to order g3 zero-field case here for simplicity. We expand the tanh
iPk ik k
leads to D=ATAPJTAPC , where in the max-likelihood equation to third order, giving
ATijAP =AnijMF[1−(1−m2i)X(JTAP)2il(1−m2l)]. (5) Din =XJikhsksni− 31XJikJilJimhskslsmsni+··· .
l k klm
(8)
These equationscannotbe solveddirectlyas inthe nMF Correlationshereareatequaltimes,exceptforD . The
in
case because ATAP depends on JTAP. However, one dominant contributions in the sum over k,l,m are those
can derive a cubic equation for the quantities Fi = withk=l,l=mandm=k. Multiplyingontherightby
(1−m2i)Pl(JTAP)2il(1−m2l): (C−1)nj, summing over n and using (3) for nMF, yields
Fi(1−Fi)2 =(1−m2i)X(JnMF)2ij(1−m2j). (6) JinjMF =Jij −XJi2kJij, (9)
j k
This determines ATAP = AnMF(1−F ), yielding (3) for withcorrectionsofrelativeorder1/N. Eq. (9)alsoyields
ij i
theTAPcase. Therelevantrootof (6)isthesmallestone the TAP-approximationcouplings found above, showing
3
that the TAP reconstruction indeed corrects the leading where B(i) = h(1−m2(t+1))C (t)i . The problem is
kj i kj t
MFerrors. To leadingorderthe sumonk is justg2, and more complex than the stationary one only because one
the asymptotic mean square MF error is has to invert a different matrix B(i) for each i.
For TAP, analogously to the stationary case, the B(i)
g6
ǫ∞ =(J −JnMF)2 = . (10) acquire an extra factor inside the time average:
nMF ij ij N
The solid curves in Fig. 1a are 1/L+g6/N; the fit is Bk(ij) =h(1−m2i(t+1))(1−Fi(t))Ckj(t)it, (15a)
evidently good. As shown in Fig. S1 [11], nMF exhibits Fi(t)=(1−m2i(t+1))X(JTAP)2il(1−m2l(t)). (15b)
a systematic error by underestimating the magnitude of
l
the couplings. The factor 1−F in TAP formula cor-
i
rects for this to relative order g2. Thus, when one is in- Exact TAP inversion requires iterative solution of (14),
withJTAP insteadofJnMF,togetherwith(15). Wehave
terested only in the presence or absence of connections, ij ij
found, however, that effective reconstruction is still pos-
there would be little difference between nMF and TAP.
sible under the simplifying approximation that F (t) in
The error for the TAP reconstruction is much lower i
Eq. (15a) can be represented by its temporal mean. In
than that of the nMF one and reaches its minimum at
this case, F ≡hF (t)i solves the cubic equation
much larger L: for N = 20 and the coupling strengths i i t
westudied, wehadto goto L∼109 toseethe errorflat-
ten(Fig.1b). Tocalculatetheasymptoticreconstruction Fi(1−Fi)2 =X(JnMF)2ijh(1−m2i(t+1))(1−m2j(t))it.
j
errorfor TAP,we expandthe tanh to 5thorderandpro-
ceed to evaluate the averages as we did for nMF. The Solving it and using it in Eq. (15a), one can calculate
nMFerrorterms analyzedabovearecompensatedforby JiTjAP = JinjMF/(1−Fi). Similar to the stationary case,
the TAP equations, as N → ∞, leading to an asymp- afterinferringthe couplings,onecanusetheforwarddy-
totic ǫ∞ = 4g10/N. For N ≫ 1/g2 this is the leading namical nMF and TAP equations Eqns. (11) and (12)
TAP
term in the asymptotic TAP error. Outside this regime, to infer the time-varying external field. The result of
afinite-sizeeffectshouldalsobetakenintoaccount. This
isbecauseinmakingthatTAPcorrection,thetermin(8) a) b)
0.1 0.1
(fwr9oi)tmhinktsht=eisaldo=voefmrcoohnuacnest.binTegehniescmo(u2e/na3tne)2dsJqti6hujra=eree(te2irm0rgoe6rs)i/tnh(3aoNtbtr3ae)isnuainnltdgs nMFJ(nonstat.)ij−00..00055 nMFJ(stat.)ij−00..00055
should be added to the 4g10/N term. −0.1 −0.1
Non-stationary case.— The magnetizations, −0.1−0.05 0 0.05 0.1 Jtrue −0.1−0.05 0 0.05 0.1 Jtrue
ij ij
m (t)=hsr(t)i , are now time-dependent and, for nMF,
i i r
solve
C)0.5
mi(t+1)=tanh[hi(t)+XJinjMFmj(t)]. (11)
j h 0
We have also proved [9] that the TAP equations hold −0.5 time
100 150 200 250 300
even in a nonstationary state, in the form
mi(t+1)=tanh[hi(t)+XJiTjAPmj(t) FIG. 2. (Color online) Inference in the nonstationary case.
(a) Couplings of a network of N =20 driven by a sinusoidal
j
external field inferred using the nonstationary nMF, and (b)
−mi(t+1)X(JTAP)2ij(1−m2j(t))]. (12) thestationarynMF.(c)Twoperiodsoftheexternalfield(thin
j bluefullcurve)anditsreconstructionusingthenonstationary
nMFcouplings(reddashedcurve)andstationarynMF(thick
Thus, we can extend both our inversion algorithms to
black full curve).
nonstationary systems, as we show in the following.
We start by defining time-dependent correlation ma-
reconstructing the couplings of a network driven by a
trices D (t) ≡ hδsr(t + 1)δsr(t)i and C (t) ≡
ij i j r ij common sinusoidal external field to all spins is shown in
hδsr(t)δsr(t)i . FornMF,usingthesameprocedurethat
i j r Fig. 2. Fig. 2a shows how well the couplings are inferred
lead to (4), we find
by nonstationary MF using L=105 and R=100. Non-
hDij(t)it =XJinkMFh(1−m2i(t+1))Ckj(t)it. (13) sstqautaiorendareyrrToAr:P6c.o7u×pl1in0g−s7(nveortssuhso1w0n−)6hfaovrenaMloFw.eIrnmFeiagn.
k
2b, we also plot the couplings inferred using stationary
One can still solve for J by simple matrix algebra:
nMF inversionfor each of the 100 repeats and averaging
JinjMF =XhDik(t)it[(B(i))−1]kj, (14) over them. Not surprisingly, the stationary nMF per-
forms poorly on this nonstationary data. Importantly,
k
4
there is a systematic overestimation of the couplings in resulting distributions are completely overlapping when
thiscase,becausethestationarymethodaccountsforcor- and equilibrium TAP learning is used. When using a
relations induced by the common, time-varying external model like (1) to infer connectivity in a system with a
field through adjusting the couplings. Correspondingly, different dynamics, or when faced with data limitation,
if one uses the couplings inferred by stationary nMF in including prior knowledge about the network could be
(11) to infer h (t), the amplitude of this field is underes- very beneficial. In particular, taking into account spar-
i
timated, while the use of nonstationary nMF couplings sity of the connections via a l-1 regularizer added to the
yields a very good reconstruction of h (t); see Fig. 2c. likelihood has been shown to be very useful [13]. It is
i
Discussion.— We have shown how to infer interac- easy to show that adding an l-1 regularizer to the likeli-
tions in a simple but nontrivial nonequilibrium system: hoodofthedataunder(1)wouldmodify (3)byaddinga
akineticIsingmodelwithrandomandpotentiallyasym- term proportionalto AnMF/TAP−1sgn(J)C−1 to the right
metric interactions. The model is the maximum entropy hand side. How this improves inferring connections in
modelfor eachtimestep,givenmeanmagnetizationsand biological networks will be discussed elsewhere.
one step separatedcorrelations. We havedescribed both
Asimpleextensionof (1)isitscontinuoustimeversion.
an exact iterative algorithm and two approximate ones,
Asshownin[15],forthismodel,too,ameanfieldtheory
based ondynamical nMF and TAP equations,which are
can be developed using the approach presented here. In
correct up to corrections of order 1/N. We calculated
otherrecentkinetic approachestoproblemslikethis,the
analytically the errors of these approximations for weak
equilibrium maximum-entropy approach [1] is extended
coupling. The method shows particular promise when
toincludenon-equal-timecorrelations[5]andanapprox-
applied to nonstationary states, where it separates true
imate scheme for fitting an integrate-and-firenetwork to
interactions from the apparentones found by applying a
data was developed in [4]. There has also been work [6],
stationary theory to a nonstationary state.
closelyconnectedto(1),inwhichs (t+1)dependsonlin-
i
ear combinations of h(t′) and s(t′), for t′ ≤t. Given the
a)x 104 b)x 104
advantage of these nonequilibrium models over the equi-
6
8
libriumonesfordescribingspiketrainstatistics,amean-
6 4 field theory for inferring their parameters would be of
4
great theoretical and practical benefit. For such models,
2
2 weexpectthatitwillbepossibletousethetechqniuesin
0 −0.2 −0.1 0 0.1 0.2 0−0.2 −0.1 0 0.1 0.2 [8] or [10, 16] to derive dynamical nMF and TAP equa-
JTijA(Pnonequib.) JTijA(Pequib.) tions. Employing the approach developed here one can
then build approximate mean field inversion techniques
FIG. 3. (Color online) Finding connections in a cortical net- based on these equations.
work model. (a) The histogram of the couplings inferred us-
ing the stationary nonequilibrium TAP for pairs of neurons
thatwereconnected(bluefullbars),andthosethatwerenot
(red empty bars). The separation between the histograms
shows that one can use the TAP approximation to separate
[1] E. Schneidman, et al, Nature, 440, 1007 (2005);
connected and disconnected pairs. (b) same as (a) for equi-
J.Shlens,etal,J.Neurosci.,26,8254(2006); G.Tkacik,
librium TAP.
et al,arXiv:0912.5409v1 [q-bio.NC], (2009).
[2] M. Weigt, et al, PNAS, 106, 67 (2009); T. Lezon, et al,
AkineticIsingmodelwillshowanintrinsicerrorwhen
ibid., 103, 19033 (2006).
appliedtodatafromadifferentkindofsystem. However, [3] Y. Roudi, J. Tyrcha,J. A. Hertz, Phys. Rev. E, 79,
even when applied to data from a realistic network, the 051915 (2009).
simple approximate learning rules developed here iden- [4] S. Cocco, et al, PNAS,106, 14058 (2009).
tify the connections much better than their equilibrium [5] O. Marre, et al, Phys. Rev.Lett., 102, 138101 (2009).
[6] W. Truccolo, et al,J. Neurophys.,93, 1074 (2005).
counterparts. Fig. 3 showsthe distribution ofcouplings
[7] H. J. Kappen and F. B. Rodriguez, Neur. Comp., 10,
found by applying the nonequilibrium TAP learning to
1137 (1998); T. Tanaka, Phys.Rev. E, 58, 2302 (1998).
data from a simulated model cortical column with in-
[8] H.J.KappenandJ.J.Spanjers,Phys.Rev.E,61,5658
hibitory and excitatory neurons [12]. The connections (2000).
in the model were dilute with 10% probability of con- [9] Y. Roudiand J. Hertz, In preparation (2010).
nection. When there is no synapse from neuron j to i, [10] A. C. C. Coolen, arXiv:cond-mat, 0006011 (2000).
the inferred J follows a zero mean distribution, while [11] See EPAPS Document No. [] for details about the com-
ij
putersimulations and Fig. S1.
if there is an excitatory/inhibitory synapse, it follows a
[12] J. A.Hertz, Neural Comp., 22 (2010).
positive/negativemeandistribution,wellseparatedfrom
[13] P. Ravikumar,et al,Annals. Statistics, 38 (2010).
the first one. One can thus easily use the distribution of
[14] J. A.Hertz, et al, BMC Neuroscience, 10 (2010).
inferredcouplingstoinferthepresence,absenceandsign [15] H.-L Zenget al.,arXiv:1011.6216v1 [stat.CO], (2010) .
oftheconnections;see[11]and[14]. Onthecontrary,the [16] G. Biroli, J. Phys.A: Math. and Gen., 32 (1999).