Nonequilibrium Green’s Functions for Functional Connectivity in the Brain

A theoretical framework describing the set of interactions between neurons in the brain, or functional connectivity, should include dynamical functions representing the propagation of signal from one neuron to another. Green’s functions and response functions are natural candidates for this but, while they are conceptually very useful, they are usually defined only for linear time-translationally invariant systems. The brain, instead, behaves nonlinearly and in a time-dependent way. Here, we use nonequilibrium Green’s functions to describe the time-dependent functional connectivity of a continuous-variable network of neurons. We show how the connectivity is related to the measurable response functions, and provide two illustrative examples via numerical calculations, inspired from Caenorhabditis elegans.

A theoretical framework describing the set of interactions between neurons in the brain, or functional connectivity, should include dynamical functions representing the propagation of signal from one neuron to another. Green's functions and response functions are natural candidates for this but, while they are conceptually very useful, they are usually defined only for linear time-translationally invariant systems. The brain, instead, behaves nonlinearly and in a time-dependent way. Here, we use nonequilibrium Green's functions to describe the time-dependent functional connectivity of a continuous-variable network of neurons. We show how the connectivity is related to the measurable response functions, and provide two illustrative examples via numerical calculations, inspired from C. elegans.

NOTATION
g is for Green's functions, f is for response functions. Lower case is for the direct paths, upper case for the effective ones. Subscript 0 stands for time-translationally invariant. The superscript on G 0 denotes the indices that are excluded from the summation in the definition of G 0 .
Green's functions: • g 0,ij (t − t ): equilibrium direct Green's function for the pair of neurons i ← j • G j 0,ij (t − t ): equilibrium effective Green's function for neurons i ← j in a network • g ij (t, t ) = g 0,ij (t − t ) + π ij [ψ](t, t ): nonequilibrium direct Green's function written as the sum of the equilibrium g 0,ij and a nonequilibrium term π ij that depends on the state ψ of the system.
• G ij (t, t ): nonequilibrium effective Green's functions for neurons in a network Response functions • f 0,ij (t − t ): equilibrium direct response function for an isolated pair of neurons (equal to g 0,ij ) • F 0,ij (t − t ): equilibrium effective response function for two neurons in a network (equal to G j 0,ij ) • f ij (t, t ) = f 0,ij (t − t ) + χ ij [ψ](t, t ): nonequilibrium direct response function for an isolated pair of neurons, written as the sum of the equilibrium response function f 0,ij and a nonequilibrium term χ ij which depends on the state ψ of the system.
• In the linear and time-translationally invariant (TTI) case, the activity ψ i (t) = ψ i,eq + ∆ψ i (t) of neuron i is given by ψ i (t) = ψ i,eq + j∈ all g 0,ij * ∆ψ j (t) + g ext 0,i * I ext i (t), (1) where the ψ eq are the resting activity levels of the neurons (which depend on the equilibrium inputs from the other neurons) and the ∆ψ are the changes of the activity with respect to those resting values. * denotes a convolution such that g * ∆ψ (t) = dt 1 g(t, t 1 )∆ψ(t 1 ) g(t, t ) = g(t − t ) for TTI case (2) g 0,ij (t − t ) is the Green's function, or transfer function, describing the direct interaction i ← j from neuron j to neuron i (the subscript 0 denotes their time-translational invariance). I ext i is an external input to j, and g ext 0,i is the transfer function needed to calculate the activity induced by I ext i . Note that we use ∆ψ to denote the deviation of the current state from the equilibrium activity, and δψ to denote perturbations induced in order to probe the properties of the system (see below).
If we consider an isolated pair of neurons i and j, not connected to any network and with a one-way connection i ← j, then the right-hand side of Eq. (1) for neuron i contains a convolution only for the single neuron j. If we know ∆ψ j , and ∆ψ j entirely sets the boundary conditions of the problem because, e.g. it is the only neuron being perturbed, the ψ i (t) calculated from Eq. (1) is exact. In this situation, j would not have any input but I ext j , and would be ∆ψ j (t) = ψ j,eq + g ext 0,j * I ext j (t) A strategy to probe the dynamical properties of the two-neurons system in a given state ψ = ψ eq + ∆ψ is to measure its response functions. ∆ψ, here, has been produced in an unspecified way: for example by an external perturbation, a sensory stimulus, or a change in the cellular properties that has determined a change in the resting activity of the neuron with respect to the usual equilibrium. To measure the response function, we can induce in j an additional perturbation δψ j via an external current I ext j and measure the response that this produces in i. This response can be expressed as (3) is the response function of the edge (i, j), and in the linear and time-translationally invariant case is equal to the Green's function g 0,ij (t − t ). Moreover, since the system is linear, Eq. (3) is valid for δψs of any magnitude.

