Inferring a network from dynamical signals at its nodes

We give an approximate solution to the difficult inverse problem of inferring the topology of an unknown network from given time-dependent signals at the nodes. For example, we measure signals from individual neurons in the brain, and infer how they are inter-connected. We use Maximum Caliber as an inference principle. The combinatorial challenge of high-dimensional data is handled using two different approximations to the pairwise couplings. We show two proofs of principle: in a nonlinear genetic toggle switch circuit, and in a toy neural network.


I. LEARNING THE PROPERTIES OF A NETWORK FROM MEASUREMENTS AT ITS NODES
We are interested in solving the following 'inverse problem': you measure time-dependent signals from individual agents.Those agents behave in a correlated way.That is, they are connected in some network that is unknown to you.The goal is to infer the interactions between these agents, from their correlations.For example, measure the protein concentrations that are produced from an unknown gene network, and infer the degree to which the proteins activate or inhibit each other.Or measure the firings of individual neurons and infer the neuron-neuron connection strengths in the brain.These problems are the 'inverse' compared to the common situation of knowing a network and computing the flows through it.
The method of choice for inferring dynamical processes from limited information is the Principle of Maximum Caliber (Max Cal) [1][2][3][4][5][6][7].Maximum Caliber (Max Cal) is to dynamics what Maximum Entropy inference (Max Ent) is to equilibrium.Max Cal is a procedure that predicts full stochastic distributions by maximizing the route entropy subject to the constraint of a few known average rates.
The challenge here is that the number of possible interactions (here the node-node couplings) grows rapidly with system size (the number of nodes in the network and length of time of observing the signals).So, direct implementation of Max Cal is limited to small or simplified systems [8][9][10][11].For larger and more realistic situations, this Max Cal inference procedure becomes computationally intractable.In other matters of physics, the dimensionality of the problem is reduced by various approximations, including variational approximation and perturbation theory [12][13][14][15].These techniques have been used to reduce the dimensionality in other high-dimensional inference problems [16][17][18][19][20][21][22][23][24][25][26][27].
Here, we adapt such methods for inferring high-  dimensional, heterogeneous dynamical interrelationships from limited information.Related generalizations have been previously used to infer dynamical interactions in continuous-time Markovian networks [28,29].However, these approaches make strong assumptions either about the dynamics or about unknown transition rates.Here instead, with Max Cal, we can infer both the dynamics and interactions within arbitrarily complex, nonequilibrium systems, albeit in an approximate way.We describe two different levels of approximation: Uncoupled and Linear Coupling.

II. THE PROBLEM
Suppose you run an experiment and record the activity of N arbitrarily interacting agents (the nodes, i = 1, 2, . . ., N of a network) over some time T .
The data arrives as a time series: v(t) = {v 1 (t), v 2 (t), . . ., v i (t), . . ., v N (t)}, also called a trajectory, Γ (from t = 0 to T ).From the signals, we aim to predict the coupling strengths between the nodes.Our model should reliably predict certain averages over the data, with otherwise the least possible bias.Such problems are the purview of the principle of Maximum Entropy or its dynamical extension, Maximum Caliber (Max Cal) [1][2][3][4][5][6].The principle of Max Cal chooses the unique distribution, P Γ , that maximizes the path entropy, or Caliber, over all acceptable distributions {P Γ }, while respecting the observed constraints.The constraints here are the mean value, M i (t) over all possible paths, and the correlations, χ ij (t, s): for all times t and s over all agents i and j.The Caliber is expressed as where the sum over Γ is a sum over all the possible realizations of the full time series.Here h i (t) and K i,j (t, s) are the Lagrange multipliers that each enforce the constraints in Eq.1.The other Lagrange multiplier, µ, ensures the distribution is normalized (the probabilities sum to one).
Maximizing the Caliber with respect to all possible distributions {P Γ } gives (3) Z is the dynamical partition function, the quantity that normalizes P Γ .By analogy with the Ising model for equilibrium systems, h i (t) represents the strength of the external fields to which the system is coupled, whereas K ij (t, s) are the couplings between the components of the system.

