# Mutual multilinearity of nonequilibrium network currents

Sara Dal Cengio, Pedro E Harunari, Vivien Lecomte, Matteo Polettini

► To cite this version:

Sara Dal Cengio, Pedro E Harunari, Vivien Lecomte, Matteo Polettini. Mutual multilinearity of nonequilibrium network currents. SciPost Physics, 2025, 19 (4), pp.111. 10.21468/SciPost-Phys.19.4.111. [https://arxiv.org/abs/2004.04106]

HAL Id: hal-05410632

https://hal.science/hal-05410632v1

Submitted on 11 Dec 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mutual multilinearity of nonequilibrium network currents

Sara Dal Cengio¹,²*, Pedro E. Harunari³,⁴†, Vivien Lecomte² and Matteo Polettini

1 Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA

2 Université Grenoble Alpes, CNRS, LIPhy, FR-38000 Grenoble, France

3 Aix Marseille Université, CNRS, CINAM, Turing Center
for Living Systems, 13288 Marseille, France

4 Complex Systems and Statistical Mechanics, Department of Physics and Materials Science,
University of Luxembourg, 30 Avenue des Hauts-Fourneaux,
L-4362 Esch-sur-Alzette, Luxembourg

∗ saradc@mit.edu, † pedroharunari@gmail.com

## Abstract

Continuous-time Markov chains have been successful in modelling systems across numerous fields, with currents being fundamental entities that describe the flows of energy, particles, individuals, chemical species, information, or other quantities. They apply to systems described by agents transitioning between vertices along the edges of a network (at some rate in each direction). It has recently been shown by the authors that, at stationarity, a hidden linearity exists between currents that flow along edges: if one controls the current of a specific “input” edge (by tuning transition rates along it), any other current is a linear-affine function of the input current [Phys. Rev. Lett. 133, 047401 (2024)]. In this paper, we extend this result to the situation where one controls the currents of several edges, and prove that other currents are in linear-affine relation with the input ones. Two proofs with distinct insights are provided: the first relies on Kirchhoff’s current law and reduces the input set inductively through graph analysis, while the second utilizes the resolvent approach via a Laplace transform in time. We obtain explicit expressions for the current-to-current susceptibilities, which allow one to map current dependencies through the network. We also verify from our expression that Kirchhoff’s current law is recovered as a limiting case of our mutual linearity. Last, we uncover that susceptibilities can be obtained from fluctuations when the reference system is originally at equilibrium.

<figure><img src="image_2.png" /></figure>

Copyright S. Dal Cengio *et al*.

This work is licensed under the Creative Commons Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 2025-03-05

Accepted 2025-09-23

Published 2025-10-27

doi:10.21468/SciPostPhys.19.4.111

<figure><img src="image_3.png" /></figure>

Check for
updates

## Contents

1 **Introduction** 2

2 **Setup and statement of the main result** 5

2.1 **Illustration** 7
<table><tbody><tr><td>3</td><td>Graph-theoretic steady-state proof</td><td>8</td></tr><tr><td>3.1</td><td>Proof of inductive seed</td><td>8</td></tr><tr><td>3.2</td><td>Proof of inductive step</td><td>10</td></tr><tr><td>4</td><td>Laplace-space finite-time proof</td><td>11</td></tr><tr><td>4.1</td><td>Settings</td><td>12</td></tr><tr><td>4.2</td><td>Proof of mutual multilinearity</td><td>12</td></tr><tr><td>4.3</td><td>Taking the limit &zeta; &rarr; 0 and mapping to the one-input-current case</td><td>14</td></tr><tr><td>5</td><td>Further developments</td><td>15</td></tr><tr><td>5.1</td><td>Recovering Kirchhoff's current law</td><td>15</td></tr><tr><td>5.2</td><td>Connection to equilibrium fluctuations</td><td>16</td></tr><tr><td>6</td><td>Conclusion and perspectives</td><td>18</td></tr><tr><td>A</td><td>Construction of potentials and affinities from transition rates</td><td>19</td></tr><tr><td>B</td><td>Linear independence of admissible currents</td><td>19</td></tr><tr><td>C</td><td>Proof that 1 &minus; &lambda;<sup>I<sub>1</sub></sup><sub>e<sub>2</sub>&larr;e<sub>1</sub></sub>(r<sub>E\I<sub>1</sub></sub>) &lambda;<sup>I<sub>2</sub></sup><sub>e<sub>1</sub>&larr;e<sub>2</sub></sub>(r<sub>E\I<sub>2</sub></sub>) cannot cancel</td><td>20</td></tr><tr><td>D</td><td>A fluctuation-dissipation relation</td><td>21</td></tr><tr><td></td><td>References</td><td>22</td></tr></tbody></table>

# 1 Introduction

In nonequilibrium Markov chains, the net number of times an edge is traversed is called a current, constituting a key observable. Close to equilibrium, currents (of any kind) flow proportionally to conjugated forces and therefore respond linearly to one another. This is the tenet of linear thermodynamics, which grounds on the notion of conjugated observables [1] and the varied fluctuation-dissipation relations that come with it [2–4]. Far from equilibrium, currents are generally nonlinear functions of the conjugated forces and therefore nonlinear functions of one another, and fluctuation-dissipation relations break down.

Historically, only a few relations are known to hold arbitrarily far from equilibrium, such as the fluctuation theorems [5–9], the thermodynamic uncertainty relation [10, 11], and extensions of the fluctuation-dissipation relation [12–16]. Recently, new results on the response of observables have started to surface [17–29], which include fluctuation-response relations [30, 31], the connection between responses and correlations of dynamical events [32], and the mutual linearity of currents [33]. In this latter study [33], we have shown that any two stationary currents always satisfy a linear-affine relation between each other, even arbitrarily far from equilibrium, when the transition rates of a single edge are controlled. It has been conjectured that a similar property might hold when more edges are controlled [34]. In the present work, we prove that a linear-affine relation among stationary currents can indeed be established even when the transition rates along multiple edges are simultaneously controlled. We explore the details and conditions for its validity, as well as the recovery of Kirchhoff’s current law and connections to linear response theory.
<figure><img src="image_4.png" /><figcaption>Figure 1: Illustration of central graph-theoretical notions. Starting from a connected graph (left panel), a set of edges is identified in red (middle panel) whose removal results in a spanning tree (right panel), a connected graph without cycles that contains all vertices. Adding back one chord to the spanning tree creates exactly one cycle.</figcaption></figure>

A special case of linear relation among stationary currents holds in the form of the Kirchhoff current law (KCL), which is central in Schnakenberg's geometrical theory of Markov chains [35-37]. Therein, a continuous-time Markov chain is represented as a network (see Fig. 1), where each vertex corresponds to a state of the chain and an edge is drawn between two vertices whenever a transition between the corresponding states is possible. Using geometrical [36] or algebraic tools [38], it is possible to identify a subset of edges whose removal from the network results in a spanning tree, i.e. a tree that connects all vertices (see Fig. 1). The subset of edges that defines a spanning tree is not unique; it is referred to as a fundamental set, and once it is fixed, each element of this set is called a chord. The number of chords in a fundamental set is fixed by the network's topology and is known as the cyclomatic number $n_c$. Notably, for a given fundamental set, each chord can be uniquely associated with a specific cycle in the network.

As a direct consequence of Kirchhoff current law, any stationary edge current $J_e$ can be expressed as a superposition of the stationary currents $J_i$ flowing along the chords of a fundamental set, as:

$$J_e = \sum_{i=1}^{n_c} \gamma_{e,i} J_i, \qquad (1)$$

where $\gamma_{e,i} = \pm 1, 0$ are algebraic coefficients specifying whether edge $e$ belongs to cycle $i$ and the relative orientation within it. Eq. (1) is the result of the conservation of probability, which imposes constraints on the value that currents can take at steady-state. It tells that if an observer were to control or fix the value of all $n_c$ chord currents, all other currents would be determined by these values through the linear topological constraints of Eq. (1). From a thermodynamic viewpoint, the chord currents $J_i$'s can be viewed as the independent degrees of freedom at steady state associated with thermodynamic forces (or affinities) $a_i$ which contribute independently to the stationary entropy production. As such, Eq. (1) underpins a complete thermodynamic description of stationary currents, which have proved valuable both from a practical and a conceptual viewpoint [38-44].

Recently, a different form of linearity among stationary currents has been established by us in Ref. [33]. Therein, we ask the question of how stationary currents respond to manipulations of the transition rates of a single edge current, which we term input current. We proved that any current of the Markov process $J_e$ can be expressed as a linear-affine function of a single input current $J_i$, as:

$$J_e = J_e^i + \lambda_{e \leftarrow i}^i J_i. \qquad (2)$$
<figure><img src="image_5.png" /><figcaption>Figure 2: Summarizing two types of mutual linearity among nonequilibrium stationary currents. (Left) Any stationary current $J_e$ (blue) in the graph can be represented as a linear-affine function of an “input” current $J_i$ (red) when the input edge $i$ is controlled; the scalar coefficients $J_e^{\leftarrow i}$ and $\lambda_{e \leftarrow i}^i$ depend both on dynamical and topological properties of the Markov process. This result expresses the mutual linearity of currents, derived in Ref. [33]. (Right) Kirchhoff’s current law in its global form provided by Eq. (6): A fundamental set of edges (red) whose removal leaves a connected graph without cycles (a spanning tree) is fixed. The stationary current of any other edge (blue) can always be represented as a linear combination of the currents along the red edges (including when these edges are controlled). These choices of edges are not unique, and the coefficients $\gamma_{e,i}$ are either 0 or $\pm 1$ depend solely on the graph topology. In this paper, we address the intermediate situation in which the number of controlled edges is greater than one but less than the cyclomatic number $n_c$.</figcaption></figure>

In contrast to KCL, here the scalar coefficients $J_e^{\leftarrow i}$ and $\lambda_{e \leftarrow i}^i$ are functions on the transition rates of the continuous-time Markov chain, excluding the transition rates of the input current itself. As such, they depend both on topological and dynamical properties of the Markov process. In particular, $J_e^{\leftarrow i}$ is the value of the current $J_e$ in the process where the input edge is removed. Interestingly, the dependencies of $J_e^{\leftarrow i}$ and $\lambda_{e \leftarrow i}^i$ on the transition rates can be made explicit using the technique of spanning tree polynomials [33, 45]. For simplicity, here and in the rest of the paper, we use the term linearity for linear-affine relations like Eq. (2). The affine coefficient is assumed to be present unless otherwise stated.

These two forms of mutual linearity among currents, both valid far from equilibrium, are utterly different in nature. On one hand, KCL states that any stationary current can be expressed as a linear combination of a full set of chord currents, with simple coefficients $\in \{-1, 0, +1\}$, and is the consequence of probability conservation. On the other hand, Eq. (2) reveals a linear-affine relation between any stationary current and one input current, when only the rates of the latter are controlled.

In the present paper, we address the question of what is the fate of the linearity in Eq. (2) when several transitions rates, belonging to different edges, are simultaneously controlled. We demonstrate that linearity persists even when multiple $n \le n_c$ chord currents are controlled, and obtain explicit expressions for the current-current susceptibilities. We check that, as expected, such expressions yield back KCL, Eq. (1), in the case $n = n_c$. In this regard, Eqs. (1)-(2) can be seen as two limiting cases of a more general form of mutual multilinearity (MML) governing the system’s behavior to changes in transition rates (see Fig. 2). We will provide two different proofs, one based on a direct graph-theoretic analysis and one based on the Laplace
transform analysis. The strategy for the graph-theoretic proof is a “top-down” induction with a seed treating the case of one input edge less than the cyclomatic number, and a step treating the case of one input edge less than in the previous step. The above result holds for arbitrary finite graph geometry, arbitrary choices of transition rates and arbitrary variations of the input rates.

The plan of the paper is as follows. In Sec. 2 we introduce the framework of Markov processes and state our main result. Sections 3 and 4 are devoted to two different proofs, which provide different insights. The graph-theoretical proof explicitly employs the Kirchhoff current law, Eq. (1), while the Laplace transform proof allows for a generalization at finite time, and provides connection with the notion of current-to-current susceptibility [33]. Finally, in Sec. 5, we show how to recover KCL as a limiting case from MML, and make connection with the framework of linear response by introducing the notion of affinities. Remarkably, we show that the transport coefficients commonly defined close to equilibrium in fact describe the current-current relation arbitrarily far from equilibrium, provided that the reference transition rates fulfill detailed balance. Finally, in Sec. 6 we draw our conclusions and outline future perspectives.

*Remark:* Note that during the writing of this work, F. Khodabandehlou and coworkers have obtained analogous results for multigraphs [46] using different methods; our contribution sheds light on the conditions ensuring that the susceptibilities are well-defined and unique (and on the time-dependent case, in Laplace space).

## 2 Setup and statement of the main result

Consider a finite connected graph $\mathcal{G}$ with vertices $x \in \mathcal{X}$ and edges $e \in \mathcal{E}$. We assume no multiple edges between two vertices. To each edge we assign a reference forward orientation $+e$ (or simply $e$) from a source to a target vertex. The incidence matrix $\mathbf{S}$ has entries $\mathbf{S}_{xe}$ equal to $-1$ if $x$ is the source of $e$, $+1$ if $x$ is the target of $e$ and $0$ otherwise, for $x \in \mathcal{X}$ and $e \in \mathcal{E}$. We denote $-e$ the edge with opposite orientation and $\tilde{\mathcal{E}}$ the set of edges oriented backwards. We assign positive weights $\mathbf{r}: \mathcal{E} \cup \tilde{\mathcal{E}} \to \mathbb{R}_+$ to forward and backward directions of an edge. For notational simplicity we now define $\tilde{\mathcal{E}} = \mathcal{E} \cup \tilde{\mathcal{E}}$. Let $\mathbf{r}_{\mathcal{E}}$ be the rates of a subset $\mathcal{A} \cup \tilde{\mathcal{A}} \subseteq \tilde{\mathcal{E}}$.

Consider the continuous-time Markov chain having $\mathbf{r}$ as transition rates. Let $\mathbf{R}$ be the rate matrix with entries $\mathbf{R}_{xy} = r_{x \leftarrow y} - \delta_{x,y} \sum_{z \in \mathcal{X}} r_{z \leftarrow x}$. The probability $\mathbf{p}(t) = (p_x(t))_{x \in \mathcal{X}}$ of being at given vertices at time $t$, given an initial probability $\mathbf{p}(0)$, evolves according to the master equation

$$ \dot{\mathbf{p}}(t) = \mathbf{R} \mathbf{p}(t). \qquad (3) $$

The absence of source or sink terms in Eq. (3) ensures probability conservation. Let $\pi = (\pi_x)_{x \in \mathcal{X}}$ be the normalized stationary distribution satisfying $\mathbf{R}\pi = 0$, which by the above assumptions is unique. For simplicity, we also assume $r_e \neq 0 \Leftrightarrow r_{-e} \neq 0$ (weak reversibility).

The main object in our study are the stationary currents, which for the oriented edge $e = x \leftarrow y$ are given by

$$ J_e = J_{x \leftarrow y} = r_{x \leftarrow y} \pi_y - r_{y \leftarrow x} \pi_x. \qquad (4) $$

A straightforward consequence of the master equation is Kirchhoff's current law (KCL),

$$ \sum_{e \in \mathcal{E}} \mathbf{S}_{x,e} J_e = 0, \quad \forall x, \qquad (5) $$

stating that the in-and-out stationary currents balance at any vertex $x$, locally. As we now detail, there also exists a global version of KCL, equivalent to Eq. (5), expressed by means of
spanning trees, a central notion of graph theory. A spanning tree is a minimal set of edges that connects all vertices (and thus contains no cycles). A consequence of Eq. (5) is that for any spanning tree $\mathcal{T}$ there exist unique coefficients $\gamma_{e,i} \in \{0, -1, +1\}$, not dependent on the rates, such that<fnref n="1" />

$$J_e = \sum_{i \in \mathcal{E} \setminus \mathcal{T}} \gamma_{e,i} J_i, \qquad (6)$$

with the sum running over the set of edges not belonging to $\mathcal{T}$. Notice this is equivalent to Eq. (1). These edges are known as chords and they form what is known as the fundamental set. Its cardinality $n_c$, the cyclomatic number, is the dimension of the kernel of the incidence matrix $S$. The vectors $\gamma_i = (\gamma_{e,i})_{e \in \mathcal{E}}$ are null vectors of the incidence matrix, and they can be represented as simple oriented cycles [36]. In the following, we refer to Eq. (6) as KCL, as it directly follows from Kirchhoff's current law and is sometimes used as an alternative representation of it in the literature.

In this work, we generalize Eq. (6) to the case where the $n < n_c$ currents are a subset of some fundamental set, a condition that we term *admissibility* (see Fig. 3). These $n$ currents $J_i$ are termed *input currents*, while any other current $J_e$ in the graph can be regarded as *output*.

Now consider $n$ input currents along admissible edges $\mathcal{I} = (i_1, \dots, i_n)$. Our main result is the following. For every edge $e \in \mathcal{E}$ there exist, and are uniquely defined, a coefficient $J_e^{\mathcal{I}}$ and $n$ coefficients $\lambda_{e \leftarrow i}^{\mathcal{I}}$ (for $i \in \mathcal{I}$), independent of the $2n$ rates in $r_{\mathcal{I}}$, such that when varying such rates

$$J_e(r_{\mathcal{I}}) = J_e^{\mathcal{I}} + \sum_{i \in \mathcal{I}} \lambda_{e \leftarrow i}^{\mathcal{I}} J_i(r_{\mathcal{I}}). \qquad (7)$$

Notice that in general the above coefficients depend on the rest of the rates $r_{\mathcal{E} \setminus \mathcal{I}}$.

<figure><img src="image_6.png" /><figcaption>Figure 3: Illustration of the admissibility concept. Top: The removal of the input edges results in a connected graph, satisfying the criterion for admissibility. Bottom: The resulting graph is disconnected thus the input set is not admissible.</figcaption></figure>

<footnote n="1">¹Eq. (5) tells that the vector $J = (J_e)_{e \in \mathcal{E}}$ belongs to ker $S$, the right nullspace of $S$; but the cycles associated to $\mathcal{T}$, represented as vectors $(\gamma_i)_{i \in \mathcal{E} \setminus \mathcal{T}}$ of components $(\gamma_i)_e = \gamma_{e,i}$ form a basis of ker $S$; it is then easy to check that the component of $J$ along $\gamma_i$ in this basis is given by $J_i$, thus yielding Eq. (6). See e.g. [38] for an algebraic proof.</footnote>
Before proceeding, we discuss some relevant properties of this relation and its constituents. Both coefficients are unique. The affine coefficient $J_e^{\mathcal{I}}$ represents the current along edge $e$ when every edge in $\mathcal{I}$ is removed from the graph. The linear coefficient $\lambda_{e \leftarrow i}^{\mathcal{I}}$ of $J_i(\mathbf{r}_{\mathcal{I}})$ represents a current-current susceptibility when considering the $n$ input currents (associated to the edges in $\mathcal{I}$). We show that $\lambda_{e \leftarrow i}^{\mathcal{I}}$ is equal to the one-input-current susceptibility in the graph obtained after removing edges $\mathcal{I} \setminus \{i\}$ (whose expression in terms of spanning tree polynomials was obtained in [33]). Furthermore, when modifying rates starting from a detailed-balance reference state, linear coefficients can be sampled from equilibrium fluctuations, as discussed in Sec. 5.2. The mutual multilinearity (MML) of Eq. (7) reduces to Eq. (2) in the case of one input current and to KCL when $n = n_c$. In this latter case, the current-current susceptibility $\lambda_{e \leftarrow i}^{\mathcal{I}}$, which can be obtained explicitly using graph-theoretical analysis, reduces to an algebraic coefficient consistent with KCL, as shown in Sec. 5.1. We stress that Eq. (7) is a property of stationary states, likewise is KCL; a generalization of Eq. (7) to finite-time in the Laplace domain is discussed in Sec. 4.

## 2.1 Illustration

Equipped with the main result, we briefly discuss an illustration of the mutual multilinearity relation before its proof is presented in the next sections. We start with two input currents as a natural extension of the previous result [33], keeping in mind that arbitrary sized sets of input currents can be considered if they are admissible.

Consider the graph in Fig. 4 with two input currents $J_1$ and $J_2$, respectively supported by edges $i_1$ and $i_2$, and output current $J_e$ supported by $e$. This graph, which has 12 states and 5 cycles, remains connected upon removal of $\mathcal{I} = \{i_1, i_2\}$, satisfying the admissibility criterion for the set of input currents. Therefore, MML establishes that

$$J_e(\mathbf{r}_{\mathcal{I}}) = J_e^{\mathcal{I}} + \lambda_{e \leftarrow i_1}^{\mathcal{I}} J_{i_1}(\mathbf{r}_{\mathcal{I}}) + \lambda_{e \leftarrow i_2}^{\mathcal{I}} J_{i_2}(\mathbf{r}_{\mathcal{I}}). \quad (8)$$

Notice that this is the equation of a plane in the space $\{J_1, J_2, J_e\}$, identified by the coefficients $J_e^{\mathcal{I}}$, $\lambda_{e \leftarrow i_1}^{\mathcal{I}}$ and $\lambda_{e \leftarrow i_2}^{\mathcal{I}}$, which only depend on the graph topology and rates $r_{\mathcal{E} \setminus \mathcal{I}}$. Therefore, any choice of input rates $\mathbf{r}_{\mathcal{I}}$ will yield a point on the plane.

In Fig. 4, we plot this plane and observe that all points indeed fall on the plane for random choices of input rates, illustrating the validity of the equation.

<figure><img src="image_7.png" /><figcaption>Figure 4: For the graph on the right, and its highlighted choices of input and output currents, the plot shows the plane obtained from Eq. (7) by calculating the affine and linear coefficients. Also, red spheres represent the values $\{J_1, J_2, J_e\}$ calculated by solving the master equation for 200 values of input rates randomly sampled from the uniform distribution $U([0, 20]^4)$.</figcaption></figure>
### 3 Graph-theoretic steady-state proof

We proceed by induction by proving the result for $n = n_c - 1$ and then by decreasing recursively the cardinality $n$ of the admissible set of edges by steps of 1.

Notice that an admissible set of input currents is any set such that removal of the corresponding edges does not disconnect the graph (otherwise their complement could not contain a spanning tree, which by definition connects all vertices); vice-versa, the complement of any set of edges whose removal does not disconnect the graph contains a spanning tree, whose complement is a chord set. See Figure 3.

#### 3.1 Proof of inductive seed

Consider an admissible set $\mathcal{I}$ of $n_c - 1$ edges. By admissibility, one can find an edge $\alpha \in \mathcal{E} \setminus \mathcal{I}$ such that $\mathcal{I}_\alpha := \mathcal{I} \cup \{\alpha\}$ is a fundamental set of chords, i.e. such that $\mathcal{T} = \mathcal{E} \setminus \mathcal{I}_\alpha$ is a spanning tree. From KCL, we know that any edge currents $j_e$ with $e \in \mathcal{T}$ can be expressed as a linear combination (with coefficients in $\{0, \pm 1\}$) of the currents $(j_i)_{i \in \mathcal{I}_\alpha}$. Our goal is to show that currents $(j_e)_{e \in \mathcal{E} \setminus \mathcal{I}}$ are linear-affine functions of the input currents $(j_i)_{i \in \mathcal{I}}$, with coefficients that depend only on the rates $(r_e)_{e \in \mathcal{E} \setminus \mathcal{I}}$. From KCL, to show this, it is sufficient to prove that $j_\alpha$ can itself be represented as a linear-affine combination of the currents $(j_i)_{i \in \mathcal{I}}$ with coefficients that depend only on the rates $(r_e)_{e \in \mathcal{E} \setminus \mathcal{I}}$.

Our strategy is to produce two linear relations between the stationary probability on one particular vertex $x_0$ and the stationary currents, taking care that the coefficients depend appropriately on the rates. We can then use such relations to eliminate the stationary probability and one of the currents, and deploy Kirchhoff's current law.

A useful observation to simplify the analysis below is that, for a fixed spanning tree $\mathcal{T}$ and a chord $\alpha$ (that we pick as defined above), there always exists a parametrization of the rates in terms of:

* a so-called potential $u = (u_x)_{x \in \mathcal{X}}$ (defined up to an arbitrary ground potential), that depends only on the rates $(r_e)_{e \in \mathcal{T}}$,

* a so-called cycle affinity $a_\alpha$ (associated to the oriented cycle $\mathcal{C}_\alpha$ generated when adding $\alpha$ to $\mathcal{T}$), that depends only on the rates $(r_e)_{e \in \mathcal{T}}$ and on the two rates $r_\alpha$,

such that (see App. A):

$$ \frac{r_{x \leftarrow y}}{r_{y \leftarrow x}} = e^{u_y - u_x}, \quad \text{for } x \leftarrow y \in \bar{\mathcal{T}}, \qquad (9) $$

$$ \frac{r_{x \leftarrow y}}{r_{y \leftarrow x}} = e^{u_y - u_x \pm a_\alpha}, \quad \text{for } x \leftarrow y = \pm \alpha. \qquad (10) $$

With this in mind, we turn to the current-probability relation that we will use below and that follows straightforwardly from the definition of stationary current Eq. (4). For $x \leftarrow y \in \bar{\mathcal{E}}$:

$$ \pi_x = \frac{r_{x \leftarrow y} \pi_y + j_{y \leftarrow x}}{r_{y \leftarrow x}}. \qquad (11) $$

This suggests to define the following pseudo-potential $\mu$ and pseudo-affinity $h$

$$ \mu_x = e^{u_x} \pi_x, \quad \text{for } x \in \mathcal{X}, \quad \text{and} \quad h_{y \leftarrow x} = e^{u_x} \frac{J_{y \leftarrow x}}{r_{y \leftarrow x}} = e^{u_y} \frac{J_{y \leftarrow x}}{r_{x \leftarrow y}}, \quad \text{for } x \leftarrow y \in \bar{\mathcal{T}}, \quad (12) $$

such that, from Eqs. (9) and (11),

$$ \mu_x - \mu_y = h_{y \leftarrow x}, \quad \forall x \leftarrow y \in \bar{\mathcal{T}}. \qquad (13) $$
Note that $h_{x \leftarrow y} = -h_{y \leftarrow x}$ as easily checked, provided $x \leftarrow y$ belongs to $\mathcal{T}$. The addition of the chord $\alpha$ to the spanning tree $\mathcal{T}$ produces an oriented cycle $\mathcal{C}_\alpha$ with $k+1$ vertices ($x_0, x_1, \dots, x_k$) and $k+1$ edges ($x_0 \leftarrow x_1, x_1 \leftarrow x_2, \dots, x_{k-1} \leftarrow x_k, \alpha = x_k \leftarrow x_0$). On the chord $\alpha$, the application of Eqs. (10) and (11) yields

$$ \mu_{x_0} - e^{-a_\alpha} \mu_{x_k} = h_\alpha, \quad \text{with} \quad h_\alpha = e^{u_{x_0}} \frac{j_\alpha}{r_\alpha}. \qquad (14) $$

The field $h$ is a pseudo-affinity in the sense that, for a chord like $\alpha$, we have an extra prefactor $e^{-a_\alpha}$ in Eq. (14) compared to Eq. (13), and we also have $h_{-\alpha} = -e^{a_\alpha}h_\alpha$ as easily checked.

The first linear relation is found as follows. Starting from edge $x_0 \leftarrow x_1$, we propagate Eq. (13) around the cycle $\mathcal{C}_\alpha$ using a telescopic sum:

$$ \mu_{x_k} - \mu_{x_0} = \sum_{0 \le p < k} (\mu_{x_{p+1}} - \mu_{x_p}) = \sum_{0 \le p < k} h_{x_p \leftarrow x_{p+1}}. \qquad (15) $$

On the last edge $\alpha = x_k \leftarrow x_0$ of the cycle we use, from Eq. (14), $\mu_{x_k} = e^{a_\alpha} \mu_{x_0} - e^{a_\alpha} h_\alpha$ to get

$$ (e^{a_\alpha} - 1) \mu_{x_0} = e^{a_\alpha} h_\alpha + \sum_{0 \le p < k} h_{x_p \leftarrow x_{p+1}}. \qquad (16) $$

Notice that if the affinity vanishes ($a_\alpha = 0$), the above expression already yields a linear relation among the currents. We investigate this further below, but for now we assume that $a_\alpha \neq 0$.

The second linear relation comes from the normalization $1 = \sum_{y \in \mathcal{X}} \pi_y = \sum_{y \in \mathcal{X}} e^{-u_y} \mu_y$. We root the spanning tree $\mathcal{T}$ in $x_0$ by orienting all its edges towards $x_0$. For any vertex $y \neq x_0$, Eq. (13) implies that $\mu_y$ can be found by integrating the pseudo-affinity $h$ along the unique path along $\mathcal{T}$ between $y$ and $x_0$. Defining the vector $h = (h_e)_{e \in \mathcal{T}}$ (with $|X| - 1$ entries), this integration procedure is encoded by a $(|X| - 1) \times (|X| - 1)$ path matrix $\mathbf{P}$, with entries $\in \{0, 1\}$, such that:

$$ \mu_y = (\mathbf{Ph})_y + \mu_{x_0}, \quad \text{for all } y \neq x_0. \qquad (17) $$

The matrix $\mathbf{P}$ can be found by inverting a “core” submatrix of the incidence matrix of the graph (see for instance Ref. [38]). Inserting Eq. (17) in the normalization condition, we arrive at

$$ \mu_{x_0} = \frac{1 - \sum_{\substack{y \in \mathcal{X} \\ y \neq x_0}} e^{-u_y} (\mathbf{Ph})_y}{\sum_{y \in \mathcal{X}} e^{-u_y}}. \qquad (18) $$

From KCL, the currents $(j_e)_{e \in \mathcal{T}}$ are a linear combination of the chord currents $(j_i)_{i \in \mathcal{I}}$, with coefficients in $\{0, \pm 1\}$. From Eq. (12), this implies that there exists a $(|\mathcal{X}| - 1) \times (n_c - 1)$ matrix $\mathbf{K}$ and a vector $\mathbf{k}^\alpha$ with $|\mathcal{X}| - 1$ entries, both depending on the rates $r_{\mathcal{T}}$, such that

$$ \mathbf{h} = \mathbf{K} \mathbf{J}_{\mathcal{I}} + \mathbf{J}_{\alpha} \mathbf{k}^{\alpha}, \qquad (19) $$

where $\mathbf{J}_{\mathcal{I}} = (j_i)_{i \in \mathcal{I}}$ is the vector of input currents. From this and Eq. (14), we insert Eq. (18) into Eq. (16) so as to eliminate $\mu_{x_0}$. We arrive at a linear-affine relation between $\mathbf{J}_{\alpha}$ and the input currents $\mathbf{J}_{\mathcal{I}}$:

$$ \frac{e^{a_\alpha} - 1}{\sum_{y \in \mathcal{X}} e^{-u_y}} \left[ 1 - \sum_{\substack{y \in \mathcal{X} \\ y \neq x_0}} e^{-u_y} (\mathbf{PK}\mathbf{J}_{\mathcal{I}})_y \right] - \mathbf{c} \cdot \mathbf{K}\mathbf{J}_{\mathcal{I}} = \left[ \frac{e^{a_\alpha + u_{x_0}}}{r_\alpha} + \mathbf{c} \cdot \mathbf{k}^\alpha + \frac{e^{a_\alpha} - 1}{\sum_{y \in \mathcal{X}} e^{-u_y}} \sum_{\substack{y \in \mathcal{X} \\ y \neq x_0}} e^{-u_y} (\mathbf{PK}^\alpha)_y \right] \mathbf{J}_\alpha, \qquad (20) $$
where c is the vector of entries in {0, 1} such that $\sum_{0 \le p < k} h_{x_p} \leftarrow x_{p+1} = c \cdot h$. The prefactor of $J_\alpha$ in the r.h.s. of Eq. (20), which is a function of the rates $r_{\mathcal{E} \setminus \mathcal{I}}$, cannot cancel. Indeed, if it were so, the l.h.s. of Eq. (20) would be zero $\forall r_\mathcal{I}$, implying a linear-affine interdependence of the input currents $J_\mathcal{I}$, with coefficients depending on $r_{\mathcal{E} \setminus \mathcal{I}} = r_{\mathcal{E} \cup \{\alpha\}}$ — but, from the independence lemma proved in App. B, this is impossible, since then the affine contribution $(e^{a_\alpha} - 1)/\sum_{y \in \mathcal{X}} e^{-u_y}$ would be 0; which is absurd.

To conclude, we infer from Eq. (20) that the chord current $J_\alpha$ is a linear-affine combination of the input currents $J_\mathcal{I}$ with coefficients that are well-defined functions of the rates $r_{\mathcal{E} \setminus \mathcal{I}}$. Using KCL, this shows that any current $(J_e)_{e \in \mathcal{E} \setminus \mathcal{I}}$ is a linear-affine combination of the input currents $J_\mathcal{I}$ with coefficients that are functions of the rates $r_{\mathcal{E} \setminus \mathcal{I}}$, which ends the proof of the initialization step of the recursion.

**Remark:** The “cycle equilibrium” case when $a_\alpha = 0$ is even simpler, and does not involve the normalization condition. Eq. (16) indeed then reads $h_\alpha = -c \cdot h$ which, using Eq. (19), yields

$$ \left( \frac{e^{u_{x_0}}}{r_{\alpha}} + c \cdot k^{\alpha} \right) J_{\alpha} = -c \cdot K J_{\mathcal{I}}, \qquad (21) $$

which, following the same line of arguments as above, allows one to conclude.

## 3.2 Proof of inductive step

By assuming that mutual linearity holds for a set of $n < n_c$ input edges, we want to prove that it also holds for $n - 1$.

Consider $\mathcal{I}$ a set of $n-1$ input edges under the single assumption that its removal does not disconnect the network, i.e. it is admissible. Our aim is to show that the current $J_e$ of any edge $e \in \mathcal{E} \setminus \mathcal{I}$ can be decomposed in a linear-affine way on the $J_{i \in \mathcal{I}}$ with coefficients independent on the rates $r_\mathcal{I}$. To show this, we distinguish two cases.

Let us first assume that $e$ is such that $\mathcal{I} \cup \{e\}$ is also admissible, and define $e_2 := e$. Since $|\mathcal{I}| \le n_c - 2$, it is always possible to find another edge $e_1 \notin \mathcal{I} \cup \{e_2\}$ such that both $\mathcal{I}_1 := \mathcal{I} \cup \{e_1\}$ and $\mathcal{I}_2 := \mathcal{I} \cup \{e_2\}$ are admissible. Since we assumed mutual linearity holds at the level of $n$ input edges, the (output) current supported by an edge $e_2 \notin \mathcal{I}_1$ can be expressed as

$$ J_{e_2} = J_{e_2}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) + \sum_{i_1 \in \mathcal{I}_1} \lambda_{e_2 \leftarrow i_1}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) J_{i_1}, \qquad (22) $$