Neurons in a network
Eq. (1) only describes direct connections between the neurons. If the pair (i, j) is embedded in a network, the signal can travel from j to i both on direct and indirect paths. Therefore, if the boundary conditions are set by ∆ψ j as above, in this case Eq. (1) yields a result that is correct only up to the first perturbative order. To obtain the exact result, Eq. (1) should be solved for each neuron in the network, and the values of ψ(t) updated at each time-step t for each neuron.
However, the direct Green's functions g 0,ij of a linear system can be condensed in effective, or connected, Green's functions G j 0,ij (t − t ), which in the context of Volterra integral equations are the so-called resolvent kernels [1] (for the meaning of the superscript j see below). If we know ∆ψ j and if we want to calculate the activity ∆ψ i produced in i without calculating ∆ψ for the whole network (and if the boundary conditions are set by ∆ψ j ), we can recursively insert the contribution from all the neurons (and therefore all the paths) in Eq. (1) to obtain such connected Green's function. The goal is to write Focusing on the ∆ψ i in Eq. (1), dropping for the moment the external inputs, we have The ∆ψ µ s are themselves only determined by ∆ψ j , so that we can substitute ∆ψ µ with g 0,µj * ∆ψ j + ν g 0,µν * ∆ψ ν ∆ψ i (t) = g 0,ij * ∆ψ j + µ =j g 0,iµ * g 0,µj * ∆ψ j + µ =j ν g 0,iµ * g 0,µν * ∆ψ ν .
Recursively repeating this step, we obtain G j 0,ij can be read off the above equation comparing it with Eq. (4). The superscript j denotes the fact that j is excluded in the sums.
Additionally, noting that the terms inside the square brackets are G j 0,µj , we obtain an equation relating the direct and connected Green's functions: In a linear system, one would write the above equation in Fourier space, where the convolutions are transformed in simple products, like in Ref. [2]. But this cannot be done when the functions are not time-translationally invariant, as is the case in the main discussion of this work.
Because we already take into account indirect paths via G 0 , in Eq. (4) the convolution is summed only over the neurons that set the boundary conditions. In the simple example in which an external current is injected in a single neuron j, the convolution should be performed only with ∆ψ j . If an external current is injected in two neurons, the convolutions in Eq. (4) should be with the activities produced by those external currents separately.
As in the previous section, we can characterize the system measuring its response functions. This can be done experimentally by inducing a perturbation δψ j (t) and measuring the produced δψ i (t). Also for neurons in a network, in the linear and time-translationally invariant case the response functions F 0,ij are equal to the Green's functions G j 0,ij , such that one can write The activity produced by an external current I ext j in j also contains a network term if the neurons are in a network. With a recursive procedure similar to the one leading to Eq. (8), one obtains that the ∆ψ j induced by I ext j is ∆ψ j (t) = (δ(t, t ) + G 0,jj ) * g ext 0,j * I ext j .
G 0,jj does not have a superscript j, which means that it should be calculated as The next section explains why in some cases some indices have to be excluded from the sums and in others they don't.
Self loops, measured activities, and alternative formulation The goal of writing ∆ψ i = G j 0,ij * ∆ψ j , and δψ i = F 0,ij * δψ j , is to have expressions that relate measured activities. If we know ∆ψ j (t) because we measured it, ∆ψ j already contains all the corrections due to the presence of the network. In the summation in Eq. (8), µ cannot be equal to j because otherwise we would be considering paths that start from j and return to j before coming to i. These paths affect ∆ψ j because they bring input to j, so their contribution has already been accounted for in the known ∆ψ j . Excluding j from the summation in Eq. (8) is equivalent to considering paths that start from j and end on i without ever coming back to j.
When we want to calculate, instead, the activity ∆ψ j produced by an external current injected into j (Eq. (10)), we need to consider those network effects on j itself that we discarded above. This is why in Eq. (10), G 0,jj appears without a superscript j.
Alternative formulation An alternative formulation in which no index is excluded from all G 0,ij s is possible, but comes with a crucial disadvantage. Let us call ∆ψ j (t) the activity that an external stimulation I ext The activity of that same neuron j considering the feedback from the network, would be given by ∆ψ j = (δ(t, t ) + G 0,jj ) * ∆ψ j (the same as Eq. (10)). The activity of other neurons i due to ∆ψ j would now be given by ∆ψ i (t) = G 0,ij * ∆ψ j (t), with the summation in G 0,ij running over all indices. However, the relation now is between the measured ∆ψ i and ∆ψ j , the activity of neuron j if it were isolated from the network, which cannot be measured. To make use of this relation to obtain information about the network, we need to make assumptions about ∆ψ j , which translate into assumptions about g ext 0,j and I ext j . In some cases, g ext and I ext can be estimated or separately characterized, like in the case of optogenetic stimulations with light-gated ion channels, or channelrhodopsins (I ext depends, for example, on the gating kinetics of the channelrhodopsin). Given that the stimulus is an electrical current, the coupling to the membrane potential of a neuron is via a first-order equation and g ext ∝ e −t/τ . This characterization is not general, however, because in the same nervous system stimuli might be applied to the neurons also via other types of molecules, like Gprotein coupled receptors. With these receptors, the coupling to the neuron's activity is likely to be at least via a second-order equation because there are intermediate steps involved, so that g ext is not a simple exponential anymore.
Depending on the context and the goal of the calculations, be they for numerical simulations or data analysis, either of the two formulations (with ∆ψ or ∆ψ) can be used. In the following, we will keep working with the measured ∆ψ.