III. THE UNCOUPLED MAX CAL APPROXIMATION
We aim to compute h i and K ij for every time point.This presents a combinatorial challenge for large networks or long trajectories.We describe two levels of approximation to overcome this challenge.In the present section, we describe our simplest approximation, representing a mean-field or Uncoupled approach, which allows us to solve even large systems [23,29].This method works by breaking the full inference problem into simpler, independent subproblems.For our application, this suggests that we try uncoupling the trajectories of each object (i) which we denote Γ i .The approximate trajectory distribution Q Γ then factorizes into the product: Eq. 3 shows that this approximation is exact when all of the coupling coefficients K ij (t, s), i = j, are 0. We can force this condition by temporarily ignoring all pairwise constraints corresponding to K ij (t, s) and satisfying the remaining, i = j, Max Cal constraints (from Eq. 1).The now uncoupled distributions are given by: (5) This then gives a new set of effective Lagrange multipliers, hi (t) and Kii (t, s), which absorb the average effects of the neglected pairwise interactions.
In summary, this Uncoupling approximation reduces the inference problem to solving independent single-node problems for each i.These single-node inference problems are readily solved [5,9].Clearly, however, this naive approximation fails to capture any pairwise correlations between agents (i = j).The following section describes a next better approximation, based on Linear Response Theory [19].