where the right-hand side is linear-affine on the input currents with coefficients that depend on the rates of edges in $\mathcal{E} \setminus \mathcal{I}_1$. Similarly, for the edge $e_1$, because $\mathcal{I}_2$ is admissible, we can write

$$ J_{e_1} = J_{e_1}^{\mathcal{I}_2}(r_{\mathcal{E} \setminus \mathcal{I}_2}) + \sum_{i_2 \in \mathcal{I}_2} \lambda_{e_1 \leftarrow i_2}^{\mathcal{I}_2}(r_{\mathcal{E} \setminus \mathcal{I}_2}) J_{i_2}. \qquad (23) $$

We can now replace Eq. (23) on the right-hand side of (22) and collect $J_{e_2}$, yielding

$$ \begin{aligned} [1 - \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) \lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}(r_{\mathcal{E} \setminus \mathcal{I}_2})] J_{e_2} = & J_{e_2}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) J_{e_1}^{\mathcal{I}_2}(r_{\mathcal{E} \setminus \mathcal{I}_2}) \\ & + \sum_{i \in \mathcal{I}} [\lambda_{e_2 \leftarrow i}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2}(r_{\mathcal{E} \setminus \mathcal{I}_2})] J_i. \end{aligned} \qquad (24) $$

The multiplicative factor $1 - \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(r_{\mathcal{E} \setminus \mathcal{I}_1}) \lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}(r_{\mathcal{E} \setminus \mathcal{I}_2})$ cannot vanish (for any value of the rates), as shown in App. C, which proves the multilinearity of $J_{e_2}$ for the input set $\mathcal{I}$. Indeed from Eq. (24), $J_{e_2}$ is expressed as a linear-affine combination of the currents $J_{i \in \mathcal{I}}$ with coefficients depending only on the rates of the transitions in $(\mathcal{E} \setminus \mathcal{I}_1) \cup (\mathcal{E} \setminus \mathcal{I}_2) = \mathcal{E} \setminus \mathcal{I}$.
The argument works for any $e_2$ such that $\mathcal{I} \cup \{e_2\}$ is admissible. Consider now instead an edge $e$ such that removing $\mathcal{I} \cup \{e\}$ disconnects the network. We can fix a full set of chords that includes $\mathcal{I}$ (which is always possible since $\mathcal{I}$ is admissible). Since the removal of $e$ and $\mathcal{I}$ disconnect the network, by definition, $e$ belongs to a cocycle of the network [38]. If it is the only component of the cocycle, it means that $e$ is a bridge and the current $j_e$ is always zero. If instead the cocycle contains other edges, these are chords forming a subset $\mathcal{I}' \subseteq \mathcal{I}$ (because removing $\mathcal{I} \cup \{e\}$ disconnects the network). Then as a property of cocycles, the current $j_e$ can be written as a linear combination of the $j_{i' \in \mathcal{I}'}$, with coefficients $\pm 1$.