Isolated pair of neurons
Also in the nonlinear case we would like to formally write the state of the system as a convolution of a Green's function and the past activities of the neurons: (13) In contrast to the linear and time-translationally invariant case, the Green's function g ij in Eq. (13) is not a time-invariant kernel (hence, nonequilibrium). Instead it is, in principle, functionally dependent on the whole state of the system ψ. We can write g ij [ψ] as the sum of a linear and time-translationally invariant term g 0,ij (t − t ) and a nonequilibrium term π ij (t, t ): There is no general recipe for calculating π ij [ψ]. Different models describing the nonlinear or time-dependent nature of the edge (i, j) will involve their own specific solutions. However, Eq. (13) and (14) are useful as formal devices because they will allow us to express how a given π ij [ψ], once it is known, should be used to calculate the modified properties of the rest of the network.

Response function of isolated pair
In the nonequilibrium case, the response function f ij is, in general, not equal to the Green's function g ij . If the system is in the state ψ = ψ eq + ∆ψ (where ∆ψ has been produced in an unspecified way), and we induce a perturbation δψ β in β to probe the system, the activity δψ α produced in α can still be written as but here f αβ is given by where χ αβ [ψ] is the nonequilibrium part of the response function and is given by (17) Also here, f ij is not a time-translationally invariant transfer function, and depends on the state of the system ψ via χ αβ [ψ]. Therefore, it has to be calculated according to its full expressions (or approximations to it). But, provided that δψ j is "small" (see below), once χ αβ and f ij are computed, Eq. (15) can be used to to compute the response δψ α to any perturbation δψ β .
In neurons, we can expect the linear-response regime to be valid also for δψ of significant magnitude. In the most general nonlinear system, linear responses around an arbitrary state can be very limited approximations. But nonlinearities like sigmoidal thresholds at the synapses have large activity intervals where they are well-approximated by a linear function. In the examples in the main text, at equilibrium the synapse (α, β) "sits" at the bottom of φ(V β ), the sigmoidal function describing the dependence of the synaptic vesicle release on the calcium influx into the presynaptic site. As the system is brought out of equilibrium, the synapse reaches the central region of the sigmoid φ(V β ), where the sigmoid is again almost linear and where the linear-response approximation works well.
Two cases can be discussed starting from Eq. (17). π(t, t ) can be nonequilibrium in two ways: it can be time-dependent or explicitly nonlinear (or a combination of the two). If π is purely time-dependent and does not depend on ψ β , the only nonzero term in the right-hand side of Eq. (17) is π. In this case, the response function is equal to the Green's function: f αβ = g 0,αβ + π αβ . If π depends on δψ β , then also the second term is nonzero and the f αβ = g αβ .