IV. THE LINEAR COUPLING MAX CAL APPROXIMATION
Here, we go beyond the uncoupling assumption and take the first-order perturbation term around our Uncoupled approximation.We call this the Linear Coupling Max Cal approximation.The first-order approximation for the Lagrange multipliers for each agent i are given by (see Appendix B): Eq 6 -analogous to familiar mean-field models in physics -attempts to recover the true Lagrange multipliers (with ' denoting the Linear Coupling approximation) from the effective, uncoupled Lagrange multipliers (denoted by ∼) [23,29].Thus our only remaining unknowns are the values of the pairwise couplings K ij (t, s).For our first order approximation, the Linear Coupling estimates for these Lagrange multipliers have a closed-form solution (see Methods, eq. 9) given by: where is the matrix of covariances.Once these estimates are known, the remaining Lagrange multipliers are easily found from Eq. 6.
Below, we give two examples to illustrate two points.First, we show that even the Uncoupled approximation can give a fairly accurate closed-form solution for a network with nonlinear interactions.We show this for a genetic toggle switch [8,30].Second, we show how the  Linear Coupling approximation readily handles a highdimensional heterogeneous system, namely a toy network of neurons, which is otherwise computationally intractable.

V. FINDING STABLE STATES IN THE GENETIC TOGGLE SWITCH
Collins et al. engineered a synthetic circuit into a single-celled organism called a genetic toggle switch [30].It consists of two genes (A and B, blue and yellow in Figure 2a).Each gene produces a protein that inhibits the other.So, in analogy with an electrical toggle switch, either A is being produced while B is turned off, or vice versa.This network has previously been computationally simulated using Max Cal [8], so it provides a 'Ground Truth' for comparison with our approximation here.The present exercise shows that Uncoupled Max Cal, which can be solved analytically, gives an excellent approximation to the nonlinearity and the phase diagram in this known system.Importantly, beyond this proof of principle, Uncoupled Max Cal is readily applicable to bigger more complex systems.
Here, our input data takes the form of the stochastic numbers of protein molecules, totaling N A (t) and N B (t) on nodes A and B at time t.Our trajectories are the counts of the numbers (l α and l β ) of newly produced proteins of types A and B respectively over each new short time interval δt.From these trajectories, we use Max Cal to compute three quantities: the production rate of each protein, the survival rates (the count of proteins that are not degraded), and the strength of the negative feedback.To keep the model simple, we suppose that both proteins have the same production rate, and both have the same survival rate.From the data, we obtain the average production and survival rates, h P and h S , which are enforced in the Max Cal modeling as Lagrange multipliers.And, from the data, we obtain the correlation between production of A and survival of B, l α l B (and vice-versa); these are enforced by a third Lagrange multiplier, K, the coupling coefficient [8] (see Appendix C for details).
The behavior of this network is known from the Ground Truth simulations; see Fig 2b .There is a critical value, K c < 0 of the coupling parameter (or repression strength).When the repression is weak (K > K c ), the circuit has a single stable state, producing equivalent amounts of A and B (top).Below the critical point, however, this circuit becomes a bistable toggle switch, either producing A and inhibiting B or vice versa (bottom).This transition corresponds to the bifurcation, from one to two stable points, in the phase diagram of the system (Fig 2c).While this phase diagram (green) is known from previous simulations, no analytical relationship was found, particularly for K c , the critical point.
Here we have modeled this system using Uncoupled Max Cal (Eq. 6) to find accurate ( + F q 0 5 J 5 s 5 h j 9 w P n 8 A 5 g W O p w = = < / l a t e x i t > K i,k (t, t) C Eq. C12):

VI. LEARNING THE DYNAMICAL WIRING OF A NETWORK OF NEURONS
Here we consider a brain-like neural network problem to illustrate how Linear Coupling Max Cal can infer a large network from limited information.Consider a network of N neurons (N = 40 here).Taking the stochastic signals from the neurons, we want to infer the neuronal connectivities, and activity patterns.We illustrate how Linear Coupling Max Cal handles this even when we don't measure signals of all of the neurons together.
At any given time t, a neuron either fires (+1) or is silent (−1) in a time interval δt.The state of each neuron i, v i (t), is dependent on both the present and past states of other connected neurons.We assume only limited information is available, namely the mean values and correlations, as in Eq. 1.The probabilities of different activity patterns {v 1 (t), v 2 (t), . . ., v N (t)} are computed using Eq.
3).This model resembles an Ising model of heterogeneous agents, which is often found effective in capturing observed neural activity [11,31,32].In the context of neural activity, h i (t) (bias) controls the probability that neuron i fires, while K ij (t, s) (connection strength) controls the probability that two neurons (i and j) fire together.The challenge here for learning the dynamics is the large number of neurons [11,33].
We test our predictions against a biologically plausible Ground Truth model of this network [11,31] using time-independent Lagrange multipliers h i (t) = h 0 and K ij (t, s) = K ij (τ ), with τ = |t − s| (Fig. 3; see Methods B for the parameters of the model).h 0 0 was chosen to reflect the tendency of real neurons towards silence, while K ij (τ ) was chosen from a normal distribution to reflect the heterogeneity between neurons [31].A realistic assumption is that for τ > 3, K ij (τ ) ≈ 0 [11].In addition, although real neurons are usually silent, occasionally random firing of a few neurons can trigger a large cascade, or "avalanche" of activity [34].These events can only occur when the wiring strengths between neurons (here our Ground Truth model) are tuned near a critical point, where the wiring strengths are weak enough to allow spontaneous neural activity but strong enough to force other neurons to entrain [31,35].
Linear Coupling Max Cal (Eqs.6 and 7) recovers accurately the key features of neural activity present in the Ground Truth model (Fig. 4).It requires input of only the means and correlations between the neurons (Eq.1).In sharp contrast to the uncoupled approximation K ij (τ ) = 0, Linear Coupling Max Cal correctly recovers the dynamical connections between neurons (Fig. 4a).We then took all three of these models and simulated (see Methods B) how average activity, or neural synchrony, s(t) = i v i (t)/N (N = 40) evolved over time [36].In particular, the Linearly Coupled model correctly captures neural avalanches, where s suddenly spikes and many neurons simultaneously fire, whereas the Uncoupled model does not (Fig. 4b).It also correctly captures the spike frequencies (probabilities of s > 0; Fig. 4c).
Linear Coupling Max Cal is just a first-order approximation, valid in the limit of weak interactions.Here, we also tested how this approximation errors changes as interactions are strengthened.Acting like an inverse temperature β ∼ T −1 , we can modulate the average correlation strength between neurons by multiplying each Lagrange multiplier by β: h i → βh i , K ij → βK ij .When β > 1, connections are stronger; when β < 1, they are weaker.Fig. 4d shows how well Linear Coupling Max Cal captures the features of neural synchrony, P (s), over a wide range of β.As expected, both methods accurately capture the mean s value of synchrony, but only Linear Coupling reasonably captures the fluctuations, or variance V ar(s).In addition, the error is maximal near β = 1 (our original model), suggesting that our method gives reasonable results even in the worst-case (i.e.near critical points).
In principle, higher-order approximations can also be treated, as follows.We could employ the Plefka expansion, which has been fruitfully applied to equilibria [22].Another option would be the Bethe approximation, starting from two-body, rather than one-body terms [25][26][27].Or, mean-field variational inference could be used to constrain arbitrary marginal and joint distributions [23,37,38], rather than means and variances.And deep learning methods could be used to learn higher-order in- teractions [19,39,40].