To summarize, in the second case where $\mathcal{I} \cup \{e\}$ is not admissible, the standard KCL directly gives the result (with rate-independent linear-affine coefficients), while in the first case ($\mathcal{I} \cup \{e\}$ admissible) the coefficients do depend on the rates and the derivation did require the recursion hypothesis at step $n$, completing the inductive step.

We conclude with a remark. Eq. (24) implies a recursive relation among susceptibilities associated with different input sets. Consider, for instance, the case where the input set consists of a single edge, $\mathcal{I} = \{e_3\}$. In this case we read from Eq. (24):

$$ \lambda_{e_2 \leftarrow e_3}^{\mathcal{I}}(\mathbf{r}_{\mathcal{E} \leftarrow e_3}) = \frac{\lambda_{e_2 \leftarrow e_3}^{\mathcal{I}_1}(\mathbf{r}_{\mathcal{E} \leftarrow e_1, e_3}) + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(\mathbf{r}_{\mathcal{E} \leftarrow e_1, e_3})\lambda_{e_1 \leftarrow e_3}^{\mathcal{I}_2}(\mathbf{r}_{\mathcal{E} \leftarrow e_2, e_3})}{1 - \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(\mathbf{r}_{\mathcal{E} \leftarrow e_1, e_3})\lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}(\mathbf{r}_{\mathcal{E} \leftarrow e_2, e_3})} \quad (25) $$