Neurons in a network: nonequilibrium Green's functions
The connected Green's functions for neurons in a network will be determined also by the nonequilibrium terms. To derive an equation for the nonequilibrium Green's functions G ij we start from the following experimental consideration related to C. elegans that can be taken as a general starting point. Characterizations of some synpases in the worm have shown that they are linear throughout a large part of the physiological range of membrane potentials [3][4][5]. However, we know that nonlinearities and time-dependences are critically important in the C. elegans nervous system and in nervous systems generally, because they allow the network to perform computations, including for example responding to sensory stimuli in a context dependent manner [6,7].
How does a network have many linear edges but also show widespread nonlinear behaviors? In the integral formulation with nonequilibrium Green's functions it is straightforward to show how these two observations can coexist.
We, therefore, consider a network in which only one of the edges, (α, β), displays a significant nonlinearity. This is in contrast to an approach in [8] which assumes nonlinearities that are homogeneous over the network and then proceeds with their systematic expansion.
The connected nonequilibrium Green's function for the effective edge (i, j) can be calculated with a recursive procedure similar to the one leading from Eq. (1) to Eqs. (4) and (8), but using the nonequilibrium Green's function in Eq. (14) for the edge (α, β). The result for i = α, j = β is where (A * B)(t, t ) = dt 1 A(t, t 1 )B(t 1 , t) and the activity ∆ψ i can be written as with, as usual, a summation on the ∆ψs that set the boundary conditions (b.c.). Differently from Eq. (8), here G ij is not a time-translationally invariant transfer function. The nonequilibrium π αβ [ψ] and G ij have to be calculated for each time-step taking the nonlinear and time-dependent nature of π αβ [ψ] into account, and the calculation depends, therefore, on the specific state ψ considered. Once it is calculated, however, it can be plugged in in the expression for all the nonequilibrium G ij (Eq. (18)), and can be used to calculate ∆ψ i via convolutions only, without running a nonlinear calculation.
The full nonlinear calculation of π αβ [ψ] does not involve, however, all the neurons of the network, in common situations. π αβ [ψ] will depend likely on a limited subsets of neurons only. A nonlinearity determined by threshold or saturation at a synapse will depend only on the presynaptic activity ψ β . If neuromodulation is involved, it will be determined by a defined set of neurons.
Additionally, in the calculation of G ij the effect of the network recurrence is condensed in only a single additional effective edge (β, j). To calculate G ij one needs, therefore, to simultaneously calculate only G β,j .
There are some special cases for Eq. (18).
If i = α and j = α If i = j and i = α, β Neurons in a network: nonequilibrium response functions The response function of an effective edge (i, j) is modified by the network in a similar way as the nonequilibrium G ij in Eq. (18). We assume again one single nonequilibrium edge (α, β) and insert recursively the contributions from all the neurons in Eq. (3), using Eq. (16) for the edge (α, β). The result for i = α, j = β is The response δψ i in i produced by a small perturbation δψ j in j can be written as Also here F ij is not a time-translationally invariant transfer function, and depends on the state of the system ψ via χ αβ [ψ]. But, provided that δψ j is "small" (see discussion in the section about the isolated-pair response function), once χ αβ and F ij are computed in the full nonlinear calculation, Eq. (28) can be used to to compute the response δψ i to any perturbation δψ j . Also for the response functions F ij , as for the Green's functions G ij , there are some special cases of Eq. (27).
If i = α, β and j = β If i = α and j = α If i = j and i = α, β Experimental characterization F ij are the response functions that can be obtained in experiments on networks of neuron. Importantly, they are always well defined, for two reasons. First, regardless of how nonlinear the system is, an experiment can always be designed to probe the system by perturbing it and measuring its response. Second, their role is always clear even if we are not considering any other neuron in addition to i and j, both in an experiment and in theoretical calculations, even if in the latter the response function would be difficult to calculate without making assumptions on the rest of the network.
In an experiment, F ij can be measured as responses to impulsive perturbations. One could object that, by applying an impulsive external stimulation I ext j to neuron j, the δψ j induced in j would be given by Eq. (10) and would not be itself a delta function. Therefore, the response of neuron i would be given not by F ij alone, but by the convolution of (1 + G 0,jj ) * g ext 0,j with F ij . However, if the coupling between the external stimulation I ext j and ψ j is via a first-order differential equation, e.g. for current in an equation for membrane potential, for which g ext 0,j ∝ e −(t−t )/τ , this limit is of no concern. The timescale τ is a property of neuron j and applies to any input to that neuron, so that there are no response functions with a path through j that are "faster" than g ext j .