VII. CONCLUSIONS
We describe here a way to compute the coupling strengths among nodes in a network, given signals taken from the individual nodes.We use an inference principle for dynamics and networks called Maximum Caliber.We give two levels of approximation -Uncoupling and Linear Coupling, which can render computations feasible even for networks that are large or have nonlinearities and feedback.The present approach is simple and computationally efficient.

ACKNOWLEDGMENTS
The research was funded by the WM Keck Foundation (LRMP, KD), the NSF BRAIN Initiative (LRMP, KD: ECCS1533257, ncs-fr 1926781), and the Stony Brook University Laufer Center for Physical and Quantitative Biology (KD).

A. Linear Response Theory
Here we show how to estimate the pairwise interactions, K ij (t, s) using our Linear Coupling approximation.Our approach naturally follows from Linear Response Theory [19].We first recognize that, from the properties of Maximum Caliber distributions, ∂Mi(t) ∂hj (s) = C ij (t, s).
Thus from Eq. 6: Since we already have estimates for our single trajectory Lagrange multipliers, we only use Eq. 9 for the pairwise terms i = j; since our Uncoupled estimates depend on single-trajectory (i = j) terms only, their derivative is 0.
Here the relationship is approximate because we neglect the derivatives of K with respect to M ; we assume that these terms are small, but their inclusion would lead to higher order corrections [41].Due to our Linear Coupling approximation, our couplings are only approximate, K .These results directly imply Eqs.6 and 7.