On the left hand side, we have the susceptibility onto edge $e_2$ of edge $e_3$ when taken as single input edge. Notice that this generally depends on $\mathbf{r}_1$ and we can trace this dependence through the terms $\lambda_{e_1 \leftarrow e_3}^{\mathcal{I}_2}(\mathbf{r}_{\mathcal{E} \leftarrow e_2, e_3})$ and $\lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}(\mathbf{r}_{\mathcal{E} \leftarrow e_2, e_3})$ on the right hand side. Consider now pruning the graph from edge $e_1$ (or equivalently setting rates $\mathbf{r}_1$ identically to zero). In this case, no current flows through edge $e_1$ and thus $\lambda_{e_1 \leftarrow e_3}^{\mathcal{I}_2}(\mathbf{r}_{\mathcal{E} \leftarrow e_2, e_3})$, $\lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}(\mathbf{r}_{\mathcal{E} \leftarrow e_2, e_3}) \to 0$. As a result, Eq. (25) reduces to:

$$ \lambda_{e_2 \leftarrow e_3}^{\mathcal{I}}(\mathbf{r}_{\mathcal{E} \leftarrow e_3}) \Big|_{r_1=0} = \lambda_{e_2 \leftarrow e_3}^{\mathcal{I}_1}(\mathbf{r}_{\mathcal{E} \leftarrow e_1, e_3}), \quad (26) $$

where the right hand side is unaffected by the pruning as it does not depend on $\mathbf{r}_1$. Interestingly, Eq. (26) tells that the susceptibility $\lambda_{e_2 \leftarrow e_3}$ with input set $\mathcal{I}_1 = \mathcal{I} \cup \{e_1\}$ on the full graph can be obtained from the corresponding susceptibility with smaller input set $\mathcal{I}$ on a smaller graph (pruned by edge $e_1$). As we will prove in Sec. 4.3, this recursive structure generalizes to arbitrary input set sizes. Thus, susceptibilities can be evaluated using the case of a single input current—where simpler analytical expressions have been established [33]—by considering the graph with all other input edges removed.

## 4 Laplace-space finite-time proof

MML in Eq. (7) is a property of steady state: already with one input current, the study of simple examples shows that it does not hold if we were to replace stationary currents by finite-time ones [33]. Similarly, with $n_c$ input currents MML reduces to KCL (see Sec. 5.1), which only holds for stationary currents.

In this section, we show how one can generalize Eq. (7) by going to the Laplace domain, introducing a frequency variable $\zeta$ conjugated to time $t$. By sending $\zeta$ to zero (equivalent to sending $t$ to infinity), this provides another proof of Eq. (7) with an explicit expression of the current-current susceptibility $\lambda_{e \leftarrow i}^{\mathcal{I}}$. Interestingly, the proof offers a physical representation of such a susceptibility: $\lambda_{e \leftarrow i}^{\mathcal{I}}$ is equal to the susceptibility $\tilde{\lambda}_{e \leftarrow i}^{i}$ for one input current (along $i$) and one output current (along $e$) in a modified graph where every edge of $\mathcal{I}$ (apart from $i$) is removed. Such a representation is also valid at non-zero $\zeta$. This allows one to represent $\lambda_{e \leftarrow i}^{\mathcal{I}}$ using tree ensembles, using results that we derived in Ref. [33].
## 4.1 Settings

The case of non-stationary currents can be understood by turning to the frequency domain, as follows. The probability distribution at time $t$ is the solution $\mathbf{p}(t) = \exp(t\mathbf{R})\mathbf{p}(0)$ to the master equation, given an initial probability distribution $\mathbf{p}(0)$. Defining the Laplace transform $\hat{\mathbf{p}}(\zeta) = \int_0^\infty dt e^{-\zeta t} \mathbf{p}(t)$ (and similarly for other functions of time), one arrives at the expression $\hat{\mathbf{p}}(\zeta) = (\zeta\mathbf{1} - \mathbf{R})^{-1}\mathbf{p}(0)$. Both the Laplace transform $\hat{\mathbf{p}}(\zeta)$ and the resolvent $(\zeta\mathbf{1} - \mathbf{R})^{-1}$ are defined for $\zeta \in \mathbb{C} \setminus \text{Sp} \mathbf{R}$, i.e. for complex numbers not belonging to the spectrum of $\mathbf{R}$.

We first prove a Laplace-domain generalization of MML:

$$ \hat{j}_e(\mathbf{r}_I, \zeta) = \hat{j}_e^I(\zeta) + \sum_{i \in \mathcal{I}} \hat{\lambda}_{e \leftarrow i}^I(\zeta) \hat{j}_i(\mathbf{r}_I, \zeta), \quad (27) $$

and then analyze the $\zeta \to 0$ limit in Sec. 4.3 where we recover Eq. (7) (holding at stationarity $t \to \infty$).

## 4.2 Proof of mutual multilinearity

Proving Eq. (27) is equivalent to show the existence of constants $\hat{\lambda}_{e \leftarrow i}^I(\zeta)$'s such that the gradient of both sides w.r.t. the $2n$ variables $\mathbf{r}_I$ is the same. In other words, denoting $\hat{\lambda}^I(\zeta)$ the vector of $n$ components $(\hat{\lambda}^I(\zeta))_i = \hat{\lambda}_{e \leftarrow i}^I(\zeta)$ for $i \in \mathcal{I}$, we have to show that

$$ \exists \hat{\lambda}^I(\zeta) : \forall k \in \mathcal{I}, \forall \epsilon \in \{-1, +1\}, \quad \partial_{r_{ek}} \hat{j}_e(\mathbf{r}_I, \zeta) = \sum_{i \in \mathcal{I}} \hat{\lambda}_{e \leftarrow i}^I(\zeta) \partial_{r_{ek}} \hat{j}_i(\mathbf{r}_I, \zeta). \quad (28) $$

The key technical ingredient is the expression, obtained in Ref. [33], of the partial derivative of currents with respect to rates $r_{ek}$ (with $\epsilon \in \{-1, +1\}$), as proportional to a matrix element of the resolvent of $\mathbf{V}^\top \mathbf{S}$:

$$ \partial_{r_{ek}} \hat{j}_i(\zeta) = \epsilon \zeta \langle \mathbf{s}(\epsilon k) | \hat{\mathbf{p}}(\zeta) \rangle \langle i | (\zeta \mathbf{1} + \mathbf{V}^\top \mathbf{S})^{-1} | k \rangle. \quad (29) $$

Here $\mathbf{s}(\pm k)$ (resp. $\mathbf{s}(-k)$) is the source (resp. target) of edge $k$, and $\mathbf{S}$ and $\mathbf{V}$ are respectively the incidence and the current matrices of dimensions $X \times E$ and entries<fnref n="2" />

$$ \mathbf{S}_{x,e} = \delta_{x,s(-e)} - \delta_{x,s(+e)}, \quad (30) $$

$$ \mathbf{V}_{x,e} = r_{-e} \delta_{x,s(-e)} - r_{+e} \delta_{x,s(+e)}, \quad (31) $$

where in this section we define $X = |\mathcal{X}|$ and $E = |\mathcal{E}|$ for compactness. We also use the bracket notation: $|k\rangle$ (resp. $\langle k|$) denotes the canonical column (resp. line) vector along the $k$-th direction.

The interest of Eq. (29) is that the prefactors of the matrix element of the resolvent often vanish by compensation. And indeed since Perron-Frobenius ensures $\langle \mathbf{s}(\epsilon k) | \hat{\mathbf{p}}(\zeta) \rangle > 0$, we see from Eq. (28) that we have to prove

$$ \exists \hat{\lambda}^I(\zeta) : \forall k \in \mathcal{I}, \quad \langle e | (\zeta \mathbf{1} + \mathbf{V}^\top \mathbf{S})^{-1} | k \rangle = \sum_{i \in \mathcal{I}} \hat{\lambda}_{e \leftarrow i}^I(\zeta) \langle i | (\zeta \mathbf{1} + \mathbf{V}^\top \mathbf{S})^{-1} | k \rangle \quad (32) $$

(notice the essential disappearance of the sign $\epsilon$), where $\hat{\lambda}^I(\zeta)$ has to be independent on the $2n$ rates $r_I$. Now, denoting $\mathbf{M} = \zeta\mathbf{1} + \mathbf{V}^\top \mathbf{S}$ and $(\mathbf{M}^{-1})_{I,J}$ the submatrix of $\mathbf{M}^{-1}$ where only lines and columns of tuples of distinct indices $I$ and $J$ are kept, we see that the condition of Eq. (32) transforms into a matrix equation for the vector $\hat{\lambda}^I(\zeta)$:

$$ \mathbf{M}^{-1} \mathbf{1}_{\mathcal{I},e} = \mathbf{M}^{-1} \mathbf{1}_{\mathcal{I},\mathcal{I}} \hat{\lambda}^I(\zeta). \quad (33) $$

<footnote n="2">²The current writes as $j(t) = -\mathbf{V}^\top \mathbf{p}(t)$ and the rate matrix as $\mathbf{R} = -\mathbf{S}\mathbf{V}^\top$. See Refs [33, 45, 47] for example uses.</footnote>
We now prove that this equation can be solved by inverting the $n \times n$ matrix $\mathbf{M}^{-1} {}_{\mathcal{I},\mathcal{I}}^{\top}$, and that the result does not depend on the rates $\mathbf{r}_{\mathcal{I}}$ of the input currents. A central identity to do so is Jacobi's formula, which states that for an invertible matrix $\mathbf{B}$, one has

$$ \det(\mathbf{B}^{-1})_{I,J} = \varepsilon_{IJ} \frac{\det \mathbf{B}_{(J,I)}}{\det \mathbf{B}}, \quad (34) $$

with $\mathbf{B}_{(J,I)}$ the submatrix of $\mathbf{B}$ where lines $J$ and columns $I$ are removed, and $\varepsilon_{IJ} \in \{-1, +1\}$ is a sign that depends on the tuples $I$ and $J$:

$$ \varepsilon_{IJ} = (-1)^{\sum_{i \in I} i + \sum_{j \in J} j} \sigma(I) \sigma(J). \quad (35) $$

Here $\sigma(I)$ is the signature of $I$, seen as the permutation of its elements sorted by increasing order (for instance $I = (1, 4, 2)$ is seen as the permutation of $(1, 2, 4)$ so that $\sigma(I) = -1$). We thus have<fnref n="3" />

$$ \det \mathbf{M}^{-1} {}_{\mathcal{I},\mathcal{I}}^{\top} = \frac{\det \mathbf{M}_{(\mathcal{I},\mathcal{I})}}{\det \mathbf{M}} = \frac{\det (\zeta \mathbf{1} + (\mathbf{V}_{(\cdot,\mathcal{I})})^{\top} \mathbf{S}_{(\cdot,\mathcal{I})})}{\det (\zeta \mathbf{1} + \mathbf{V}^{\top} \mathbf{S})} = \zeta^{-n} \frac{\det (\zeta \mathbf{1} - \mathbf{R}_{\setminus \mathcal{I}})}{\det (\zeta \mathbf{1} - \mathbf{R})}, \quad (36) $$

where $\mathbf{R}_{\setminus \mathcal{I}}$ is the rate matrix of a reduced network where the $n$ input edges $\mathcal{I}$ are removed. The matrix $\mathbf{M}^{-1} {}_{\mathcal{I},\mathcal{I}}^{\top}$ is thus invertible for complex $\zeta$ not belonging to the spectra of $\mathbf{R}$ or $\mathbf{R}_{\setminus \mathcal{I}}$ and thus, since these matrices are both stochastic, in a finite vicinity of $0$ in $\mathbb{C} \setminus \{0\}$ (which will allow us to study the $\zeta \to 0$ limit safely).

Eq. (33) is a linear equation of the form $\mathbf{y} = \mathbf{B}\mathbf{x}$ with $\mathbf{x}$ and $\mathbf{y}$ two vectors of length $n$ and $\mathbf{B}$ a $n \times n$ invertible matrix. Cramer's rule gives its solution in a determinantal form: the component $i$ of $\mathbf{x}$ writes as $x_i = \det \mathbf{B}_i / \det \mathbf{B}$ where $\mathbf{B}_i$ is obtained replacing column $i$ of $\mathbf{B}$ by the vector $\mathbf{y}$. In the case of Eq. (33), the vector on the l.h.s. is equal to the column labelled by $e$ of the matrix on the r.h.s. Cramer's rule thus offers an explicit expression of the components of the vector $\hat{\lambda}_{\mathcal{I}}^{\mathcal{I}}(\zeta)$ as

$$ \hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}(\zeta) = \frac{\det[\mathbf{M}^{-1}_{\mathcal{I}|_{i \to e}, \mathcal{I}}]}{\det[\mathbf{M}^{-1}_{\mathcal{I}, \mathcal{I}}]} \quad (37) $$

$$ = (-1)^{i+e} \sigma(\mathcal{I}) \sigma(\mathcal{I}|_{i \to e}) \frac{\det(\zeta \mathbf{1} + \mathbf{V}^{\top} \mathbf{S})_{(\mathcal{I}, \mathcal{I}|_{i \to e})}}{\det(\zeta \mathbf{1} + \mathbf{V}^{\top} \mathbf{S})_{(\mathcal{I}, \mathcal{I})}}, \quad (38) $$