RELATION OF NONEQUILIBRIUM GREEN'S FUNCTIONS TO OTHER METHODS
Eqs. (16) and (27) (corresponding to Eqs. 2 and 4 of the main text) describe an approximately linear regime on top of an arbitrary state of the system. Switching linear dynamical systems (SLDS) models [9,10] assume that such locally linear regimes exist. They describe the dynamics of a nonlinear system as a temporal sequence of linear systems each with different parameters, and have been applied to C. elegans. In existing SLDS models, the time-dependent switching between parameters is entirely phenomenological. In our approach, Eqs. (16) (main text Eq. 2), (27) (main text Eq. 4), and (18), explicitly govern how nonlinearities in the network produce timedependent changes to a linear system. Here the response functions contain the time-dependent parameters of the SLDS.
To draw a parallel between the numerical examples in the main text and SLDS [9,10], the state induced by the perturbation I µ corresponds to switching from the equilibrium linear system to another one with different parameters.

Equilibrium Green's functions
the equilibrium Green's functions σ 0,ij , g g 0,ij , and g s where θ(t) is the Heaviside step function.

Nonequilibrium Green's functions and nonlinear terms
There are three sources of nonlinearities that can be added back to the above linear and TTI Green's functions. The first is the saturation of the postsynaptic current when the postsynaptic membrane potential approaches the reversal potential of the ionotropic receptor through which the current flows (Eq. 5 of the main text). This nonlinearity changes g s ij . The second is the saturation of the synaptic activity (which can range between 0 and 1) due to the finite number of receptors, and the third is the sigmoidal dependence φ(V j ) of the vescicle release on the presynaptic potential (Eq. 6 of the main text). These two nonlinearities change σ ij and are the ones considered in the numerical examples.
Derivation of NEGF To derive the explicit expression for the nonequilibrium Green's functions we can proceed as follows. As a reminder, in Eq. 6 of the main text electrical synapses linearly couple V i ← V j via the equilibrium Green's function g g ij in Eq. (39). Chemical synapses, instead, couple V i ← V j via an intermediate synaptic activity variable s ij , governed by Eq. 7 (main text). Both the steps V i ← s ij and s ij ← V j involve nonlinear terms, so we derive two separate nonequilibrium Green's functions for the two steps that can be then combined, σ ij (t, t ) and g s ij (t, t ), that account for different nonlinearities as described in the previous paragraph.
In both cases, we consider the equations for the deviations ∆s ij and ∆V i from the equilibrium values, and then write the equations as On the left-hand side, x is either V i or s ij and D t is a linear differential operator that contains all the linear terms of the equations. On the right hand-side, A(t) is the sum of a linear forcing term ∆y and all the nonlinear terms of the differential equation. ∆y(t) corresponds to ∆s ij for the step ∆V i ← ∆s ij , and to ∆V j for ∆s ij ← ∆V j . If there are no nonlinear terms, then we have an equilibrium Green's function g 0 that is obtained from Given the Green's function g 0 , we can write ∆x(t) = g 0 * A)(t). We can write the latter expression for ∆x also if A(t) is a modified forcing term containing also nonlinear terms. In the nonlinear case, we should read off the nonequilibrium Green's function g(t, t ) from the expression of ∆x(t), using the definition ∆x(t) ≡ g * ∆y)(t) = dt 1 g(t, t 1 )∆y(t 1 ).
While this reading-off step is just a regrouping of terms, it yields the nonequilibrium Green's function, which is a useful formal device to compute other properties of the network, like the response functions. The nonlinear differential equation for ∆s ij is (from Eq. 7 of the main text) By linearizing it around equilibrium we obtain the Green's function σ 0,ij in Eq. (37), while using the above procedure and the definition ∆s ij (t) ≡ dt 1 σ ij (t, t 1 )∆V j (t 1 ), we obtain the expression for the nonlinear σ ij (t, t ) (Eq. 8 in the main text).
The nonlinear equation for ∆V i is (from Eq. 6 of the main text) Again, its linearization leads to the equilibrium Green's functions in Eqs. (39) and (40). Focusing on one edge ij to obtain the nonequilibrium Green's function, we get ∆V i (t) = (g g 0,ij * ∆V j )(t) From the above equation we can read off the nonequilibrium Green's function g ij in Eq. 9 of the main text. Finally, the nonequilibrium response function χ αβ (t, t ) in Eq. 10 of the main text is obtained applying Eqs. 3 and 4 (main text) to the nonequilibrium Green's function σ αβ (t, t ):