B. Ground Truth Toy Brain Model
Here we provide additional mathematical details of our method.In particular, we discuss how we chose our Ground Truth, brain model and how we, in practice, back-infer the dynamical couplings between our synthetic network of neurons.We chose our couplings to capture the key properties of the experimental observations described for static [31] and dynamic [11] clusters of real neurons.First, the heterogeneity neural interactions (i = j) can be captured by choosing K ij (τ ) (τ = |t − s|) from a normal distribution mean (K 0 a −τ and standard deviation (K δ a −τ ) [31].For simplicity, we choose K 0 = K δ .Here a > 1 describes the rate that correlations between neurons decay with time.Second, in weakly-interacting systems, such as networks of neurons, self-interactions (i = j) are much stronger than pairinteractions (i = j).Except at τ = 0 (since K ii (0) = 0 for the Ising model), we choose, without loss of gener-ality, K ii (τ ) = 20K 0 a −τ .Third, since neurons have a strong tendency towards silence, we chose h i (t) = h 0 (h 0 < 0) for all neurons.Fourth and most importantly, experimentally observed, neuronal avalanches (see main text) can only occur when pairwise couplings are tuned near a critical point [11]; below this point, neural activity is uncorrelated and random, while below this point it is strongly correlated and perpetually silent.Up to a change of scale, we can choose K 0 = .015and a = 4 for our convenience.To tune our network near criticality, we choose h 0 = −.1;here the very weak correlations between our synthetic neurons (≈ .02 on average) can occasionally sum together to create a cascade of neural activation ("avalanche").

IX. IMPLEMENTATION DETAILS
Here we describe how we computed our Linear Coupling estimates of the Ground Truth Lagrange multipli-ers described for our toy brain example.First, we used a standard Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm (5×10 5 iterations) to compute the means and correlations between our synthetic neurons.From these constraints, we computed our estimates for the pairwise couplings, K ij (τ ) using Eq. 7. When τ ≥ 4, couplings are, on average, greater than 4 4 ≈ 100 fold weaker than at τ = 0 and were safely neglected.Second, each of the uncoupled problems is simply a 4-spin Ising model, constrained by the Ground Truth means and autocorrelations, and was solved exactly for each of our N = 40 synthetic neurons.Finally, Eq. 6 was used to reconstruct the remaining Ground Truth Lagrange multipliers from our Uncoupled estimates.

1 A
< l a t e x i t s h a 1 _ b a s e 6 4 = "I B Q i v 9 g d k D + b 9 u r q E M h G f v I d N O o = " > A A A C W H i c b Z F B a 9 s w F M d l r 1 0 b d 1 u 9 9 t i L a e j I L s F O C 9 u x s M t g l x a a J h C H I C v P i W J Z N t L z I B h / y U I P 2 1 f Z Z X L s l j b p A / H + / J 7 e k / R X l A u u 0 f f / W P a 7 v f 3 3 B 4 c d 5 + j D x 0 / H 7 u e T e 5 0 V i s G Q Z S J T 4 4 h q E F z C E D k K G O c K a B o J G E X J j 7 o + + g 1 K 8 0 z e 4 T q H a U o X k s e c U T R o 5 m a h g B h 7 Y Q Q L L k u q F F 1 X J W O s c n 7 N S s 6 r L 3 V a N S m p w r D G q w a v G p w 8 4 a T B y T N 2 Q p D z d q Y T K r 5 Y 4 t e Z 2 / X 7 / i a 8 X R G 0 o k va u J m 5 D + E 8 Y 0 U K E p m g W k 8 C P 8 e p m Y q t e x i t s h a 1 _ b a s e 6 4 = " c 1 P c n m 5 1

FIG. 2 .
FIG. 2. The toggle switch 2-gene network.(a) A and B inhibit each other.(b) Time trajectories of protein number NA (blue) and NB (yellow).(c) Its phase diagram, from the Uncoupled Max Cal approximation.(green).Solid green: stable phases.Dashed green: unstable phase.Red dots (top to bottom and corresponding to Ground Truth simulations): Supercritical, single state (K > Kc).Critical point (Kc), big fluctuations.Subcritical, toggle switch, two bistable states (K < Kc).(d) The histograms of populations P (N ) in each case, comparing the Grand Truth (blue) to the Uncoupled result (magenta).The red lines correspond to the markers in the phase diagram.Uncoupled Max Cal captures the distribution correctly except at the critical point (center, truncated power law).
Fig 2c, red), analytical relationships for the phase diagram of the toggle switch (See Appendix C Eqs. C10 and C11).In addition, Uncoupled Max Cal captures not only the mean, but the full protein distributions away from the critical point (Fig 2d).And even when this method should fail (Fig 2d, middle), Uncoupled Max Cal allows us to calculate analytically the correct critical point (see Appendix time h j (t) < l a t e x i t s h a 1 _ b a s e 6 4 = " D P h 0 FIG. 3. Simplified neuron wiring diagram.Blue circles are neurons each having bias hj.green and purple edges are connections between neurons (signals separated by a time τ ) with strengths Kij(τ ).

FIG. 4 .
FIG. 4. The Neural network example.(a) The Linear Coupling approximation (K and h ) recovers neuron-neuron connection strengths (Kact) and biases (hact, inset).The Uncoupled approximation would estimate Kij = 0.The black diagonal line represents perfect accuracy (b) Average neural activity (or synchrony, s) from the uncoupled (blue), Linearly Coupled (purple), and true (orange) networks.Like the Ground Truth model, the Linearly Coupled model exhibits avalanches (spikes), an important feature of neural activity.(c) The histogram of s for each model.The Linear Coupling model is much more accurate than the uncoupled approximation alone.(d).Model predictions for different connection strengths (β).(a-c) reflect β = 1, the maximum error.While all methods capture mean activity s , only the Linearly Coupled model captures the fluctuations V ar(s).(See Methods B for the details of our Ground Truth as well as the implementation details).