where $\mathcal{I}|_{i \to e}$ denotes the tuple obtained from $\mathcal{I}$ by replacing label $i$ by $e$. For the second line we make use of Jacobi's identity, Eq. (34). From Eq. (38), we see explicitly that the susceptibility $\hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}(\zeta)$ is independent of the rates of the input edges $\mathcal{I}$, since $\mathbf{V}^{\top} \mathbf{S}$ depends on the two rates $\mathbf{r}_i$ only through its line $i$. (This is not obvious to read from Eq. (37) since the matrices involved on its r.h.s. depend explicitly on the rates $\mathbf{r}_i$.) This ends the proof of the Laplace-domain multilinearity, Eq. (27).

Note that the affine constant $\hat{j}_{e}^{\mathcal{I}}(\zeta)$ in Eq. (27) is simply given by the (Laplace) current for a reduced graph with the edges in $\mathcal{I}$ removed, since $\hat{j}_i(\mathbf{r}_{\mathcal{I}}, \zeta)$ is 0 when the rates $\mathbf{r}_{\mathcal{I}}$ are sent to 0.

<footnote n="3">³We use Sylvester's determinant theorem to write: $\det(\zeta \mathbf{1} + (\mathbf{V}_{(\cdot,\mathcal{I})})^{\top} \mathbf{S}_{(\cdot,\mathcal{I})}) = \zeta^{E-n-X} \det(\zeta \mathbf{1} - \mathbf{R}_{\setminus \mathcal{I}})$, and $\det \mathbf{M} = \zeta^{E-X} \det(\zeta \mathbf{1} - \mathbf{R})$.</footnote>
## 4.3 Taking the limit $\zeta \to 0$ and mapping to the one-input-current case

Using Sylvester's determinant theorem, the denominator of the expression of Eq. (38) of the Laplace-domain susceptibility $\hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}(\zeta)$ rewrites as

$$ \det(\zeta \mathbf{1} + \mathbf{V}^{\top} \mathbf{S})_{\setminus (\mathcal{I}, \mathcal{I})} = \zeta^{E-n-X} \det(\zeta \mathbf{1} - \mathbf{R}_{\setminus \mathcal{I}}) \quad (39) \\ \zeta \xrightarrow{\sim} \zeta^{E-n-X+1} \operatorname{tr} \operatorname{adj}(-\mathbf{R}_{\setminus \mathcal{I}}), $$

where the trace of the adjugate of $-\mathbf{R}_{\setminus \mathcal{I}}$ can be expressed as a sum of weights of rooted spanning trees of $\mathcal{G}_{\setminus \mathcal{I}}$, the Markov chain graph where edges from $\mathcal{I}$ are removed (this is an example of relation between determinants and spanning tree expansions [48], and see e.g. Refs. [33, 45, 49] for recent applications in the field of continuous-time Markov chains).

Understanding the $\zeta \to 0$ behavior of the numerator of Eq. (38) proves to be a little trickier. Consider the graph $\tilde{\mathcal{G}}$ obtained from $\mathcal{G}$ by removing the edges in $\mathcal{I} \setminus \{i\}$, where we fix for a moment an input edge $i \in \mathcal{I}$. This defines an auxiliary one-input-current problem, with input edge $i$ and output edge $e$ in the reduced graph $\tilde{\mathcal{G}}$. The graph $\tilde{\mathcal{G}}$ has $X$ vertices and $\tilde{E} = E - n + 1$ edges, and is connected (since removing $\mathcal{I}$ does not disconnect $\mathcal{G}$). As a special case of Eq. (38), the corresponding Laplace-domain current-current susceptibility, that we denote $\tilde{\lambda}_{e \leftarrow i}^{i}(\zeta)$, is expressed as

$$ \tilde{\lambda}_{e \leftarrow i}^{i}(\zeta) = (-1)^{\tilde{\iota}+\tilde{\bar{e}}} \frac{\det(\zeta \mathbf{1}_{\tilde{E}} + \tilde{\mathbf{V}}^{\top} \tilde{\mathbf{S}})_{\setminus (\tilde{\iota}, \tilde{\bar{e}})}}{\det(\zeta \mathbf{1}_{\tilde{E}} + \tilde{\mathbf{V}}^{\top} \tilde{\mathbf{S}})_{\setminus (\tilde{\iota}, \tilde{\bar{i}})}}, \quad (40) $$

where $\tilde{\mathbf{S}}$ is the incidence matrix of $\tilde{\mathcal{G}}$ and $\tilde{\mathbf{V}}$ its corresponding weighted version, defined similarly to Eqs. (30)-(31) (also, to facilitate the reading until the end of this section, we make explicit the dimension of identity matrices as indices). The integers $\tilde{\iota}$ and $\tilde{e}$ are the indices of input and output edges in the reduced set of edges of $\tilde{\mathcal{G}}$, so that$^4(-1)^{\tilde{\iota}+\tilde{e}} = (-1)^{i+e} \sigma(\mathcal{I}) \sigma(\mathcal{I}|_{i \leftarrow e})$. Using Sylvester's determinant theorem as previously, we have:

$$ \det(\zeta \mathbf{1}_{\tilde{E}} + \tilde{\mathbf{V}}^{\top} \tilde{\mathbf{S}})_{\setminus (\tilde{\iota}, \tilde{\bar{\iota}})} = \zeta^{\tilde{E}-1-X} \det(\zeta \mathbf{1}_{X} - \tilde{\mathbf{R}}_{\setminus \tilde{\iota}}) \stackrel{(39)}{=} \det(\zeta \mathbf{1}_{E} + \mathbf{V}^{\top} \mathbf{S})_{\setminus (\mathcal{I}, \mathcal{I})}, \quad (41) $$

where we used $\tilde{\mathbf{R}}_{\setminus \tilde{\iota}} = \mathbf{R}_{\setminus \mathcal{I}}$. Also, by direct inspection of definitions, we have:

$$ (\zeta \mathbf{1}_{\tilde{E}} + \tilde{\mathbf{V}}^{\top} \tilde{\mathbf{S}})_{\setminus (\tilde{\iota}, \tilde{\bar{e}})} = (\zeta \mathbf{1}_{E} + \mathbf{V}^{\top} \mathbf{S})_{\setminus (\mathcal{I}, \mathcal{I}|_{i \leftarrow e})}. \quad (42) $$

Collecting these results, we finally obtain

$$ \hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}(\zeta) = \tilde{\lambda}_{e \leftarrow i}^{i}(\zeta). \quad (43) $$

In other words, the susceptibility of the output current $e$ in the input current $i \in \mathcal{I}$ in the graph $\mathcal{G}$ is equal to the (one-input-current) susceptibility for the same input and output edges but in the reduced graph $\tilde{\mathcal{G}}$ where input edges of $\mathcal{I}$ other than $i$ are removed. Eq. (43) thus provides a interesting physical representation of the susceptibility in the Laplace-space mutual linearity of Eq. (27): the susceptibility $\hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}(\zeta)$ of the full problem can be computed from that of a simpler problem (one input and one output current, in a reduced graph).

Coming back to the $\zeta \to 0$ asymptotics, Eq. (43) implies that $\hat{\lambda}_{e \leftarrow i}^{\mathcal{I}} = \lim_{\zeta \to 0} \hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}(\zeta)$ exists and can be expressed, following Ref. [33], as a sum of weights of trees constructed from the reduced graph $\tilde{\mathcal{G}}$. This ends the Laplace-space proof of the stationary MML [Eq. (7)], and since Eq. (43) extends to the $\zeta \to 0$ limit, also provides a representation of the stationary susceptibility $\hat{\lambda}_{e \leftarrow i}^{\mathcal{I}}$ as the one, $\tilde{\lambda}_{e \leftarrow i}^{i}$, of a simpler problem.

<footnote>$^4$From the definition of $\tilde{\iota}$ and $\tilde{e}$, we see that $(-1)^{\tilde{\iota}-\tilde{e}} = (-1)^{i-\tilde{e}}(-1)^{\tilde{\iota}}$ where $\tilde{\iota}$ is the cardinality of $\mathcal{J} = \mathcal{I} \cap \{i\}$, the set of elements of $\mathcal{I}$ strictly comprised between the two integers $i$ and $e$. Also, the composition of $\mathcal{I}|_{i \leftarrow e}$ and of the inverse of $\mathcal{I}$ (here, as discussed after Eq. (35), a sequence of integers is seen as the permutation of its elements sorted by increasing order) is a cycle permutation of length $\tilde{\iota} + 1$, because it shifts $\{e\} \cup \mathcal{J}$ to $\mathcal{J} \cup \{i\}$. Its signature $\sigma(\mathcal{I}^{-1} \circ \mathcal{I}|_{i \leftarrow e}) = \sigma(\mathcal{I})\sigma(\mathcal{I}|_{i \leftarrow e})$ is thus equal to $(-1)^{\tilde{\iota}}$, which proves the announced result.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 5: Resulting graph after removal of $\mathcal{I}$ and $e$. Left: $e$ and $i$ do not belong to the same cycle $\mathcal{C}$, notice that $i$ is fully contained by $\mathcal{T}^-$. Right: Since $e$ and $i$ belong to the same cycle, $i$ is in between both trees $\mathcal{T}^\pm$.</figcaption></figure>

# 5 Further developments

## 5.1 Recovering Kirchhoff's current law

MML in Eq. (7) generalizes KCL in the sense that it extends the validity of Eq. (6) to subsets of currents that do not constitute a fundamental set, cf. Fig. 2. Importantly, we now show that Eq. (7) yields back Eq. (6) when the input set $\mathcal{I}$ is fundamental as a result of the terms $J_e^{\mathcal{I}}$ vanishing and the terms $\lambda_{e \leftarrow i}^{\mathcal{I}}$ reducing to 0 or $\pm 1$.

Let $\mathcal{I}$ be a fundamental, and therefore also admissible, set. It is immediate that $J_e^{\mathcal{I}} = 0$ for all edges $e$ since there are no stationary currents on a tree.

Regarding the susceptibility $\lambda_{e \leftarrow i}^{\mathcal{I}}$ of a current over $e$ with respect to $i \in \mathcal{I}$, we have shown in Sec. 4.3 that it is equivalent to the susceptibility $e \leftarrow i$ when $i$ is a single input current in the reduced graph obtained from $\mathcal{G}$ by removing the edges $\mathcal{I} \setminus \{i\}$. Then, using the expression derived in Ref. [33] of such a single-input-current susceptibility in terms of spanning tree polynomials, we arrive at:

$$ \lambda_{e \leftarrow i}^{\mathcal{I}} = r_{s(-e) \leftarrow s(+e)} \frac{\tau_{\mathcal{I},e}^{\prime s(+e) \to s(+i)} - \tau_{\mathcal{I},e}^{\prime s(+e) \to s(-i)}}{\tau_{\mathcal{I}}} - r_{s(+e) \leftarrow s(-e)} \frac{\tau_{\mathcal{I},e}^{\prime s(-e) \to s(+i)} - \tau_{\mathcal{I},e}^{\prime s(-e) \to s(-i)}}{\tau_{\mathcal{I}}}, \quad (44) $$

where

$$ \tau_{\mathcal{I}} := \sum_x \sum_{\mathcal{T}_x \subseteq \mathcal{G}} \omega(\mathcal{T}_x), \quad (45) $$

is the sum over all vertices $x$ and over all rooted spanning trees $\mathcal{T}_x$ in $\mathcal{G} := \mathcal{G} \setminus \mathcal{I}$, of the product of the transition rates $\omega(\mathcal{T}_x)$ belonging to $\mathcal{T}_x$; also,