Additional example and steps for calculations
We consider a more complex network with recurrence ( Fig. 1a) and provide a detailed description of the steps involved in calculating that network's response functions. We consider the nonequilibrium state of the network induced by a 0.5 pA current pulse I β (t) injected into neuron β, which starts at t = 0 and lasts 1 s (the network parameters are described in the next section).
The response functions are calculated in two main steps. First, the nonequilibrium response function f αβ of the nonlinear edge α ⇐ β is calculated. This part of the calculation explicitly involves the nonlinearities, but it only involves the nonlinear edge α ⇐ β, and the effective feedback edges β ← α and β ← β. Second, the nonequilibrium response function of any other effective edge is calculated via simple convolutions, without the need to re-run the nonlinear calculation.
In the first part we must compute χ αβ (t, t ). For our C. elegans specific model, χ αβ (t, t ) can be obtained using main text Eq. 3 and considering that f αβ = g s αβ * χ αβ , where χ αβ is given by main text Eq. 10. To compute χ αβ , we need to know ∆V β , which can be calculated by convolving the input current I β (t) with the nonequilibrium Green's function G ext β (t, t ), which needs to be included in the explicitly nonlinear calculation using Eq. (24).
In the second part of the calculation, we use χ αβ (t, t ) to compute any other response function of the network via Eq. (27) (Eq. 5 of the main text). This step does not involve any explicit knowledge of the form of the nonlinearities, because they are already captured by χ αβ (t, t ), computed in the first step. As an example, we compute  is calculated for neurons µ and ν (blue). b Nonequilibrium and equilibrium F (0,)νµ (t, t ) for selected times t (colored and black curves, respectively). c δVν←µ induced by Iµ,1 (red) and Iµ,2 (green), calculated as ∆Vν I β +Iµ − ∆Vν I β (solid lines) and Fνµ * δVµ Iµ (red crosses; Since Iµ,2 = 3000Iµ,1, the responses to the two currents calculated via Fνµ are simply proportional). δVν←µ calculated with the equilibrium response function F0,νµ(t − t ) (blue). inset ∆V β produced by I β only (black), by I β plus Iµ,1 (red dots, almost perfectly overlapping with the black line), and I β plus Iµ,2 (green curve). The gray curve is Iµ,i (scale not shown).
the nonequilibrium response function F νµ (t, t ), plotted in Fig. 1b for selected times (colored lines) together with the equilibrium F 0,νµ (t−t ) (black line). To compute any response function of an edge (i, j) that does not end on β, we need to know F βj (t, t ), the last factor of Eq. (27) (Eq. 5 of the main text). To compute F βj (t, t ), however, we do not need to know any other response function beyond f αβ = f 0,αβ + χ αβ (because the left-hand side and the last factor of Eq. (27) are equal). Therefore, knowing χ αβ (t, t ), we can compute any nonequilibrium F νµ with just the additional step of computing F βµ . With the response functions, we can then compute the response of the system to arbitrary additional stimuli, for example a current I µ (t) injected into neuron µ. In general, this is possible only with additional stimuli that are sufficiently small, because it corresponds to a linearization of the system. However, as discussed above, because of the typical sigmoidal nonlinearities we find in neural networks, we can expect the nonequilibrium linearization to work over stimuli with a large range of amplitudes. We show this in Fig. 1c.
To check that the response function yields the correct result, we plot the responses of the system to a small, 0.05 s current pulse I µ,1 , injected in neuron µ at time 0.5 s in addition to the current I β . This current pulse into µ produces a very small additional change of the membrane potential V β of the downstream neuron β, as can be seen in the inset of Fig. 1c. There the red dotted curve, computed including I µ,1 , overlaps almost perfectly with the black curve, computed without I µ,1 . Therefore, we expect the linear approximation to work well. We can check this by comparing the response δV ν←µ of neuron ν to I µ,1 computed using the nonequilibrium response function (red crosses in Fig. 1c) with the one computed including I µ,1 in the explicitly nonlinear calculation (solid red line). As expected, they overlap very well.
The linear regime captured by the nonequilibrium response function holds also with large stimuli, in certain cases. To show this, we inject a much larger current I µ,2 into µ, which produces a downstream change in V β that is of the same order as the one produced by the main current pulse I β (green and black curves, respectively, in the inset of Fig. 1c). In such condition, it is not a given that the system can be linearized. Nevertheless, the response induced in ν by I µ,2 is still largely captured by the response function. In fact, δV ν←µ calculated including I µ,2 in the explicit nonlinear calculation (Fig. 1c green curve), is very close to δV ν←µ calculated via the response function (red curve, which, except for the appropriate rescaling by the ratio of I µ,2 and I µ,1 , is the same as the response induced by I µ,1 ). As a comparison, the blue curve is the response δV ν←µ calculated with the equilibrium response function, which produces a qualitatively wrong result.  the synapses, γ s ij = 10 S/F, E syn,ij = 0 mV, a r,ij = 5, a d,ij = 5, β ij = 125 V −1 . V th,ij = V eq,j for all the edges as in [12], except for edge (α, β), for which we set V th,αβ = −10 mV. For the second example, the additional inhibitory edge 3 1 has γ s 31 = 20 S/F, E syn,31 = −70 mV, and a r,ij = 1, while all the other parameters of the synapse are equal to the ones above.
All the convolutions are performed using the Gregory integration scheme (of order 8), adapted from Ref. [13].