$$ \tau_{\mathcal{I},e}^{\prime x \to y} := \sum_x \sum_{\mathcal{T}_x \subseteq \mathcal{G}': \{x \to y\} \subseteq \mathcal{T}_x} \omega(\mathcal{T}_x), \quad (46) $$

is a similar sum of polynomials for the graph $\mathcal{G}' := (\mathcal{G} \cup \{x \to y\}) \setminus \{e\}$ where an auxiliary edge $x \to y$ of rate 1 is introduced while all edges in $\mathcal{I}$ and edge $e$ are removed, and only rooted spanning trees that contain this auxiliary edge are considered. We will use Eq. (44) to explicitly obtain the susceptibilities in this case.

Since $\mathcal{I}$ constitutes a fundamental set, the first step of removing $\mathcal{I} \setminus \{i\}$ leads to a single cycle $\mathcal{C}$ and possibly branching leaves. Secondly, removing $i$ results in a tree $\mathcal{T} = \mathcal{E} \setminus \mathcal{I}$. Lastly, two disconnected trees sprout from the removal of $e$, $\mathcal{T}^+$ and $\mathcal{T}^-$. The former contains $s(+e) \subseteq \mathcal{T}^+$ while $s(-e) \subseteq \mathcal{T}^-$, see Fig. 5.
We first consider the case $e \notin \mathcal{C}$ (left panel of Fig. 5), in which case $s(\pm i)$ belong to one of the two trees. If $i \subseteq \mathcal{T}^{-}$,

$$ \tau'_{\backslash \mathcal{I},e} s(+e) \to s(+i) = \tau'_{\backslash \mathcal{I},e} s(-i) = \omega(\mathcal{T}_{s(+e)}^{+}) \sum_{x \in \mathcal{T}^{-}} \omega(\mathcal{T}_{x}^{-}), \quad (47) $$

because the added edges $s(+e) \to s(\pm i)$ act as bridges between the tree $\mathcal{T}^{+}$ and $\mathcal{T}^{-}$, while

$$ \tau'_{\backslash e, \mathcal{I}} s(-e) \to s(+i) = \tau'_{\backslash e, \mathcal{I}} s(-i) = 0, \quad (48) $$

because $s(\pm i)$ and $s(-e)$ belong all to $\mathcal{T}^{-}$ and the resulting graph is disconnected (no edge connects $\mathcal{T}^{+}$ and $\mathcal{T}^{-}$) and therefore has no spanning trees. Consequently, $\lambda_{e \leftarrow i}^{\mathcal{I}} = 0$, and similarly for $i \subseteq \mathcal{T}^{+}$.

The second case, where $e \in \mathcal{C}$, results in the source and target of $+i$ belonging to distinct trees. If $s(+i) \subseteq \mathcal{T}^{-}$,

$$ \tau'_{\backslash \mathcal{I},e} s(+e) \to s(+i) = \omega(\mathcal{T}_{s(+e)}^{+}) \sum_{x \in \mathcal{T}^{-}} \omega(\mathcal{T}_{x}^{-}), \quad (49) $$

$$ \tau'_{\backslash \mathcal{I},e} s(-e) \to s(-i) = \omega(\mathcal{T}_{s(-e)}^{-}) \sum_{x \in \mathcal{T}^{+}} \omega(\mathcal{T}_{x}^{+}), \quad (50) $$

and

$$ \tau'_{\backslash e, \mathcal{I}} s(+e) \to s(-i) = \tau'_{\backslash e, \mathcal{I}} s(-i) = 0. \quad (50) $$

Plugging these values in Eq. (44), we identify the numerator as equivalent to $\tau_{\backslash \mathcal{I}}$ and therefore $\lambda_{e \leftarrow i}^{\mathcal{I}} = 1$. Finally, if $s(+i) \subseteq \mathcal{T}^{+}$, the same procedure leads to $\lambda_{e \leftarrow i}^{\mathcal{I}} = -1$.

In summary, for an input set that is also fundamental, any output current can be written as a linear combination of input currents with coefficient 0 when $e$ does not lie in the same cycle of $i$ in the graph $\mathcal{E} \setminus (\mathcal{I} \setminus \{i\})$, and $\pm 1$ otherwise, with sign depending on the convention directions, which is Eq. (6).

## 5.2 Connection to equilibrium fluctuations

Here we show that when currents are controlled starting from an equilibrium state, the susceptibility can be interpreted in terms of currents’ covariances at equilibrium. This brings with and extends known results in so-called linear regime theory, but in a fully nonequilibrium context.

First we need this:

$$ \begin{aligned} J_e &= \lim_{t \to \infty} \frac{\langle J_e(t) \rangle}{t}, \\ c_{ee'} &= \lim_{t \to \infty} \frac{\langle J_e(t) J_{e'}(t) \rangle_c}{t} = \lim_{t \to \infty} \frac{\langle J_e(t) J_{e'}(t) \rangle - \langle J_e(t) \rangle \langle J_{e'}(t) \rangle}{t}, \end{aligned} \quad (51) $$

where $J_e(t)$ are the cumulated currents up to time $t$ along the stochastic process and $\langle \cdot \rangle$ denote sample averages. When the sample mean is over an equilibrium process we have $J_e = 0$ and denote $c_{ee'} = c_{ee'}^{\text{eq}}$. Notice that by definition covariances are symmetric: $c_{ee'} = c_{e'e}$.

The main result of this section is that, starting from equilibrium (where $J_e^{\text{eq}} = 0$ for all $e$, thus the rates satisfy detailed balance) and perturbing the rates $r_i$ of an input current arbitrarily far from equilibrium, the following linear response relation among currents holds

$$ J_e = \frac{c_{ei}^{\text{eq}}}{c_{ii}^{\text{eq}}} J_i, \quad (52) $$

where $c_{ei}^{\text{eq}}$ are the currents’ steady covariances at equilibrium. Therefore, properties of response arbitrarily far from equilibrium can be inferred from equilibrium fluctuations.
The proof of the result is straightforward given the results of Refs. [33, 38] and linear response theory (see e.g. [50] for a review and [51] for closely related results). For a moment, let us parametrize the rates according to Eqs. (9)-(10), where we fix the potentials $u_{x \in \mathcal{X}}$ and take the cycle affinities $\mathbf{a} = (a_\alpha)$ (where $\alpha$ runs through chords) as control parameters to drive the dynamics out of the equilibrium reference state $\mathbf{a} = \mathbf{0}$. The following fluctuation-dissipation relation holds (coming from linear response theory, see App. D for a full derivation):

$$ \left. \frac{\partial J_e}{\partial a_\alpha} \right|_{\text{eq.}} = \frac{1}{2} c_{e\alpha}^{\text{eq}}. \qquad (53) $$

Consider first the case of one admissible input current $i$ (that can thus be taken as a chord). From the results of Ref. [33], there exist two constants $J_e^i$ and $\lambda_{e \leftarrow i}^i$ such that for all values of the rates $\mathbf{r}_i$, we have $J_e(\mathbf{r}_i) = J_e^i + \lambda_{e \leftarrow i}^i J_i(\mathbf{r}_i)$. If we parametrize the rates by fixing the potentials $u_{x \in \mathcal{X}}$ and affinities $a_{\alpha \neq i}$, while controlling the currents through $a_i$, this also means that we have (with the same constants)

$$ J_e(a_i) = J_e^i + \lambda_{e \leftarrow i}^i J_i(a_i). \qquad (54) $$

For rates satisfying detailed balance ($\mathbf{a} = \mathbf{0}$) we have $J_e^i|_{\text{eq}} = 0$. Taking the derivative of Eq. (54) with respect to $a_i$, using Eq. (53) and rearranging terms we obtain $\lambda_{e \leftarrow i}^{\text{eq}} = c_{ei}^{\text{eq}}/c_{ii}^{\text{eq}}$ and thus Eq. (52).

The generalization to an arbitrary number of admissible input rates being modified is straightforward. Multilinearity reads

$$ J_e(\mathbf{r}_I) = J_e^I + \sum_{i \in I} \lambda_{e \leftarrow i}^I J_i(\mathbf{r}_I). \qquad (55) $$

We parametrize the rates with a potential $u_{x \in \mathcal{X}}$ and cycle affinities $\mathbf{a}$, and we consider equilibrium reference rates ($\mathbf{a} = \mathbf{0}$). We define the vectors $\lambda_e^I = (\lambda_{e \leftarrow i}^{\text{eq}})_{i \in I}$, $c_e^I = (c_{e,i}^{\text{eq}})_{i \in I}$, $J_I = (J_i)_{i \in I}$, $\mathbf{a}_I = (a_i)_{i \in I}$, and the matrix $C_I = (c_{k,i}^{\text{eq}})_{k,i \in I}$. Controlling currents through $\mathbf{a}_I$, Eq. (55) implies<fnref n="5" />

$$ J_e(\mathbf{a}_I) = \sum_{i \in I} \lambda_{e \leftarrow i}^{I \text{ eq}} J_i(\mathbf{a}_I). \qquad (56) $$

Since $I$ is admissible, it can be taken as a subset of chords, that we fix and this allows one to use Eq. (53) in this expression, which yields $c_e^I = C_I \lambda_e^I$. Now, we remark that $C_I$ is a principal submatrix of the full (i.e. all-chord) response matrix between cycle affinities and chord currents, which is symmetric definite positive (see e.g. [38]). Hence $C_I$ has the same property and is invertible. This allows one to write $\lambda_e^I = C_I^{-1} c_e^I$ and we obtain finally

$$ J_e(\mathbf{r}_I) = c_e^I \cdot C_I^{-1} J_I(\mathbf{r}_I), \qquad (57) $$

valid when controlling input currents starting from an equilibrium reference point.

The physical relevance of our result is that it is operationally accessible: the currents' susceptibilities do not have to be obtained analytically, but can be sampled from a realization of the process.

<footnote n="5">⁵In this expression, explicitly, $\lambda_{e \leftarrow i}^{\text{eq}}$ is given by $\lambda_{e \leftarrow i}^I$ evaluated for cycle affinities $a_\alpha = 0, \forall \alpha \notin I$.</footnote>
# 6 Conclusion and perspectives

Control over multiple transition rates is a scenario that emerges across disciplines. For example, tuning the concentration of chemicals in biochemical networks often affects several transition rates simultaneously. The same occurs when adjusting the temperature of a reservoir in a heat engine, modifying energy barriers in activated processes such as protein folding or ion transport, or altering transcription factors in gene expression dynamics. Understanding how stationary currents reorganize in response to changes in several transition rates is therefore central.

In this work, we have uncovered a mutual linear-affine relation between stationary currents when controlling multiple transition rates. Borrowing approaches and terminology from graph theory, we have considered a subset of chord currents as input currents, whose rates can be arbitrarily and independently manipulated. In this context, the Kirchhoff current law is a tenet of the theory of continuous-time Markov chains, which states that a complete set of chord currents fully specifies the currents at stationarity. Our result generalizes Kirchhoff current law by drawing a (linear-affine) relation between any stationary current and an incomplete set of chords. As such, it also generalized the recent result put forward by us in Ref. [33] for the case of a single input current.

Two noteworthy insights emerging from the ideas herein merit further discussion. First, as seen in the illustration of Sec. 2.1, MML bears a geometrical interpretation: it states the existence of a hyperplane in the space of input and output currents, on which all possible values of {$J_{i \in I}, J_e$} lie for any choice of input rates; nevertheless, the allowed points may form a bounded domain within the hyperplane, constrained by factors such as mean first-passage times between source and target of input edges, and its exact shape remains to be explored. Second, we have seen that susceptibilities can be calculated from the single input current case in the graph pruned from all other input currents; indicating that the influence of an isolated input current on the output is preserved once more edges are added to the graph and the corresponding currents are treated as input.

Our result is remarkably general, being valid for arbitrary finite network geometries and arbitrarily large changes of rates. Yet, we envisage a few limitations which call for future work. Extending our results to the framework of open (driven) systems is an interesting direction. In the case of a single input current we have proved that the linearity extends to the case of open systems despite the fact that global conservation of probability breaks down [33]. This suggests that such a generalization is viable also for multiple input currents and could lead to interesting applications in the field of adaptive physical networks such as flow, elastic and resistor networks [52–54]. A second, more far-fetched generalization is to incorporate non-linear interacting networks. This is the case of biochemical and metabolic networks, where interactions involve more than two chemicals (for example due to autocatalytic motifs [55]) and therefore maps to hypergraphs instead of Markov chains. Despite the apparent profound difference, recent work has uncovered similarities between noninteracting and interacting networks [38], which encourage one to pursue this path.

We conclude with a final, practical remark. Despite the result presented in this work representing a generalization from single to multiple input currents, it does not include the case of macroscopic input currents that are a linear combination of multiple edge currents. For example, suppose the case where two edge currents $J_1$ and $J_2$ are not independently measurable, but $J = J_1 + J_2$ is measurable. Our result states that

$$J_e = \lambda_e^{\{1,2\}} + \lambda_{e,1}^{\{1,2\}} J_1 + \lambda_{e,2}^{\{1,2\}} J_2 = \lambda_e^{\{1,2\}} + (\lambda_{e,1}^{\{1,2\}} - \lambda_{e,2}^{\{1,2\}}) J_1 + \lambda_{e,2}^{\{1,2\}} J, \quad (58)$$

and we immediately observe that we cannot get rid of the unmeasurable current $J_1$ unless $\lambda_{e,1}^{\{1,2\}} = \lambda_{e,2}^{\{1,2\}}$. As pointed out in Ref. [33], however, the result easily extends to the case of macroscopic output currents that are a linear combination of multiple edge currents.
Acknowledgments

We warmly thank Faezeh Khodabandehlou and Christian Maes for discussions. MP is grateful to Villa ai Colli Research Institute for hospitality.

**Funding information** PH was supported by the project INTER/FNRS/20/15074473 funded by F.R.S.-FNRS (Belgium) and FNR (Luxembourg). SDC and VL acknowledge support from IXXI, CNRS MITI, CNRS Tremplin and the ANR-18-CE30-0028-01 grant LABS.

A Construction of potentials and affinities from transition rates

Given a spanning tree $\mathcal{T}$ with an arbitrarily fixed root $x_0$, we define for $x \leftarrow y \in \overline{\mathcal{T}}$ the local affinity $A_{x \leftarrow y} := \log \frac{r_{x \leftarrow y}}{r_{y \leftarrow x}}$. For any vertex $z \in \mathcal{X}$, there is a unique path in $\mathcal{T}$ going from $z$ to $x_0$, that we denote $\mathcal{P}_z = (x_0 \leftarrow x_1, x_1 \leftarrow x_2, \dots, x_p \leftarrow z)$. The potential in $z$ is defined as the “integral” of the local affinity along this path

$$u_z = \sum_{x \leftarrow y \in \mathcal{P}_z} A_{x \leftarrow y}. \qquad (A.1)$$

This implies that the local affinity is the (discrete) “gradient” of the potential: $u_y - u_x = A_{x \leftarrow y}$ for $x \leftarrow y \in \overline{\mathcal{T}}$. This relation ceases to be true in general on a chord $x \leftarrow y = \alpha \in \mathcal{E} \setminus \mathcal{T}$ (unless the Wegscheider condition is satisfied along its associated cycle $\mathcal{C}_\alpha$). The so-called cycle affinity $a_\alpha := \log \frac{r_{x \leftarrow y}}{r_{y \leftarrow x}} - u_y + u_x$ quantifies precisely how it is broken. The potential and cycle affinity so defined verify the relations of Eqs. (9)-(10) of the main text. Note that the potential depends only on the spanning-tree rates $(r_e)_{e \in \mathcal{T}}$, while the cycle affinity depends on these rates and on the rates $r_\alpha$ of the chord (more specifically, it depends on the cycle rates $(r_e)_{e \in \mathcal{C}_\alpha}$).

B Linear independence of admissible currents

An admissible set $\mathcal{I}$ of edges can be represented as a subset of a fundamental set of chords associated to a spanning tree $\mathcal{T}$ (i.e. $\mathcal{I} \subseteq \mathcal{E} \setminus \mathcal{T}$). It is thus natural to expect that their corresponding currents $(J_i)_{i \in \mathcal{I}}$ are linearly independent when varying the transition rates of the continuous-time Markov chain. The following lemma makes this statement precise.

**Independence lemma:** Fix the rates $r_{\mathcal{E} \setminus \mathcal{I}}$ of the transitions not in the admissible set. Then, denoting by $J_i(r_\mathcal{I})$ the stationary current seen as a function of the remaining rates $r_\mathcal{I}$,

$$\exists \lambda^0, (\lambda_i)_{i \in \mathcal{I}} : \forall r_\mathcal{I}, \lambda^0 + \sum_{i \in \mathcal{I}} \lambda_i J_i(r_\mathcal{I}) = 0 \quad \Rightarrow \quad \lambda^0 = 0, \text{ and } \forall i \in \mathcal{I}, \lambda_i = 0.$$

*Proof* — Sending the rates $r_\mathcal{I}$ to 0, the stationary currents $(J_i)_{i \in \mathcal{I}}$ also go to 0, which implies $\lambda^0 = 0$. Now, fix $i \in \mathcal{I}$. This chord is associated to a cycle $\mathcal{C}_i$ which does not involve the other edges in $\mathcal{I}$. Sending $(r_i)_{i \in \mathcal{I} \setminus \{i\}}$ to 0 sends $(J_i)_{i \in \mathcal{I} \setminus \{i\}}$ to 0, and by hypothesis, we are left with

$$\forall r_i, \lambda_i J_i(r_\mathcal{I})|_{r_i \neq i=0} = 0. \qquad (B.1)$$

For any choice of $r_{\mathcal{E} \setminus \mathcal{I}}$, one can find two rates $r_i$ such that the Wegscheider condition along cycle $\mathcal{C}_i$ is not satisfied. For such a choice, one has $J_i \neq 0$. This implies that $\lambda_i = 0$, and since this reasoning can be done for every $i \in \mathcal{I}$, this ends the proof. □
Remark — This independence lemma implies that the linear-affine decomposition of Eq. (7) of a current $J_e$ on a set of admissible currents $J_I$ with coefficients depending on the rates $r_{\mathcal{E}\setminus I}$ is unique, and this for any value of $r_{\mathcal{E}\setminus I}$.

### C Proof that $1 - \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}(r_{\mathcal{E}\setminus I_1}) \lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}(r_{\mathcal{E}\setminus I_2})$ cannot cancel

We first rewrite Eqs. (22)-(23) as

$$J_{e_1} = J_{e_1}^{\mathcal{I}_2} + \lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2} J_{e_2} + \sum_{i \in \mathcal{I}} \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2} J_i, \quad (C.1)$$

$$J_{e_2} = J_{e_2}^{\mathcal{I}_1} + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} J_{e_1} + \sum_{i \in \mathcal{I}} \lambda_{e_2 \leftarrow i}^{\mathcal{I}_1} J_i, \quad (C.2)$$

where in this Appendix, for short, we omit the arguments of functions, i.e. : $J_{e_2}^{\mathcal{I}_1} = J_{e_2}^{\mathcal{I}_1}(r_{\mathcal{E}\setminus I_1})$, $J_{e_1}^{\mathcal{I}_2} = J_{e_1}^{\mathcal{I}_2}(r_{\mathcal{E}\setminus I_2})$, $\lambda_{e_2 \leftarrow i}^{\mathcal{I}_1} = \lambda_{e_2 \leftarrow i}^{\mathcal{I}_1}(r_{\mathcal{E}\setminus I_1})$, and $\lambda_{e_1 \leftarrow i}^{\mathcal{I}_2} = \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2}(r_{\mathcal{E}\setminus I_2})$. Equivalently, one has:

$$M \begin{pmatrix} J_{e_1} \\ J_{e_2} \end{pmatrix} = b, \quad \text{with} \quad M = \begin{pmatrix} 1 & -\lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2} \\ -\lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} & 1 \end{pmatrix}, \quad \text{and} \quad b = \begin{pmatrix} J_{e_1}^{\mathcal{I}_2} + \sum_{i \in \mathcal{I}} \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2} J_i \\ J_{e_2}^{\mathcal{I}_1} + \sum_{i \in \mathcal{I}} \lambda_{e_2 \leftarrow i}^{\mathcal{I}_1} J_i \end{pmatrix}, \quad (C.3)$$

which we read as an equation for $J_{e_1}$ and $J_{e_2}$. Provided $M$ is invertible, we obtain the mutual multilinearity needed in the inductive step of Sec. 3.2 (as also directly shown in the text). We now show *ad absurdum* that $M$ is always invertible.

Let us assume that $\det M = 1 - \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} \lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2} = 0$. Then necessarily both $\lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}$ and $\lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2}$ do not cancel and we have $\lambda_{e_1 \leftarrow e_2}^{\mathcal{I}_2} = 1/\lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}$. Using this in Eq. (C.1) and multiplying by $\lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}$ we arrive at

$$-J_{e_2} = \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} J_{e_1}^{\mathcal{I}_2} - \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} J_{e_1} + \sum_{i \in \mathcal{I}} \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2} J_i. \quad (C.4)$$

Adding Eq. (C.2), this yields

$$0 = J_{e_2}^{\mathcal{I}_1} + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} J_{e_1}^{\mathcal{I}_2} + \sum_{i \in \mathcal{I}} \left( \lambda_{e_2 \leftarrow i}^{\mathcal{I}_1} + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2} \right) J_i, \quad (C.5)$$

which is a linear-affine interdependence of the admissible currents $J_I$ with coefficients in $r_{\mathcal{E}\setminus I}$ (since $(\mathcal{E}\setminus I_1) \cup (\mathcal{E}\setminus I_2) = \mathcal{E}\setminus I$). From the independence lemma of App. B, we thus have

$$J_{e_2}^{\mathcal{I}_1} + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} J_{e_1}^{\mathcal{I}_2} = 0, \quad (C.6)$$

$$\forall i \in \mathcal{I}, \quad \lambda_{e_2 \leftarrow i}^{\mathcal{I}_1} + \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1} \lambda_{e_1 \leftarrow i}^{\mathcal{I}_2} = 0. \quad (C.7)$$

We now come back to the matrix equation (C.3). By hypothesis, $M$ is singular; hence, Eq. (C.3) has solutions provided $b \in \text{Im } M$. Since $\text{Im } M = (\text{Ker } M^{\top})^{\perp}$ and $\text{Ker } M^{\top} = \text{Span}(\lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}, 1)^{\top}$ this condition is equivalent to $b \cdot (\lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1}, 1) = 0$, but this is precisely guaranteed by Eqs. (C.6)-(C.7). Now, because $\text{Ker } M = \text{Span}(1, \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1})^{\top}$, this means that Eq. (C.3) has a family of solutions of the form $(J_{e_1}, J_{e_2})^{\top} = M^{\oplus}b + \gamma(1, \lambda_{e_2 \leftarrow e_1}^{\mathcal{I}_1})^{\top}$ for any real number $\gamma$, where $M^{\oplus}$ is the Moore-Penrose pseudo-inverse of $M$. This is absurd because then we have several distinct linear-affine decompositions of $J_{e_1}$ on the admissible currents $J_I$ with coefficients depending on $r_{\mathcal{E}\setminus I}$, which is excluded by the independence lemma. □
D A fluctuation-dissipation relation

Along each directed edge *i*, we define the time-integrated current $J_i$ of a trajectory of duration *t* as the number of times transition +*i* occurs minus the number of times transition −*i* occurs. This quantity is time-additive and, upon a transition *e* ∈ $\overline{\mathcal{E}}$, evolves as

$$J_i \mapsto J_i + \delta_{e,+i} - \delta_{e,-i}. \qquad (D.1)$$

Its average and second cumulant are related to the stationary mean and covariance of the current as

$$J_i = \lim_{t \to \infty} \frac{1}{t} \langle J_i \rangle, \quad c_{i_1 i_2} = \lim_{t \to \infty} \frac{1}{t} \langle J_{i_1} J_{i_2} \rangle_c. \qquad (D.2)$$

We define similarly the entropy production $\Sigma$ of a trajectory of duration $t$ as the logarithm of the ratio between the probability of this trajectory and the probability of its time reverse. It is a time-additive observable that increases at every jump: if transition $e \in \overline{\mathcal{E}}$ occurs, one has

$$\Sigma \mapsto \Sigma + \log \frac{r_{+e}}{r_{-e}}. \qquad (D.3)$$

From the definitions and from the parametrization Eqs. (9)-(10) of the rates (see App. A), it is related to the currents $J_i$'s by

$$\Sigma = \sum_{\alpha} a_{\alpha} J_{\alpha} + \text{b.t.}, \qquad (D.4)$$

where the sum runs over the chords {$\alpha$} and the boundary term b.t. depends only on the initial and final configurations of the trajectory. We recall that $a_\alpha$ is the cycle affinity associated to chord $\alpha$. We also see from the definition of $\Sigma$, that for any observable $\mathcal{O}$ depending on the trajectory, we have

$$\langle \mathcal{O}(\text{traj}) \rangle = \langle \mathcal{O}(\text{traj}^R) e^{-\Sigma(\text{traj})+\text{b.t.}} \rangle, \qquad (D.5)$$

where traj$^R$ denotes the time reversed of trajectory traj. For instance $J_i(\text{traj}^R) = -J_i(\text{traj})$.

The cumulant generating function of the currents, defined as

$$\psi(s) = \lim_{t \to \infty} \frac{1}{t} \langle e^{\sum_{e \in \mathcal{E}} s_e J_e} \rangle, \qquad (D.6)$$

is such that $J_e = \partial_{s_e} \psi|_{s=0}$ and $c_{e_1 e_2} = \partial_{s_{e_1}} \partial_{s_{e_2}} \psi|_{s=0}$. From Eqs. (D.4)-(D.5), it satisfies

$$\psi(s^{\text{cochord}}, s^{\text{chord}}, a) = \psi(-s^{\text{cochord}}, -a - s^{\text{chord}}, a), \qquad (D.7)$$

where we split the full vector $s$ into its chord contribution $s^{\text{chord}} = (s_\alpha)$ and its cochord one, and $a = (a_\alpha)$ represent the affinity vector of the chords. Eq. (D.7) constitutes a fluctuation theorem that (slightly) generalizes the one derived in Ref. [51]. Differentiating it w.r.t. $s_e$ and $a_\alpha$, and then sending $s$ and $a$ to zero, we arrive at

$$\left. \frac{\partial^2 \psi(s, a)}{\partial s_e \partial a_\alpha} \right|_{\substack{s=0 \\ a=0}} = \frac{1}{2} \left. \frac{\partial^2 \psi(s, a)}{\partial s_e \partial s_\alpha} \right|_{\substack{s=0 \\ a=0}}. \qquad (D.8)$$

From Eq. (D.2), this implies that, at equilibrium ($a = 0$, denoted by a superscript “eq”), the response of any stationary current $J_e$ to a change in the affinity $a_\alpha$ of a chord is proportional to a current covariance:

$$\frac{\partial J_e}{\partial a_\alpha}^{\text{eq}} = \frac{1}{2} c_{e\alpha}^{\text{eq}}. \qquad (D.9)$$

This constitutes a form of fluctuation-dissipation relation at equilibrium. Because input currents can be taken as a subset of chords, this proves Eq. (53) of the main text.
References

[1] H. B. G. Casimir, *On Onsager's principle of microscopic reversibility*, Rev. Mod. Phys. **17**, 343 (1945), doi:10.1103/RevModPhys.17.343.

[2] H. B. Callen and T. A. Welton, *Irreversibility and generalized noise*, Phys. Rev. **83**, 34 (1951), doi:10.1103/PhysRev.83.34.

[3] R. Kubo, *Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems*, J. Phys. Soc. Jpn. **12**, 570 (1957), doi:10.1143/JPSJ.12.570.

[4] U. Marconi, A. Puglisi, L. Rondoni and A. Vulpiani, *Fluctuation-dissipation: Response theory in statistical physics*, Phys. Rep. **461**, 111 (2008), doi:10.1016/j.physrep.2008.02.002 [preprint doi:10.48550/arXiv.0803.0719].

[5] J. Kurchan, *Fluctuation theorem for stochastic dynamics*, J. Phys. A: Math. Gen. **31**, 3719 (1998), doi:10.1088/0305-4470/31/16/003.

[6] J. L. Lebowitz and H. Spohn, *A Gallavotti-Cohen-type symmetry in the large deviation functional for stochastic dynamics*, J. Stat. Phys. **95**, 333 (1999), doi:10.1023/A:1004589714161.

[7] D. Andrieux and P. Gaspard, *Fluctuation theorem for currents and Schnakenberg network theory*, J. Stat. Phys. **127**, 107 (2007), doi:10.1007/s10955-006-9233-5.

[8] M. Esposito and C. Van den Broeck, *Three detailed fluctuation theorems*, Phys. Rev. Lett. **104**, 090601 (2010), doi:10.1103/PhysRevLett.104.090601 [preprint doi:10.48550/arXiv.0911.2666].

[9] U. Seifert, *Entropy production along a stochastic trajectory and an integral fluctuation theorem*, Phys. Rev. Lett. **95**, 040602 (2005), doi:10.1103/PhysRevLett.95.040602 [preprint doi:10.48550/arXiv.cond-mat/0503686].

[10] A. C. Barato and U. Seifert, *Thermodynamic uncertainty relation for biomolecular processes*, Phys. Rev. Lett. **114**, 158101 (2015), doi:10.1103/PhysRevLett.114.158101 [preprint doi:10.48550/arXiv.1502.05944].

[11] T. R. Gingrich, J. M. Horowitz, N. Perunov and J. L. England, *Dissipation bounds all steady-state current fluctuations*, Phys. Rev. Lett. **116**, 120601 (2016), doi:10.1103/PhysRevLett.116.120601 [preprint doi:10.48550/arXiv.1512.02212].

[12] G. S. Agarwal, *Fluctuation-dissipation theorems for systems in non-thermal equilibrium and applications*, Z. Phys. A: Hadrons nucl. **252**, 25 (1972), doi:10.1007/BF01391621.

[13] L. F. Cugliandolo, D. S. Dean and J. Kurchan, *Fluctuation-dissipation theorems and entropy production in relaxational systems*, Phys. Rev. Lett. **79**, 2168 (1997), doi:10.1103/PhysRevLett.79.2168 [preprint doi:10.48550/arXiv.cond-mat/9705002].

[14] J. Prost, J.-F. Joanny and J. M. R. Parrondo, *Generalized fluctuation-dissipation theorem for steady-state systems*, Phys. Rev. Lett. **103**, 090601 (2009), doi:10.1103/PhysRevLett.103.090601.

[15] U. Seifert and T. Speck, *Fluctuation-dissipation theorem in nonequilibrium steady states*, Europhys. Lett. **89**, 10007 (2010), doi:10.1209/0295-5075/89/10007.
[16] B. Altaner, M. Polettini and M. Esposito, *Fluctuation-dissipation relations far from equilibrium*, Phys. Rev. Lett. **117**, 180601 (2016), doi:10.1103/PhysRevLett.117.180601 [preprint doi:10.48550/arXiv:1604.08832].

[17] J. A. Owen, T. R. Gingrich and J. M. Horowitz, *Universal thermodynamic bounds on nonequilibrium response with biochemical applications*, Phys. Rev. X **10**, 011066 (2020), doi:10.1103/PhysRevX.10.011066 [preprint doi:10.48550/arXiv:1905.07449].

[18] A. Dechant and S.-i. Sasa, *Fluctuation-response inequality out of equilibrium*, Proc. Natl. Acad. Sci. **117**, 6430 (2020), doi:10.1073/pnas.1918386117.

[19] G. Falasco, M. Esposito and J.-C. Delvenne, *Beyond thermodynamic uncertainty relations: Nonlinear response, error-dissipation trade-offs, and speed limits*, J. Phys. A: Math. Theor. **55**, 124002 (2022), doi:10.1088/1751-8121/ac52e2.

[20] G. Fernandes Martins and J. M. Horowitz, *Topologically constrained fluctuations and thermodynamics regulate nonequilibrium response*, Phys. Rev. E **108**, 044113 (2023), doi:10.1103/PhysRevE.108.044113 [preprint doi:10.48550/arXiv:2305.19348].

[21] T. Aslyamov and M. Esposito, *General theory of static response for Markov jump processes*, Phys. Rev. Lett. **133**, 107103 (2024), doi:10.1103/PhysRevLett.133.107103.

[22] T. Aslyamov and M. Esposito, *Nonequilibrium response for Markov jump processes: Exact results and tight bounds*, Phys. Rev. Lett. **132**, 037101 (2024), doi:10.1103/PhysRevLett.132.037101.

[23] K. Ptaszyński, T. Aslyamov and M. Esposito, *Dissipation bounds precision of current response to kinetic perturbations*, Phys. Rev. Lett. **133**, 227101 (2024), doi:10.1103/PhysRevLett.133.227101 [preprint doi:10.48550/arXiv:2406.08361].

[24] J. Zheng and Z. Lu, *Universal response inequalities beyond steady states via trajectory information geometry*, Phys. Rev. E **112**, L012103 (2025), doi:10.1103/scg2-qkxv [preprint doi:10.48550/arXiv:2403.10952].

[25] K. Liu and J. Gu, *Dynamical activity universally bounds precision of response in Markovian nonequilibrium systems*, Commun. Phys. **8**, 62 (2025), doi:10.1038/s42005-025-01982-w [preprint doi:10.48550/arXiv:2410.20800].

[26] H.-M. Chun and J. M. Horowitz, *Trade-offs between number fluctuations and response in nonequilibrium chemical reaction networks*, J. Chem. Phys. **158**, 174115 (2023), doi:10.1063/5.0148662 [preprint doi:10.48550/arXiv:2304.04961].

[27] E. Kwon, H.-M. Chun, H. Park and J. S. Lee, *Fluctuation-response inequalities for kinetic and entropic perturbations*, Phys. Rev. Lett. **135**, 097101 (2025), doi:10.1103/h45s-4118 [preprint doi:10.48550/arXiv:2411.18108].

[28] R. Bao and S. Liang, *Nonlinear response identities and bounds for nonequilibrium steady states*, (arXiv preprint) doi:10.48550/arXiv:2412.19602.

[29] T. Van Vu, *Fundamental bounds on precision and response for quantum trajectory observables*, PRX Quantum **6**, 010343 (2025), doi:10.1103/PRXQuantum.6.010343 [preprint doi:10.48550/arXiv:2411.19546].

[30] T. Aslyamov, K. Ptaszyński and M. Esposito, *Nonequilibrium fluctuation-response relations: From identities to bounds*, Phys. Rev. Lett. **134**, 157101 (2025), doi:10.1103/PhysRevLett.134.157101 [preprint doi:10.48550/arXiv:2410.17140].
[31] K. Ptaszynski, T. Aslyamov and M. Esposito, *Nonequilibrium fluctuation-response relations for state observables*, (arXiv preprint) doi:10.48550/arXiv.2412.10233.

[32] J. Zheng and Z. Lu, *Unified linear fluctuation-response theory arbitrarily far from equilibrium*, (arXiv preprint) doi:10.48550/arXiv.2501.01050.

[33] P. E. Harunari, S. Dal Cengio, V. Lecomte and M. Polettini, *Mutual linear-ity of nonequilibrium network currents*, Phys. Rev. Lett. **133**, 047401 (2024), doi:10.1103/PhysRevLett.133.047401.

[34] M. Polettini, P. E. Harunari, S. Dal Cengio and V. Lecomte, *Coplanarity of rooted spanning-tree vectors*, (arXiv preprint) doi:10.48550/arXiv.2407.16093.

[35] T. L. Hill, *Studies in irreversible thermodynamics IV. Diagrammatic representation of steady state fluxes for unimolecular systems*, J. Theor. Biol. **10**, 442 (1966), doi:10.1016/0022-5193(66)90137-8.

[36] J. Schnakenberg, *Network theory of microscopic and macroscopic behavior of master equation systems*, Rev. Mod. Phys. **48**, 571 (1976), doi:10.1103/RevModPhys.48.571.

[37] E. C. Kirby, R. B. Mallion, P. Pollak and P. J. Skrzyński, *What Kirchhoff actually did concerning spanning trees in electrical networks and its relationship to modern graph-theoretical work*, Croat. Chem. Acta **89**, 403 (2016), doi:10.5562/cca2995.

[38] S. Dal Cengio, V. Lecomte and M. Polettini, *Geometry of nonequilibrium reaction networks*, Phys. Rev. X **13**, 021040 (2023), doi:10.1103/PhysRevX.13.021040.

[39] M. Polettini and M. Esposito, *Irreversible thermodynamics of open chemical networks. I. Emergent cycles and broken conservation laws*, J. Chem. Phys. **141**, 024117 (2014), doi:10.1063/1.4886396 [preprint doi:10.48550/arXiv.1404.1181].

[40] P. Pietzonka, A. C. Barato and U. Seifert, *Universal bounds on current fluctuations*, Phys. Rev. E **93**, 052145 (2016), doi:10.1103/PhysRevE.93.052145 [preprint doi:10.48550/arXiv.1512.01221].

[41] M. Polettini and M. Esposito, *Effective thermodynamics for a marginal observer*, Phys. Rev. Lett. **119**, 240601 (2017), doi:10.1103/PhysRevLett.119.240601 [preprint doi:10.48550/arXiv.1703.05715].

[42] N. Freitas, J.-C. Delvenne and M. Esposito, *Stochastic thermodynamics of nonlinear electronic circuits: A realistic framework for computing around kT*, Phys. Rev. X **11**, 031064 (2021), doi:10.1103/PhysRevX.11.031064 [preprint doi:10.48550/arXiv.2008.10578].

[43] E. Ilker, Ö. Güngör, B. Kuznets-Speck, J. Chiel, S. Deffner and M. Hinczewski, *Shortcuts in stochastic systems and control of biophysical processes*, Phys. Rev. X **12**, 021048 (2022), doi:10.1103/PhysRevX.12.021048 [preprint doi:10.1101/2021.06.13.448255].

[44] N. Ohga, S. Ito and A. Kolchinsky, *Thermodynamic bound on the asymmetry of cross-correlations*, Phys. Rev. Lett. **131**, 077101 (2023), doi:10.1103/PhysRevLett.131.077101 [preprint doi:10.48550/arXiv.2303.13116].

[45] F. Avanzini et al., *Methods and conversations in (post)modern thermodynamics*, SciPost Phys. Lect. Notes 80 (2024), doi:10.21468/SciPostPhysLectNotes.80 [preprint doi:10.48550/arXiv.2311.01250].
[46] F. Khodabandehlou, C. Maes and K. Netočný, *Affine relationships between steady currents*, J. Phys. A: Math. Theor. **58**, 155002 (2025), doi:10.1088/1751-8121/adc8ea [preprint doi:10.48550/arXiv.2412.05019].

[47] C. Monthus, *Markov generators as non-Hermitian supersymmetric quantum Hamiltonians: Spectral properties via bi-orthogonal basis and singular value decompositions*, J. Stat. Mech.: Theory Exp. 083207 (2024), doi:10.1088/1742-5468/ad613a.

[48] P. Chebotarev and R. Agaev, *Forest matrices around the Laplacian matrix*, Linear Algebra Appl. **356**, 253 (2002), doi:10.1016/S0024-3795(02)00388-9 [preprint doi:10.48550/arXiv.math/0508178].

[49] F. Khodabandehlou, C. Maes and K. Netočný, *Trees and forests for nonequilibrium purposes: An introduction to graphical representations*, J. Stat. Phys. **189**, 41 (2022), doi:10.1007/s10955-022-03003-4.

[50] M. Polettini and M. Esposito, *Effective fluctuation and response theory*, J. Stat. Phys. **176**, 94 (2019), doi:10.1007/s10955-019-02291-7.

[51] D. Andrieux and P. Gaspard, *Fluctuation theorem and Onsager reciprocity relations*, J. Chem. Phys. **121**, 6167 (2004), doi:10.1063/1.1782391.

[52] J. W. Rocks, H. Ronellenfitsch, A. J. Liu, S. R. Nagel and E. Katifori, *Limits of multifunctionality in tunable networks*, Proc. Natl. Acad. Sci. **116**, 2506 (2019), doi:10.1073/pnas.1806790116 [preprint doi:10.48550/arXiv.1805.00504].

[53] M. Stern and A. Murugan, *Learning without neurons in physical systems*, Annu. Rev. Condens. Matter Phys. **14**, 417 (2023), doi:10.1146/annurev-conmatphys-040821-113439.

[54] M. Guzman, F. Martins, M. Stern and A. J. Liu, *Microscopic imprints of learned solutions in tunable networks*, Phys. Rev. X **15**, 031056 (2025), doi:10.1103/f2hb-c9s1 [preprint doi:10.48550/arXiv.2412.19356].

[55] A. Blokhuis, D. Lacoste and P. Nghe, *Universal motifs and the diversity of autocatalytic systems*, Proc. Natl. Acad. Sci. **117**, 25230 (2020), doi:10.1073/pnas.2013527117.